搬运自我(没人看)的博客:数字信号处理与音频——低通滤波器原理
因为要测试一下n久没用的博客,故公式可能会出现渲染错误,欢迎指正。
并且避免篇幅较长,我将原文拆开了几篇文章,如果要有完整的阅读体验,可以去博客看看
读本文可能需要一点信号与系统的基础和编程基础。
滤波在计算机中如何实现?
省流版:
double cutoff_Freq=432; //截止频率
double sampleRate=44100.0; //采样率
double a=cutoff_Freq/sampleRate; //滤波系数
for(int i=1;i<samples.length();i++){
samples[i]=a*samples[i]+(1.0-a)*samples[i-1];
}
公式为:
x [ i ] = a x [ i ] + ( 1 − a ) x [ i − 1 ] x[i]=ax[i]+(1-a)x[i-1] x[i]=ax[i]+(1−a)x[i−1]
这个公式怎么来的,a又是什么,下文会提到。
为什么不直接在频域做处理?
首先:傅里叶变换 的运算量极大,且涉及大量三角函数和浮点运算,经过ft(傅里叶变换)和ift(逆傅里叶变换)过程,耗时非常高;而且频谱是根据时间不断变化的,所以要用stft等进行分帧,再进行变换,而每次分帧加窗都会有误差产生,误差累积起来就会对音质有影响。
如何从时域去进行操作?
模拟电路?
我们从RC电路开始说起:
打开电路仿真网站: http://lushprojects.com/circuitjs/circuitjs.html 选择RC电路:
然后右键下方的波形图,选择properties(属性):
勾上show Spectrum,当然如果有其他想观测的信息也可以勾上;
现在可以对比一下,去掉电阻后的输出有什么不同(重点关注85.2Hz后的频谱,因为计算得出的截止频率为85.2Hz,这里请先关注现象,截止频率计算什么的将在后面提到)
这个电路就是RC滤波电路,电感以磁能存储电能,在电流变化时产生反向电压,阻止电流突变;电容存储电荷,当电压变化时电荷会流出或流入,减缓电压变化。这两种对交流电的阻碍都叫做阻抗。
电容的阻抗的公式为 Z C = 1 j w c Z_C = \frac{1}{jwc} ZC=jwc1
电压的阻抗的公式为 Z R = R Z_R = R ZR=R
根据KVL基尔霍夫电压定律: V o u t = V c = V i n − V R V_{out} = V_c = V_{in} - V_R Vout=Vc=Vin−VR
因为输入的是交流电,将其用微分方程表示 R C d v c ( t ) d t = v i n ( t ) − v c ( t ) RC \frac{dv_c(t)}{dt} = v_{in}(t) - v_c(t) RCdtdvc(t)=vin(t)−vc(t) (这里的公式大概能懂吧,有点像电路中电流的等式,可以结合高中学过的公式: C = Q U C = \frac{Q}{U} C=UQ , I = Q t = U R I = \frac{Q}{t}= \frac{U}{R} I=tQ=RU ,也可结合信号与系统——奥本海姆书P153来看)
令输入电压 V i n = e j w t V_{in} =e^{jwt} Vin=ejwt,输出电压一定为:$V_{out} =H(jw)e^{jwt} $。 其中 H ( j w ) H(jw) H(jw) 是 频率响应 ,即 H ( j w ) = v o u t v i n H(jw) = \frac{v_{out}}{v_{in}} H(jw)=vinvout
(为什么用到 e j w t e^{jwt} ejwt ?因为对 e j w t e^{jwt} ejwt 进行求导或积分,求导后的表达式仍有$ e^{jwt} $,利用复指数的性质能简化运算和方便研究系统特性)
得: R C d d t H ( j w ) e j w t = e j w t − H ( j w ) e j w t RC\frac{d}{dt}H(jw)e^{jwt} = e^{jwt} - H(jw)e^{jwt} RCdtdH(jw)ejwt=ejwt−H(jw)ejwt
可直接得出: H ( j w ) e j w t = 1 1 + R C j w e j w t H(jw)e^{jwt} = \frac{1}{1+RCjw}e^{jwt} H(jw)ejwt=1+RCjw1ejwt
其中上文提到的截止频率 w c = 1 R C w_c = \frac{1}{RC} wc=RC1
对频率响应取模: ∣ H ( j w ) ∣ = ∣ 1 1 + R C j w ∣ = 1 1 + ( w w c ) 2 |H(jw)| = |\frac{1}{1+RCjw}| = \frac{1}{ \sqrt{1 + {(\frac{w}{w_c})}^2} } ∣H(jw)∣=∣1+RCjw1∣=1+(wcw)21 ((实部平方+虚部平方)并开根号)
所以当信号的频率 w > w c w>w_c w>wc,则信号衰减;
当将$ H(jw) $推广到拉普拉斯域,即 s = σ + j w s=\sigma+jw s=σ+jw,得: H ( s ) = 1 1 + s w c H(s) = \frac{1}{1+ \frac{s}{w_c}} H(s)=1+wcs1
此时 H ( s ) H(s) H(s) 就叫做 传递函数 .
数字滤波
计算机处理的数据是离散的,音频也是离散化的。所以我们要利用z变换将传递函数离散化,利用一阶前向差分法:
s
=
z
−
1
T
s=\frac{z-1}{T}
s=Tz−1
差分法公式推导:一阶前向差分法公式推导
将上式代入到传递函数: H ( s ) = 1 1 + s w c H(s) = \frac{1}{1+ \frac{s}{w_c}} H(s)=1+wcs1,得
H ( z ) = Y [ z ] X [ z ] = 1 1 + z − 1 T w c H(z) = \frac{Y[z]}{X[z]} = \frac{1}{1+ \frac{z-1}{Tw_c}} H(z)=X[z]Y[z]=1+Twcz−11
为了照顾数学和我一样不太好的同学,我多写两步,将上式整理:
Y [ z ] + Y [ z ] z − Y [ z ] T w c = X [ z ] Y[z] +\frac{Y[z]z-Y[z]}{Tw_c} = X[z] Y[z]+TwcY[z]z−Y[z]=X[z]
根据上文提到过的 Z变换的时移性质,并从z域变换到时域:
y [ n ] + y [ n + 1 ] − y [ n ] T w c = x [ n ] y[n] +\frac{y[n+1]-y[n]}{Tw_c} = x[n] y[n]+Twcy[n+1]−y[n]=x[n]
令 T w c = a Tw_c = a Twc=a,此时a就是滤波系数,得
a y [ n ] + y [ n + 1 ] − y [ n ] = a x [ n ] ay[n] +y[n+1]-y[n] = ax[n] ay[n]+y[n+1]−y[n]=ax[n]
即
y [ n + 1 ] = a x [ n ] + ( 1 − a ) y [ n ] y[n+1] =ax[n] + (1-a)y[n] y[n+1]=ax[n]+(1−a)y[n]
公式基本推导完成。(和刚开始的滤波公式差不多对吧(不过说起来,刚开始提到的滤波公式应该采用的是后向差分进行离散化,读者可以自行证明一下,方法差不多的),另外实践中离散化一般不采用这种方法,还有一些方法如双线性变换法等)
下次将会用C++与Juce实现这个滤波器(官方示例好像也有,可以去翻翻Juce的示例)