简介：这份MATLAB脚本面向声学测试、环境噪声监测与工业设备噪声评估方向的工程师及研究人员，用于对时域噪声信号执行标准A计权处理并完成1/3倍频程频谱划分与声压级计算。脚本内置ISO 266与IEC 61672推荐的A计权滤波器系数，支持自定义采样率与频带范围，输入为列向量格式的原始时域信号，输出各1/3倍频程中心频率对应的A加权声压级，便于后续绘图或导出，无需额外工具箱，兼容MATLAB R2015b及以上版本。资源包共8个文件，约250KB，包含核心m脚本、Python对照实现、A计权曲线与NR曲线示意图、依赖说明及配置文件等，结构紧凑、开箱即用。已有106人学习下载，适合需要快速搭建合规性声学分析流程、对照标准验证算法或作为教学演示参考的读者。
1. 为什么你的声压级数据总被质疑：A计权与1/3倍频程到底在解决什么
做噪声测试的工程师大概率遇到过这种场景：采集回来一段时域声压信号，直接算总声压级，结果和第三方检测报告对不上，差了 3 到 8 dB(A)。问题往往不在设备，而在后处理——没有做 A 计权，也没有做 1/3 倍频程分解。A 计权模拟人耳对低频不敏感、对 2 到 5 kHz 敏感的听感特性，1/3 倍频程则把宽频噪声拆成中心频率从 20 Hz 到 20 kHz 的 31 个频带，让每个频带的贡献一目了然。MATLAB 实现 A 计权 1/3 倍频程声压级分析脚本，本质就是把这两步标准化、可复现地串起来：读入时域波形，加 A 计权滤波器，做频谱分解，输出各频带声压级和总 A 计权级。适合做产品噪声测试、环境噪声评估、电机与风机异响定位的从业者，也适合需要批量处理多通道数据的实验室场景。
2. A计权滤波器怎么在MATLAB里落地：从模拟传递函数到数字IIR
2.1 A计权曲线的数学来源与双线性变换
A 计权的定义来自 IEC 61672，它由三个极点两个零点的模拟传递函数描述。核心思路是：先把模拟域传递函数写成零极点形式，再用双线性变换映射到数字域，得到 IIR 滤波器系数。模拟传递函数的标准形式包含四个部分：两个位于原点的零点、两个高频极点、一个低频极点，以及一个与参考频率相关的增益因子。直接手推容易出错，常见做法是用bilinear函数配合预先算好的零极点向量。
我一般会先把模拟零极点写清楚：
% A计权模拟传递函数零极点（基于IEC 61672） f1 = 20.598997; % 低频极点 Hz f2 = 107.65265; % 高频零点 Hz f3 = 737.86223; % 高频极点 Hz f4 = 12194.217; % 高频极点 Hz % 模拟极点（rad/s） p = -2*pi*[f1; f3; f4]; % 模拟零点 z = -2*pi*[f2; f2]; % 注意：两个零点在f2处 % 增益因子，使1kHz处归一化为0 dB k = (2*pi*f4)^2 / (2*pi*f2)^2;这段代码里f1到f4是标准常数，不要随意改。p是三个极点，z是两个零点，k是增益。接下来用双线性变换：
fs = 48000; % 采样率，必须与你的数据一致 [b,a] = bilinear(z, p, k, fs);bilinear会自动完成预畸变，但要注意：采样率低于 40 kHz 时，接近奈奎斯特频率的频带会有偏差，建议采样率至少 48 kHz。如果原始数据是 44.1 kHz，也能用，但 16 kHz 以上频带的 A 计权值会略偏。
2.2 用freqz验证滤波器频响是否贴合标准
滤波器系数算出来后，别急着往数据上套。先画频响，和标准 A 计权曲线对比。这一步是很多人的血泪经验：系数算错一个符号，低频段能差 10 dB。
f = logspace(1, log10(fs/2), 2000); [H, ~] = freqz(b, a, f, fs); A_db = 20*log10(abs(H)); % 标准A计权参考值（部分频点） f_ref = [20 100 1000 10000 20000]; A_ref = [-50.5 -19.1 0 -2.5 -9.3]; semilogx(f, A_db, 'b', 'LineWidth', 1.5); hold on; semilogx(f_ref, A_ref, 'ro', 'MarkerSize', 8); grid on; xlabel('频率 (Hz)'); ylabel('增益 (dB)'); legend('MATLAB计算','标准参考');跑完这张图，1 kHz 处应该在 0 dB 附近，20 Hz 处约 -50 dB，10 kHz 处约 -2.5 dB。如果 1 kHz 偏了超过 0.1 dB，检查k的计算；如果高频段翘起来，检查f4是否写错。这一步过了，后面才有意义。
2.3 时域滤波的边界处理：filtfilt还是filter
直接filter(b,a,x)会引入相位失真，而且首尾有瞬态。声压级分析关心幅值，不关心相位，所以推荐filtfilt做零相位滤波。但filtfilt会延长信号，首尾各补一段，对于短信号要小心。
% 假设x是时域声压，单位Pa x = x(:); % 强制列向量 x_A = filtfilt(b, a, x); % 零相位A计权滤波如果信号长度小于 3 倍滤波器阶数，filtfilt会报错或结果不可靠。常见做法是信号至少 1 秒以上。另外，filtfilt默认用反射填充，对于非平稳噪声（比如撞击声）首尾可能失真，可以改用filter然后手动去掉前 0.1 秒。
注意：A 计权滤波后的信号单位仍是 Pa，但已经是“A 计权声压”，后续算级时不要再乘 A 计权修正值。
3. 1/3倍频程分解：从FFT到频带声压级的完整链路
3.1 中心频率与频带上下限的精确计算
1/3 倍频程的 31 个中心频率不是随便取的，标准规定从 20 Hz 到 20 kHz，每个中心频率是前一个的 2^(1/3) 倍。但实际标准（IEC 61260）给出了精确值，并且规定了每个频带的上下限。常见做法是直接查表，但用公式生成更灵活：
% 1/3倍频程中心频率（20 Hz到20 kHz，共31个） f_center = 20 * 2.^((0:30)/3); % 修正到标准标称值 f_nominal = [20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 ... 500 630 800 1000 1250 1600 2000 2500 3150 4000 5000 ... 6300 8000 10000 12500 16000 20000];f_center是精确计算值，f_nominal是报告里常用的标称值。两者在 1 kHz 以下略有差异，报告里一般用标称值。频带上下限按f_lower = f_center / 2^(1/6)，f_upper = f_center * 2^(1/6)计算。
3.2 用FFT做频带能量累加的正确姿势
最直接的方法：对 A 计权后的信号做 FFT，得到功率谱，然后在每个频带范围内累加功率，再转成声压级。这里有个坑：FFT 的频率分辨率必须足够细，否则低频段一个频带里只有一两个谱线，误差很大。
N = length(x_A); X = fft(x_A .* hanning(N)); % 加汉宁窗减少泄漏 P = (abs(X(1:N/2+1)).^2) / (fs * sum(hanning(N).^2)); % 功率谱密度 f_axis = (0:N/2)' * fs / N; % 每个频带累加 Lp_band = zeros(31,1); for k = 1:31 fl = f_center(k) / 2^(1/6); fu = f_center(k) * 2^(1/6); idx = f_axis >= fl & f_axis < fu; if any(idx) p_rms = sqrt(sum(P(idx)) * (f_axis(2)-f_axis(1))); % 频带声压有效值 Lp_band(k) = 20*log10(p_rms / 2e-5); % 参考声压20微帕 else Lp_band(k) = -Inf; % 该频带无谱线 end end关键参数：hanning(N)窗函数，sum(hanning(N).^2)是窗函数功率归一化因子，少了这个因子，声压级会整体偏高。2e-5是参考声压，单位 Pa。如果数据已经是 A 计权后的，这里算出来的就是各频带 A 计权声压级。
3.3 总A计权声压级的两种合成方式与验证
总 A 计权声压级有两种算法：一是直接对 A 计权后的时域信号算总 RMS，二是把 31 个频带声压级按能量相加。两者应该一致，差异超过 0.5 dB 说明频带分解有遗漏。
% 方式一：时域总RMS Lp_total_time = 20*log10(rms(x_A) / 2e-5); % 方式二：频带能量合成 Lp_total_band = 10*log10(sum(10.^(Lp_band/10))); fprintf('时域总级: %.2f dB(A)\n', Lp_total_time); fprintf('频带合成: %.2f dB(A)\n', Lp_total_band);如果两者差得多，检查 FFT 是否漏了 20 Hz 以下或 20 kHz 以上的能量，或者窗函数归一化是否写错。我一般会把这个差值打印出来，超过 0.3 dB 就回头查。
4. 脚本工程化：批量处理、参数配置与结果导出
4.1 把采样率、参考声压、频带表做成配置结构体
一次性脚本跑通不难，难的是换一组数据就要改十几处。我习惯把参数收进一个结构体：
cfg.fs = 48000; cfg.p_ref = 2e-5; cfg.f_center = 20 * 2.^((0:30)/3); cfg.f_nominal = [20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 ... 500 630 800 1000 1250 1600 2000 2500 3150 4000 5000 ... 6300 8000 10000 12500 16000 20000]; cfg.window = @hanning; cfg.filtfilt = true;这样换采样率只改cfg.fs，换参考声压只改cfg.p_ref。频带表也统一从cfg取，避免函数里硬编码。
4.2 多通道批量处理与进度反馈
实验室数据经常是 8 通道或 16 通道同步采集。用循环处理，每个通道独立算，最后拼成矩阵：
n_ch = size(data, 2); Lp_all = zeros(31, n_ch); Lp_total = zeros(1, n_ch); for ch = 1:n_ch x = data(:, ch); x_A = filtfilt(b, a, x); [Lp_all(:,ch), Lp_total(ch)] = band_analysis(x_A, cfg); fprintf('通道 %d/%d 完成，总级 %.2f dB(A)\n', ch, n_ch, Lp_total(ch)); endband_analysis是封装好的函数，输入 A 计权信号和配置，输出频带级和总级。进度打印看着简单，但批量跑 16 通道时能让你知道是不是卡在某个通道。
4.3 导出CSV与绘图：让报告直接可用
结果导出成 CSV，列是频带中心频率，行是通道。再画一张 1/3 倍频程柱状图，报告里直接贴。
T = array2table([cfg.f_nominal' Lp_all], ... 'VariableNames', [{'CenterFreq_Hz'}, ... arrayfun(@(x) sprintf('Ch%d_dBA',x), 1:n_ch, 'UniformOutput', false)]); writetable(T, 'octave_results.csv'); bar(cfg.f_nominal, Lp_all(:,1)); set(gca, 'XScale', 'log'); xlabel('中心频率 (Hz)'); ylabel('声压级 dB(A)'); grid on;array2table的列名用sprintf生成，避免手动写 16 个名字。bar配合XScale对数轴，出来的图就是标准 1/3 倍频程谱图。
5. 避坑与排查：A计权1/3倍频程分析中最容易翻车的5个点
5.1 现象：低频段声压级异常高，20 Hz处比100 Hz还大
原因：A 计权滤波器系数算错，或者根本没加 A 计权，直接对原始信号做频带分解。原始信号里低频能量通常很强，不做 A 计权就会看到 20 Hz 处很高。
解决：先画freqz频响，确认 20 Hz 处增益约 -50 dB。如果频响不对，检查bilinear的零极点输入顺序，MATLAB 的bilinear要求零点列向量、极点列向量、增益、采样率。
5.2 现象：总级和第三方报告差 2 到 3 dB
原因：参考声压用错。空气声常用 20 微帕，但水声用 1 微帕。如果数据来自水声或特殊传感器，参考声压不同。
解决：确认cfg.p_ref。另外检查窗函数归一化因子，sum(win.^2)漏掉会导致总级偏高约 1.5 dB。
5.3 现象：频带合成总级比时域总级低很多
原因：FFT 频率分辨率不够，低频段一个频带里没有谱线，被-Inf丢掉了。比如采样率 48 kHz，N=4800，分辨率 10 Hz，20 Hz 频带（下限 17.8 Hz，上限 22.4 Hz）可能只有一条谱线。
解决：增加 FFT 点数，或者对低频段单独做更细的 FFT。常见做法是 N 取 2 的整数次幂，至少 16384，保证分辨率小于 3 Hz。
5.4 现象：filtfilt报错“信号长度必须大于滤波器阶数的3倍”
原因：信号太短，或者滤波器阶数太高。A 计权滤波器通常 4 到 6 阶，但双线性变换后可能到 8 阶以上。
解决：检查信号长度，至少 1 秒。如果确实短，改用filter然后手动去掉前 0.05 秒。或者用designfilt设计更低阶的滤波器。
5.5 现象：多通道数据某一通道结果全为-Inf
原因：该通道数据全是零或全是 NaN。采集时可能某个通道没接传感器。
解决：在循环里加判断if all(x==0) || any(isnan(x))，跳过并记录。不要直接算，否则-Inf会污染后续合成。
6. 进阶技巧：用pwelch替代FFT做频带分析，以及A计权验证的快速方法
前面用 FFT 加窗做频带分解，适合稳态噪声。如果信号是非稳态的，比如风噪随时间变化，用pwelch做功率谱估计更稳。pwelch自带分段平均，方差更小。
[Pxx, f_axis] = pwelch(x_A, hanning(4096), 2048, 4096, cfg.fs); % Pxx单位是Pa^2/Hz，后续频带累加同前hanning(4096)是段长，2048是重叠，4096是 FFT 点数。段长决定频率分辨率，4096 点在 48 kHz 下分辨率约 11.7 Hz，对 20 Hz 频带仍偏粗，可以加大到 8192。pwelch的Pxx已经是功率谱密度，频带累加时乘分辨率即可。
另一个技巧：验证 A 计权滤波器是否准确，不用每次都画频响。用 1 kHz 正弦波做输入，A 计权后幅值应该几乎不变（0 dB）；用 100 Hz 正弦波，A 计权后幅值应衰减约 19.1 dB。写个快速自检函数：
function ok = check_A_weight(b, a, fs) t = (0:fs-1)'/fs; x1k = sin(2*pi*1000*t); x100 = sin(2*pi*100*t); y1k = filtfilt(b, a, x1k); y100 = filtfilt(b, a, x100); g1k = 20*log10(rms(y1k)/rms(x1k)); g100 = 20*log10(rms(y100)/rms(x100)); ok = abs(g1k) < 0.1 && abs(g100 + 19.1) < 0.5; fprintf('1kHz增益 %.2f dB, 100Hz增益 %.2f dB\n', g1k, g100); end这个自检跑一次不到一秒，但能挡住大部分系数错误。我现在的习惯是：每次换采样率或换 MATLAB 版本，先跑这个自检，过了再处理正式数据。声压级分析这行，数据量一大，一个系数错误能让你白跑一整天，后悔药没处买。希望帮到你。
本文还有配套的精品资源，点击获取
