https://www.dazhuanlan.com/2019/10/20/5dab47e6a4855/
窗函数法是一种设计FIR滤波器的方法,FIR全称为Finite impulse response,即有限冲激响应滤波器。 通带为(omega_cltpi)的理想低通滤波器的系统函数为 [ H_d(e^{jomega})=begin{cases} 1cdotp e^{jalphaomega} & |omega|leomega_e 0 & omega_eltpi end{cases} ] 反变换求它的冲激响应为 [ begin {aligned} h_d(n) & = mathcal{F}^{-1}bigl[H_d(e^{jomega})bigr] = frac{1}{2pi} int_{-pi}^pi H_d(e^{jomega})e^{jomega n}domega \ & = frac{1}{2pi} int_{-omega_e}^{omega_e} 1cdot e^{-jalphaomega} e^{jomega n}domega \ & = frac{sinleft[omega_e(n-alpha)right]}{pi(n-alpha)} end {aligned} ] 可以看出,这是一个sinc函数

可以看出,(h_d(n))的冲激响应在时域上是无限长的,不符合FIR滤波器的要求,所以需要对它截断,也就是加窗。 为了得到一个长度为M的因果的线性相位的FIR滤波器(h(n)),我们需要 [ h(n)=begin{cases} h_d(n) & 0le nle M-1 \ 0 & elsewhere end{cases} ~~ and ~~ alpha = frac{M-1}{2} ] 也就是让(h_d(n))的(0le nle M-1)之外的值全为0,这就是加窗。这个过程可以表示为 [ h(n) = h_d(n)w(n) ] (w(n))即为窗函数 [ w(n) = begin{cases} some symmetric function with respect to \ alpha over 0le nle M-1 \ 0, otherwise end{cases} ] 在这里我们将(w(n))定义为 [ w(n) = begin{cases} 1 & 0le nle M-1 \ 0 & otherwise end{cases} ] 再看看频域,使用(H(e^{jomega}))表示因果FIR滤波器的响应,(H_d(e^{jomega}))表示理想滤波器,(W(e^{jomega}))表示窗函数响应。 时域相乘,频域卷积,这里是圆周卷积。(圆周卷积以后也会介绍) [ begin {aligned} H(e^{jomega}) & = H_d(e^{jomega})bigodot W(e^{jomega}) \ & = frac{1}{2pi}int_{-pi}^pi W(e^{jomega})H(e^{j(omega - lambda)}) end {aligned} ] 我想在MATLAB中模拟这个加窗的过程,但是一直弄不出来,就把书上的图拿过来吧。

我不知道为什么我卷积出来的频域波形和窗函数的频域波形是一样的,很尴尬。 不过窗函数的原理基本还是理解了,直接生成书上卷积出来的频域波形不难,在matlab中直接生成一个FIR滤波器(相当于加了矩形窗),然后求它的频响。

代码如下: 1
2
3
4
| alpha = (M-1)/2;
n = [0:1  M-1)];
m = n - alpha + eps; % add smallest number to avoid divide by zero
hd = sin(wc*m) ./ (pi*m);
|
1
2
3
4
5
6
| wc = 0.2*pi;
M = 60;
hd = ideal_lp(wc,M);
w_rect = (rectwin(M))';
h = hd .* w_rect;
fvtool(h);
|
参考书籍: - _DigitalSignalProcessingUsingMatlabv4.0_JohnG.Proakis
- 《数字信号处理教程》程佩青
|