如何计算滤波器对应特定dB电平的幅值用于FIR2滤波?
我使用fir2函数设计滤波器,再通过filter(b,1,audio)对WAV格式音频进行滤波,目标是让音频的各频谱频段匹配指定的目标dB值。目前手头有两组dB数据:
- 音频原始频谱的dB值
- 目标频谱的dB值
我尝试通过计算两组电平的差值并做逆对数转换来得到滤波器幅值,之前用的公式是:
20^(FdB SPL - OdB SPL/20)
FdB:目标电平,OdB:原始电平
比如针对125Hz频段,计算得到20^(38-35/20)=1.995,把这些数值作为fir2的幅值参数,使用50阶滤波器,但结果要么达不到预期,要么在不同音频文件中的精度差异很大。我用Artemis软件的"1/n Octave Spectrum (FFT) (4096; 50%)"功能计算所有dB值,怀疑是公式这类基础环节出了问题,求正确的实现方法。
原MATLAB代码
[audio,fs] = audioread(append(filename,".wav")) ; f = [0 0.0004535147392 0.000566893424 0.0007256235828 0.0009070294785 0.001133786848 0.001428571429 0.001814058957 0.002267573696 0.002857142857 0.003628117914 0.004535147392 0.00566893424 0.007256235828 0.009070294785 0.01133786848 0.01428571429 0.01814058957 0.02267573696 0.02857142857 0.03628117914 0.04535147392 0.0566893424 0.07256235828 0.09070294785 0.1133786848 0.1428571429 0.1814058957 0.2267573696 0.2857142857 0.3628117914 0.4535147392 0.566893424 0.7256235828 0.9070294785 1]; m = [1 0.88180561540975 0.881805363482077 0.835826382716575 0.776642454832629 0.759436720529334 0.653244752807778 0.745167995891424 0.850117073913788 0.98912260419711 0.960201464902893 0.923199850846219 1.13897797703032 0.874874253414851 1.13953780367958 1.28584396648708 0.67694058247949 0.785826825740988 1.00665129584544 1.38588097381939 1.59663824625406 1.46013949870038 1.55967029805846 1.57082839521638 1.96165693766572 1.74722265683908 1.77693364984782 1.29579584587398 1.58136046978984 1.69990075644448 1.42024073863697 1.91114522246185 1.55168355673292 1.00086374398705 1.0301145684381 1]; b = fir2(50, f, m6); A = 1; Y = filter (b, A, audio); amp = 0.06309573445 Y_new=Y*amp/rms(Y); newname = append("f_","_",filename,".wav"); audiowrite (newname, Y_new, fs,'BitsPerSample',32);
核心问题修正
1. 幅值转换公式错误
这是最关键的问题:你用了错误的逆对数转换公式,且运算优先级错误。
正确的dB转线性幅值公式是:
linear_gain = 10^((FdB - OdB)/20);
- 错误点1:应该用10的指数,不是20的指数——dB的定义是
20*log10(线性幅值),逆运算自然是10的幂。 - 错误点2:运算优先级错误——必须先计算
FdB - OdB的差值,再除以20,作为指数。你之前的公式是20^(FdB - (OdB/20)),完全不符合dB的转换逻辑。
2. 代码笔误
原代码中b = fir2(50, f, m6);的m6是未定义的变量,应该替换为你计算好的幅值数组m,否则代码直接报错。
其他优化建议
1. 频段匹配对齐
Artemis的1/n倍频程频谱是带通平均后的结果,而fir2是基于线性频率点的FIR滤波器设计。要确保f数组中的归一化频率点,和Artemis中倍频程的中心频率严格对应(归一化到fs/2),避免频率点错位导致滤波效果偏差。
2. 滤波器阶数调整
50阶FIR对于音频的宽频段(尤其是低频段)过渡带控制可能不足,建议用firpmord函数估算合适的滤波器阶数,保证每个倍频程频段内的幅值足够接近目标值:
% 估算阶数和参数,0.01是通带/阻带的容差 [N, Fo, Ao, W] = firpmord(f, m, [0.01 0.01]); b = firpm(N, Fo, Ao, W);
3. 信号归一化优化
原代码中Y*amp/rms(Y)的固定幅值归一化不合理,建议改为基于峰值的归一化,避免滤波后的信号出现削波:
max_val = max(abs(Y)); Y_new = Y / max_val * 0.9; % 乘以0.9是为了留余量,防止信号溢出
修正后的完整代码
[audio, fs] = audioread(append(filename, ".wav")); % 确保频率点与Artemis倍频程中心频率对应(归一化到fs/2) f = [0 0.0004535147392 0.000566893424 0.0007256235828 0.0009070294785 0.001133786848 0.001428571429 0.001814058957 0.002267573696 0.002857142857 0.003628117914 0.004535147392 0.00566893424 0.007256235828 0.009070294785 0.01133786848 0.01428571429 0.01814058957 0.02267573696 0.02857142857 0.03628117914 0.04535147392 0.0566893424 0.07256235828 0.09070294785 0.1133786848 0.1428571429 0.1814058957 0.2267573696 0.2857142857 0.3628117914 0.4535147392 0.566893424 0.7256235828 0.9070294785 1]; % 替换为用正确公式计算的幅值数组:10^((FdB - OdB)/20) m = [1 0.88180561540975 0.881805363482077 0.835826382716575 0.776642454832629 0.759436720529334 0.653244752807778 0.745167995891424 0.850117073913788 0.98912260419711 0.960201464902893 0.923199850846219 1.13897797703032 0.874874253414851 1.13953780367958 1.28584396648708 0.67694058247949 0.785826825740988 1.00665129584544 1.38588097381939 1.59663824625406 1.46013949870038 1.55967029805846 1.57082839521638 1.96165693766572 1.74722265683908 1.77693364984782 1.29579584587398 1.58136046978984 1.69990075644448 1.42024073863697 1.91114522246185 1.55168355673292 1.00086374398705 1.0301145684381 1]; % 可选:用firpmord估算合适阶数,替代固定50阶 % [N, Fo, Ao, W] = firpmord(f, m, [0.01 0.01]); % b = firpm(N, Fo, Ao, W); b = fir2(50, f, m); Y = filter(b, 1, audio); % 基于峰值的归一化,避免削波 max_val = max(abs(Y)); Y_new = Y / max_val * 0.9; newname = append("f_", filename, ".wav"); audiowrite(newname, Y_new, fs, 'BitsPerSample', 32);
内容的提问来源于stack exchange,提问作者Wuchta7

