MATLAB实现三分之一倍频程计算：频带划分、滤波器组与FFT掩码

简介面向信号处理与音频分析初学者的MATLAB三分之一倍频程计算源码包可应用于环境噪声评估、音响系统设计等常见场景重点解决倍频程滤波器组构建、频带能量统计与频谱图绘制的可落地实现问题是学习1/3倍频程概念的实用参考。压缩包仅4KB体积包含3个.m脚本分别覆盖时域滤波处理、频域功率计算及示例演示结构简明紧凑适合对照学习或直接改造复用。目前已有3000余人学习下载具备较好的实用性与教学参考价值。资源精选了可直接运行的MATLAB程序完整展示了从滤波器参数设置、filter滤波、FFT功率谱计算到plot绘图输出的全流程同时结合描述中的方法讲解可帮助初学者理解1/3倍频程划分中2^(1/3)倍频率间隔、各频带能量分布等关键概念并可通过替换输入数据快速迁移到噪声监测、音响频响修正等实际项目中对需要开展音频频段识别与声学分析的用户尤为适用。1. 为什么三分之一倍频程计算我建议你亲手在 MATLAB 里写一遍打开一个环境噪声采集文件三分之一倍频程几乎总是第一个要画的频谱。它比 FFT 原始频率分辨率更贴近人耳也比倍频程 1/1 Octave 更能看出空调噪声、变压器嗡嗡声和交通噪声的差异。通常大家第一步就用poctave出图但真正要换成自定义中心频率、做主成分分析或者批量处理多通道数据时内置函数反而不自由。项目里给的one_three_time_Tim.m和one_three_time_Feq.m对应两条不同路线一路走时域滤波器组一路走 FFT 频带积分。把这两条路都撸一遍才能看懂参数为什么那么设也才能在算错时知道往哪里查。这适合信号处理/声学从业者也适合课程设计里接了 1/3 倍频壁纸任务的学生。前提是你对 MATLAB 矩阵操作和fft至少不陌生。2. 三分之一倍频程的频带划分中心频率、上下限与掩码矩阵生成在动手算任何 1/3 倍频之前先把频带边界搞清楚。三分之一倍频程的意思是每个频带的上限频率是下限频率的 2^(1/3) 倍约为 1.26 倍而中心频率是上下限的几何平均因此上下限分别约为中心频率的 0.89 倍和 1.12 倍。很多人在这里直接把中心频率上下取 ±10% 或 ±1/3 倍出来的结果和标准仪器对不上。2.1 用 10^(1/10) 还是 2^(1/3) 来定义中心频率国际标准IEC 61260 / ISO 266更常用基准频率 1000 Hz 和整数 n中心频率定义为f_c 1000 * 10^(n/10)。因为 10^(1/10) 1.2589已经很接近 2^(1/3)1.2599但在低频段和高频段会累积出一个可测量的差异。做声学评估我建议统一的中心频率序列就用基准频率和十进制指数n中心频率/Hz下限/Hz上限/Hz-5316.23281.84354.81-4398.11354.81446.68-3501.19446.68562.34-2630.96562.34707.95-1794.33707.95891.2501000.00891.251122.0211258.931122.021412.54这里的下限是fc * 10^(-1/20)上限是fc * 10^(1/20)。之所以用 1/20 而不是 1/6是为了让带宽正好为 10^(1/10) 倍标准仪器选一档归一化FFT 掩码直接做在等对数间隔上代码写起来也不容易出错。2.2 生成中心频率与频带掩码的 MATLAB 基础函数我一般先把“频率范围”转成“整数指数 n”再生成中心频率function fc oct3_fc(f_low, f_high) % 生成 1/3 倍频程中心频率序列 % 输入最低和最高分析频率输出列向量 fc f_anchor 1000; n_low ceil(10 * log10(f_low / f_anchor)); n_high floor(10 * log10(f_high / f_anchor)); n (n_low:n_high).; fc f_anchor * 10 .^ (n / 10); end这个函数的关键是log10(f_low/f_anchor)之后再乘 10因为每十倍频程包含 10 个 1/3 倍频程频带。ceil和floor保证你传入 20 Hz、20 kHz 时得到的频带不会越过边界。如果你用 2^(1/3) 体系把 10 换成3*log10(2)也可以但表格、仪器和报告都按十进制指数更顺。生成了中心频率之后掩码矩阵就建立在离散频率轴上。假设 FFT 频率轴是f_axis一个mask(k, j)表示第 k 个频带是否包含第 j 条频率线function M oct3_mask(fc, f_axis) % fc: 中心频率向量, f_axis: 频率轴行向量 % 返回逻辑矩阵 M, size [numel(fc), numel(f_axis)] f_low fc * 10^(-1/20); f_high fc * 10^(1/20); M (f_axis f_low) (f_axis f_high); end注意这里没有考虑 FFT 谱线落在频带边界时应该归到哪一侧实际应用时边界谱线能量会同时被两侧边界截断影响不大但如果你做的是高精度声级计标定建议使用三角窗或能量补插来处理落在边界上的谱线后面第 4 章会展开。这一层的选型要点是掩码矩阵更适合宽频带 FFT 批处理一次矩阵乘法就能拿到所有频带能量而时域滤波器组虽然在非稳态信号上有优势但状态变量多、批处理也不够直观。两种方法的计算量在这儿开始分化。3. one_three_time_Tim.m 的时域滤波路线滤波器组设计与能量计算第二个目录里的one_three_time_Tim.m命名已经点出这是“time domain”的算法。它处理的对象是完整时域波形核心思路是将原始信号通过一组带通滤波器得到每个 1/3 倍频程频带内的时域分量再分别计算能量。这条路在瞬态冲击、语音和机械振动分析里更稳因为滤波后的时间波形可以继续做包络、峰值因子、短时能量等分析。3.1 Butterworth 滤波器阶数与归一化参数换算MATLAB 里最省事的滤波器是butter但设计带通滤波器时参数要换算清楚。采样率Fs 44100时中心频率fc 1000上下限分别为 891.25 Hz 和 1122.02 Hz归一化边界要除以奈奎斯特频率Fs/2Fs 44100; fc 1000; f_low fc / 10^(1/20); f_high fc * 10^(1/20); Wn [f_low f_high] / (Fs/2); [b, a] butter(4, Wn, bandpass); y filtfilt(b, a, x); % 零相位滤波避免波形偏移参数说明butter的第一个参数 4 表示 4 阶带通整体等效为 8 阶阻带衰减约 48 dB/oct。对 1/3 倍频程来说4 阶 Butterworth 的带外衰减基本够用但相邻频带之间其实没有硬性正交要求噪声测试里更在意总衰减量。如果你发现相邻频带串扰明显把阶数提高到 6 阶通常能改善但低频段滤波器会变得更不稳定filtfilt也会在边界产生较长瞬态。3.2 filtfilt 与 filter 的选择和边界效应从代码逻辑上filter是因果滤波filtfilt做了正向和反向两次滤波可以做到零相位波形起始点和结束点的相位畸变更小。filtfilt的代价是边界会引入额外瞬态等效将信号向两端各自延长了约 3 倍滤波器阶数的群延迟。处理平稳噪声时用filtfilt如果信号本身很短、或者你只关心稳态段直接用filter并在前面加一段用于等待滤波器稳定的前导信号更合适。for k 1:numel(fc) Wn [fc(k)/10^(1/20), fc(k)*10^(1/20)] / (Fs/2); if Wn(1) 0 || Wn(2) 1 continue; % 超出分析范围的频带要跳过 end [b, a] butter(4, Wn, bandpass); y_band filtfilt(b, a, x); energy(k) sum(y_band.^2) / length(y_band); end这里的关键点是用sum(y.^2)/N得到平均功率而不是sum(y.^2)否则波形长度不同时无法比较。如果后续要转成 dB通常是10*log10(energy / reference_power)参考量取决于你是做声压级还是振动加速度级。工程上我还会把每个频段的时域信号存到y_band_all{k}里方便后面看包络和回放噪声。3.3 滤波器组低频段的稳定性校验低于 50 Hz 的频带很容易掉进butter精度陷阱。一个频带的下限可能只有 44.7 Hz在Fs44100下归一化后只有 0.002极点在单位圆附近高度密集双精度浮点设计出来的系数会有较大误差。遇到这种情况先检查isstable(b,a)再检查滤波器频率响应是否真的落在目标频带上[H, w] freqz(b, a, 65536, Fs); [~, idx_max] max(abs(H)); fprintf(中心频率 %.2f Hz实际峰值 %.2f Hz\n, fc(k), w(idx_max));如果实际峰值和设计中心频率偏离超过带宽的 5%不要硬调阶数考虑把采样率先降下来或者改用designfilt(bandpassiir, FilterOrder, 4, HalfPowerFrequency1, f_low, HalfPowerFrequency2, f_high, SampleRate, Fs)让 MATLAB 自动优化极点布局。时域滤波器组这种方式适合频带数量不多、信号长度较大的场景逐频带滤波的循环开销在长出几十条频带时会明显上升。4. one_three_time_Feq.m 的频域路线FFT 后按掩码积分one_three_time_Feq.m和第二份代码的明显区别是不做滤波而是直接对整段信号做 FFT再用频带掩码把频率轴上的能量加起来。它的优点是计算量小只需要一次 FFT就可以算出全部中心频率对应的 1/3 倍频程频谱。缺点是时域信息被平均掉了无法观察某个频带能量的时间波动也不适合短时瞬态分析。4.1 功率谱密度与频带能量换算MATLAB 中直接算频带能量我常用的是“单边功率谱再累加”而不是直接取abs(fft(x))的平方。假设信号长度为N采样率Fs频率分辨率df Fs/N单边功率谱Pxx的计算方法是N length(x); X fft(x); Pxx abs(X(1:floor(N/2)1)).^2 / (Fs * N); Pxx(2:end-1) 2 * Pxx(2:end-1); % 单边补偿 f_axis (0:floor(N/2)) * Fs / N;代码里除以Fs*N是把离散 FFT 结果规整成物理意义上的功率谱密度单位是功率每赫兹。这里的2 * Pxx(2:end-1)是单边谱修正因为负频率的贡献被折叠到正频率来了直流分量和奈奎斯特分量不需要乘 2。如果你只关心相对能量不除以Fs*N也行但不同信号长度之间就不具备可比性。4.2 用掩码矩阵累加频带能量有了功率谱密度Pxx和掩码矩阵M频带能量可以一行矩阵运算取到fc oct3_fc(20, 20000); M oct3_mask(fc, f_axis); % 返回逻辑矩阵 E_band M * Pxx; % [numel(fc), 1] L_band 10 * log10(E_band / 1e-12); % 声压级示例参考 20 uPa矩阵乘法的含义是把一个频带内所有谱线对应的功率密度值做求和。注意这里没有乘df因为Pxx本身就是功率谱密度累加一条谱线时相当于在频率分辨率宽度内对密度做积分。如果不想引入二义性也可以直接用E_band M * (Pxx .* df)这取决于你把Pxx视为密度还是直接视为功率。我的习惯是前者写代码时在注释里固定单位方便后续换成 1/1 倍频程时只改掩码。4.3 泄漏、窗函数和频带边界谱线的处理频域路线最容易翻车的地方是频谱泄漏。矩形窗下正弦信号的能量会泄漏到旁瓣导致 1 kHz 的信号跑到 1250 Hz 频带里尤其在 FFT 点数没有整周期截断时。两个常用解决办法一是选择汉宁窗hann(N, periodic)并做幅值修正二是把信号分段做 Welch 平均这样得到的功率谱方差更小频带能量也更稳定。xw x(:) .* hann(N, periodic); X fft(xw); Pxx abs(X(1:floor(N/2)1)).^2 / sum(hann(N, periodic).^2);这里的归一化改成除以窗函数的平方和而不是Fs*N是因为加窗后信号总能量变少需要对窗本身的能量进行补偿。如果你拿pwelch直接算它内部会做类似的事情但返回的是密度形式通常要再乘以等效噪声带宽细节容易看漏。掩码边界上的谱线也应该处理。频率分辨率低于频带宽度时直接(f_axis f_low) (f_axis f_high)会漏掉边界附近半根谱线的能量。更稳的做法是先把频带上下限向外延伸半个频率分辨率再对被跨过的谱线按比例拆分。下面这段是我常用的比例拆分逻辑适合 5 年以上经验的人对照自查df Fs / N; band_low fc * 10^(-1/20); band_high fc * 10^(1/20); i1 floor(band_low / df); i2 ceil(band_high / df); frac_low 1 - (band_low/df - i1); frac_high band_high/df - i2 1; weight zeros(size(Pxx)); weight(i11) frac_low; weight(i21) frac_high; for k 1:numel(fc) E_band(k) sum(Pxx(i1(k)1:i2(k)1) .* weight(i1(k)1:i2(k)1)); end这里的i11和i21是为了应对 MATLAB 从 1 开始的下标。权重只有首尾谱线不是 1中间谱线直接按整根能量计入。低频段频带宽度窄比如中心频率 25 Hz 的带宽只有 5.6 Hz如果df大于 2 Hz边界拆分的影响会变得非常大这时不要省略这一步。频域法本质上把每个频带当作一个矩形滤波器窗函数的主瓣拖尾会在相邻频带之间造成串扰。行业里常说的“频带间最小衰减”指标用 FFT 掩码很难做到与时域 Butterworth 滤波器组同样的选择性。工程上我的经验是只用 FFT 掩码做普查和趋势观察出正式报告时至少用poctave交叉验证一次。5. A 计权修正、多通道批量处理与 1/3 倍频程可视化拿到各频带的线性声压级之后下一步通常不是直接画频谱而是叠加 A 计权修正。环境噪声标准里最终评价量基本是 A 计权等效声级 LAeq三分之一倍频程频谱最大的用处之一就是算 A 计权总声级和频带贡献度。把 A 计权做成一个函数和中心频率一一对应落在哪个频带就减哪个 dB。5.1 A 计权公式与查表实现MATLAB 里可以直接用weightingFilter对象但要做离线数组我更愿意用经典模拟公式function A a_weight(f) % f 为标量或向量单位 Hz返回 A 计权 dB 值 f abs(f); Ra 12194^2 * f.^4 ./ ... ((f.^2 20.6^2) .* sqrt((f.^2 107.7^2) .* (f.^2 737.9^2)) .* (f.^2 12194^2)); A 20*log10(Ra) 2.00; end2.00是对 1 kHz 参考点进行归一化后得到的常数直接带入时中心频率为 1000 Hz 的 A 计权值是 0.0016 dB约等于 0。实际使用中我会对每个中心频率提前算好权重表再在频带级L_A(k) L_linear(k) A(fc(k))即可。这个公式来自 IEC 61672适用于 20 Hz 到 20 kHz超出范围的值不要外推直接置 NaN 或剔出频带。5.2 多通道批量处理时的函数封装现场采集很多是 16 或 32 通道同步录制的。逐个跑脚本会非常慢我会先把第 4 章里的频域能量计算封装成单函数再把通道维度放到循环外层function bands process_multichannel(X_all, Fs, f_low, f_high) % X_all: N x Ch 的信号矩阵N 为采样点数Ch 为通道数 f_axis (0:floor(size(X_all,1)/2)) * Fs / size(X_all,1); fc oct3_fc(f_low, f_high); M oct3_mask(fc, f_axis); bands zeros(numel(fc), size(X_all, 2)); for ch 1:size(X_all, 2) N size(X_all, 1); w hann(N, periodic); Pxx abs(fft(X_all(:,ch) .* w)).^2; Pxx Pxx(1:floor(N/2)1) / sum(w.^2); Pxx(2:end-1) 2 * Pxx(2:end-1); bands(:, ch) M * Pxx; end end这样处理有一个好处掩码矩阵M在 32 个通道之间完全复用等于把最耗时的部分只算了两次。如果通道数再多还可以把M * Pxx改成Pxx. * M.利用 BLAS 加速不过 32 通道下差别不大。容易踩的坑是不同通道长度必须一致否则N变了f_axis也要重新生成不能对多通道信号用同一个df。5.3 用 imagesc 画随时间变化的 1/3 倍频程谱时间维度上的频谱图我喜欢用imagesc搭配对数频率轴。假设你已经按上面步骤算出了t_array时间轴和SPL(t, band)矩阵figure; imagesc(t_array, log10(fc), SPL_dB); set(gca, YDir, normal, YTick, log10(fc(1:5:end)), ... YTickLabel, fc(1:5:end)); ylabel(中心频率/Hz); xlabel(时间/s); colorbar;这里把fc取对数后作为图像纵轴坐标YTickLabel又显示原始频率值视觉上让 1/3 倍频程程带在对数频率上等距分布。常见的错误是直接用线性fc作为图像纵坐标视觉上高频带会很窄低频带又宽得没有参考意义。imagesc的坐标范围和矩阵行列要一致SPL_dB的行必须对应频率列对应时间否则图会横过来但坐标标签却不变。当频带数达到几十个时图例会非常挤。我一般只标注低频、中频和高频各取 1/3 点比如 20 Hz、63 Hz、125 Hz、1 kHz、4 kHz、12.5 kHz或者在纵轴上用MinorTick显示所有倍频程线。灰度图比彩色图在打印报告里更耐看彩色图则在屏幕上找异常频带更快。6. 用合成信号快速验证你的 1/3 倍频程实现找个已知信号验证手写代码是我每次改完参数后必做的一步。验证思路很简单生成两个幅度相同的正弦波一个落在 1 kHz 中心频带内一个落在 1.25 kHz 中心频带内对比两个 1/3 倍频程频带的能量差。理想情况下它们应该进入相邻的两个中心频带而且计算出的功率应该接近。Fs 48000; Dur 10; N Fs * Dur; t (0:N-1) / Fs; x 0.6 * sin(2*pi*1000*t) 0.6 * sin(2*pi*1259*t); bands process_multichannel(x, Fs, 20, 16000); % 频带能量矩阵 dB_band 10 * log10(bands / 1e-12); % 找出 1 kHz 和 1.25 kHz 对应的中心频率索引 [~, idx_1k] min(abs(fc - 1000)); [~, idx_125k] min(abs(fc - 1259)); fprintf(1kHz 频带 %.2f dB, 1.25kHz 频带 %.2f dB\n, ... dB_band(idx_1k), dB_band(idx_125k));正常结果应当是两个频带能量非常接近偏差小于 0.2 dB且比旁边频带至少高出 30 dB。如果偏差大先检查窗函数归一化是否正确如果旁边频带串扰大再看频率分辨率df是否太高把 10 秒信号改成 20 秒后重测。第二个验证建议是直接和 MATLAB 内置函数poctave对照。poctave(x, Fs, BandsPerOctave, 3)返回的是频带功率谱密度类型是PSD输出单位可能和你手写代码的单位差一个等效噪声带宽因此比较形状而不是绝对数值。下面这段可以同时画出两条曲线[P_builtin, fc_builtin] poctave(x, Fs, BandsPerOctave, 3); plot(log10(fc_builtin), 10*log10(P_builtin)); hold on; plot(log10(fc), dB_band); legend(poctave,手写掩码);poctave的谱线和手写结果的形状如果不重合优先检查中心频率序列是否一致以及掩码边界略小于理论值时低频带差值会变大。这里还有个容易被忽略的细节poctave默认返回的是带通信号功率谱密度需要乘以频带宽度再转声压级所以我上面把它和手写能量直接画在一起只是一种快速对齐测试正式报告里不要混用单位。验证全部通过之后再回去看one_three_time_Feq.m里的参数就能看出每一步到底做了什么事也能知道为什么不能随便删掉某一行。本文还有配套的精品资源点击获取