加窗FIR滤波器设计实验【内含实际使用filter的例子，讨论了filter、fftfilt、filtfilt的差别】_filtfilt和filter

转自：http://blog.sina.com.cn/s/blog_6e8d34350100ng3f.html

题一：

利用加窗傅里叶级数法，设计一个具有如下指标的线性相位FIR低通滤波器：通带截止频率在4rad/s处，阻带截止频率在6rad/s处，最大通带衰减为0.2dB，最小阻带衰减为42dB，抽样率为18rad/s。利用下面的各个窗函数进行设计：海明窗、汉宁窗和布莱克曼窗。对于每种情况，给出冲击响应的系数并画出设计的滤波器的增益响应。分析你的结果。不要使用M文件fir1。

分析：由于题中给定的频率为实际频率，而数字滤波器一般使用归一化频率，故应该先对给定频率信息进行归一化。

%程序设计：

close all

Omegap=4;

Omegas=6;

Omegat=18;

wp=2*pi*Omegap/Omegat;%0.44pi

ws=2*pi*Omegas/Omegat;%0.667pi %归一化角频率

alphap=0.2;

alphas=42;

 

NUM_Hamming=1;

NUM_Hann=2;

NUM_Blackman=3;

c=[3.32,3.11,5.56]*pi;

wc=(ws+wp)/2;

delt_w=ws-wp;

freq_labels={'用海明窗设计的FIR频率响应','用汉宁窗设计的FIR频率响应',...

    '用布莱克曼窗设计的FIR频率响应'};

ht_labels={'海明窗FIR冲激响应','汉宁窗FIR冲激响应','布莱克曼窗FIR响应'};

for filter_kind=NUM_Hamming:NUM_Blackman,

    M=ceil(c(filter_kind)/delt_w);%向上取整

    N=2*M+1;

    switch filter_kind,

        case NUM_Hamming,

            win=hamming(N);

            display(['海明窗生成的冲激响应系数：','  (阶数N=',num2str(N),')']);

            figure

        case NUM_Hann,

            win=hann(N);          

            display(['汉宁窗生成的冲激响应系数：','  (阶数N=',num2str(N),')']);

            figure

        case NUM_Blackman,

            win=blackman(N);           

            display(['布莱克曼窗生成的冲激响应系数：','(阶数N=',num2str(N),')']);

            figure

        otherwise disp('error');

    end

    n=-M:M;

    hd=sin(wc*n)./(pi*n);

    hd(find(n==0))=wc*cos(wc*0)/pi;

        %因为n=0时公式“sin(wc*n)./(pi*n)”中分母为0，计算结果为NaN,所以应该另用

        %一种方法单独计算n=0时的hd值。这里采用洛必达法则对原hd公式进行分子分母微分

        %得hd=wc*cos(w%c*n)./pi，n->0;这样可以求得当n=0时原函数值。(也可以将n加上eps

        %—-系统精度，使n最大程度接近0而不等0）

    ht=hd.*win';

    display(['',num2str(ht)]); %打印显示冲激响应系数

    %wvtool(ht);   

    subplot(1,2,1);

    plot(n,ht,'.-') %画冲激响应系数图

    title(ht_labels(filter_kind));

    xlabel('n','FontSize',12);

    ylabel('ht','fontsize',12);

    grid on

 

    [h,w]=freqz(ht,1,512);

    W=w/pi;

    H=20*log10(abs(h));

    subplot(1,2,2); %画频响图

    hold on

    title(freq_labels(filter_kind));

    plot(W,H);

    xlabel(['\pi (\omega/\omega','s)']);

    ylabel('增益(dB)');

    grid on  

    %为了更直观的看到所设计FIR是否达到指标，可分别在图中标出相关点

    %从而可以对不同窗函数进行对比。标示过程如下：

    dotp=round(mean(find(w>wp-0.1&w<wp+0.1)));%找到点wp的坐标

    dots=round(mean(find(w>ws-0.1&w<ws+0.1)));%找到点ws的坐标

    plot(W(dotp),H(dotp),'.r','MarkerSize',25);%对应处打上显目的红点

    plot(W(dots),H(dots),'.r','MarkerSize',25);

    text(W(dotp),H(dotp),['[','',num2str(round(W(dotp)*1000)/1000),...

        ',',num2str(round(H(dotp)*1000)/1000),'dB]'],'FontSize',15);%标示wp在频响图

%中的位置，从中可以%看来相应增益

    text(W(dots),H(dots),['[','',num2str(round(W(dots)*100)/100),...

    ',',num2str(round(H(dots)*100)/100),'dB]'],'fontsize',15);  %标示ws在频响图

%的位置

    hold off

end

 

 

%产生的各种窗函数的系数：

海明窗生成的冲激响应系数：  (阶数N=31)

0.0014702  -0.0013161   -0.001885   0.0038559   0.0022981  -0.0097177 3.3611e-017   0.019275  -0.0091462   -0.031341    0.031509    0.043366   -0.083816   -0.052269     0.31032     0.55556     0.31032   -0.052269   -0.083816    0.043366   0.031509   -0.031341  -0.0091462    0.019275 3.3611e-017  -0.0097177   0.0022981  0.0038559   -0.001885  -0.0013161   0.0014702

汉宁窗生成的冲激响应系数：  (阶数N=29)

-0  -0.0001973   0.0011375   0.0010796  -0.0059013 2.3913e-017    0.015232  -0.0077763   -0.028084    0.029338    0.041522   -0.081865   -0.051739     0.30954    0.55556     0.30954   -0.051739   -0.081865    0.041522    0.029338   -0.028084 -0.0077763    0.015232 2.3913e-017  -0.0059013   0.0010796   0.0011375  -0.0001973          -0

布莱克曼窗生成的冲激响应系数：(阶数N=53)

-1.6732e-019 -5.7522e-006 -6.1629e-005  0.00011007  0.00021128 -0.00048407  -0.0003015      0.0013 -9.0902e-018   -0.002622   0.0012439   0.0042178  -0.0041222  -0.0053512  0.0092484     0.00464   -0.016847 5.1274e-017    0.026475   -0.011539   -0.036934    0.035186    0.046454   -0.087054   -0.053145      0.3116     0.55556     0.3116   -0.053145   -0.087054    0.046454    0.035186   -0.036934   -0.011539   0.026475 5.1274e-017   -0.016847     0.00464   0.0092484  -0.0053512  -0.0041222  0.0042178   0.0012439   -0.002622 -9.0902e-018      0.0013  -0.0003015 -0.00048407 0.00021128  0.00011007 -6.1629e-005 -5.7522e-006 -1.6732e-019

 

 

 

海明窗所得FIR的冲激响应图及对应频率响应图：

其中，红点标示处为wp：[0.445,-0.027dB] ;ws :[0.67,-50.79dB] 

 

 
 

汉宁窗所得FIR的冲激响应及相应频率响应：

Wp：[0.445,-0.067dB] ;ws :[0.64,-45.07dB] 

 

布莱克曼窗所得FIR的冲激响应及相应频率响应：

Wp :[0.445,-0.001dB] ;ws :[0.67,-78.18dB]  

从以上三图中可以看出，对于给定的指标[-0.2dB,42dB]，各窗函数所得FIR在[wp,ws]处的增益分别为：海明窗[-0.027dB,-50.79dB] 

                       汉宁窗[-0.067dB,-45.07dB]

                       布莱克曼窗[-0.001dB,-78.18dB]

对比原指标，可以看出，汉宁窗比较接近指标，同时其阶数（N=29）也是最少的。但就衰减性能而言，显然海明窗的-50.79dB及布莱克曼窗的-78.18dB比汉宁窗的-45.07dB要好。从以上三图还可以看出，海明窗阻带波纹较大，而布莱克曼窗阻带衰减较好但波纹较密，汉宁窗处于海明窗与布莱克曼窗的中间。

  综合上述特点，可以得出如下结论：对于加窗法FIR滤波器的设计，如果是要设计要求不高的滤波器，可以采用海明窗，因为海明窗阻带衰减较好但阻带波纹相对较大，而且滤波器阶数相对不大。如果是要设计要求较高同时又不计成本及复杂性的FIR滤波器可以采用布莱克曼加窗，因为布莱克曼加窗法可以产生较好的阻带衰减性能。虽然布莱克曼阻带波纹较密，但随频率的增加波纹衰减越厉害，这对阻带衰减有一定好处。不过，在阶数方面，布莱克曼窗的阶数又是最大的，这就会造成系统复杂性及成本的增加。综合考虑阻带效果及成本、复杂性等因素，加窗FIR设计最好采用汉宁窗，这样可以在阶数及阻带衰减上有一个较好的取舍。

 

 

 

 

 

 

 

 

题二：

使用M文件fir1，设计一个48阶FIR带通滤波器，通带边界的归一化频率为0.35和0.65。假设一个信号，其中含f1=1Hz，f2=10Hz，f3=20Hz，三种频率成分信号的采样频率为50Hz。使用fftfilt函数滤波，将原信号与通过滤波器的信号进行比较。

 

%程序设计：

close all

w1=0.35;

w2=0.65;

wn=[w1,w2];

N=48;

b=fir1(N,wn,'bandpass');%得FIR系数

fvtool(b,1); %可视化数字滤波器工具，可以同时得到频响、相位、群延时、相位函数等图形

f1=1;

f2=10;

f3=20;

fs=50;

% Omega1=2*pi*f1/fs;

% Omega2=2*pi*f2/fs;

% Omega3=2*pi*f3/fs;%归一化角频率

% t=-150:150;%频率已归一化，故间隔应为1.

% y1=cos(Omega1*t)+cos(Omega2*t)+cos(Omega3*t);

 

t=-1:1/fs:1;

y1=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t); %这里生成y1的方法也可以用上面注释掉的部分的方法，实践证明这两种方法所得波形是一样的（但前者时间轴必须进行一定的调整）

 

y2=fftfilt(b,y1); %滤波

 

subplot(2,1,1);

plot(t,y1,'.-');

title('原信号');

xlabel('t');

ylabel('y1');

grid on

subplot(2,1,2);

plot(t,y2,'.-');

title('滤波后');

xlabel('t');

ylabel('y2');

grid on

 

 

%为更直观地对滤波前后信号进行分析，下面进行频谱分析

figure

Len=length(y1);

NFFT=2^(nextpow2(Len));

Y1=fft(y1,NFFT); %快速傅立叶变换

Y2=fft(y2,NFFT);

f=((-NFFT/2:NFFT/2-1)/NFFT)*fs;%频率轴调整使之与实际频率对应

subplot(2,1,1)

semilogy(f,fftshift(abs(Y1))/NFFT,'.-r');

title('原信号y1的频谱','color','red','fontsize',15);

ylabel('频谱幅度（log）');

xlabel('频率（Hz）','color','red','fontsize',15);

set(gca,'xlim',[0,f(end)],'ylim',[0,max(abs(Y1)/NFFT)]);

grid on

subplot(2,1,2);

semilogy(f,fftshift(abs(Y2))/NFFT,'.-r');

title('滤波后的y2的频谱','color','red','fontsize',15);

ylabel('频谱幅度（log）');

xlabel('频率（Hz）','color','red','fontsize',15);

set(gca,'xlim',[0,f(end)],'ylim',[0,max(abs(Y1)/NFFT)]);

grid on

 

 

所设计的FIR滤波器的频率响应：

 

 

 

 

相位响应：（从下图可以看出FIR滤波器优越的线性相位性能）

 

 

 

 

原信号及滤波后的波形示于如下：

 

从上图可以得到如下引申及结论：

（1）、按照题目要求，Wn分别为归一化的：0.35(pi rad/s)、0.65(pi rad/s)，按照归一化公式： 

可知， 、 

ω1=0.35pi，ω2=0.65pi，代入得：F1=8.75Hz、F2=16.25Hz。综上，对于通带为[ω1, ω2]的带通滤波器，题中原信号经过滤波后理论上被保留的应为10Hz的信号。即上图中“滤波后”理论上应该为10Hz。从图中可以看出，事实也是如此的。

（2）、从上图可以看出：“滤波后”的[-1,0.6]区间，波形出现了严重失真。原本应该为近似三角波的波形却变为平坦地近乎直线。对于此现象刚开始很费解，故我们对已经接触过的三个滤波函数（filter、filtfilt、fftfilt）进行了探索，以试图找出三个滤波函数的本质区别及产生图中失真的原因。首先，我们将本次实验程序中的fftfilt分别用其他两个函数代替，发现：对于filter函数，其结果与fftfilt的结果（无论是时域还是频域上）是一样的，而当使用filtfilt函数代替时，出现如下错误信息：“Data must have length more than 3 times filter order.”意为使用filtfilt时，数据长度必须大于滤波器阶数的3倍，改变长度又可以运行了.这说明程序中数据长度是有条件限制的.为更深层地找出它们的区别，我们通过找help文档及网上搜集的资料得到如下结论：

   这三个函数的不同之处基本在于其实现算法的差异。

         (a)filter，名曰：“Generic filtering function”又称一维数字滤波器（相对于filter2而言）、直接II型滤波器，也即为通用的数字滤波器。其算法是通过如下公式“a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb) - a(2)*y(n-1) - ... - a(na+1)*y(n-na)”进行滤波处理的。使用的是递推算法，引用网上的说法如下：

“ 用filter([1,1],[1,-0.5,0.06],x);来计算系统响应时，是假设系统初始状态为全0

首先计算x(0)输入系统的响应：y(0)=x(0)
然后 y(1)=0.5y(0)+x(1)+x(0)
接着 y(2)=0.5y(1)-0.6y(0)+x(2)+x(1)
以此类推，递推的算出所有x（n）对应的系统输出y(n)
这样算出的y和x是等长度的”

这也就不难理解图中的滤波波形了。本来对加窗FIR滤波器的系数ht就近似于sinc函数，其前面几个系数（如海明窗生成的冲激响应系数：  (阶数N=31)

0.0014702  -0.0013161   -0.001885   0.0038559   0.0022981  -0.0097177 3.3611e-017   0.019275  -0.0091462   -0.031341…… ）都很小，甚至趋近于0，这样在进行滤波时，由于x(n)要与系数相乘，自然所得结果刚开始会很小，由于使用递推算法，当滤波器处理的数据增加时，相加和的波动就越来越大，这才会渐渐出现正常的滤波波形。

    (b)filtfilt，名曰：零相位数字滤波。其滤波算法是基于filter而来的。只是filtfilt实现了零相位。其基本实现过程为先让信号用filter滤波，再将信号时域反转再次通过filter滤波，由《数字信号处理》的知识可知，这样两次滤波后幅度响应为单个filter幅度的平方而相位实现了零相位。因为同样使用filter的算法，故在使用filtfilt滤波后，波形同样会出现filter滤波的失真现象，即刚开始时输出信号近似为直线。

     (c)fftfilt,名曰：基于FFT和重叠相加法的FIR滤波。英文描述为：“Filters with an FIR filter using the FFT” matlab描述为：“Overlap-add method of FIR filtering using FFT.”从网上资料来看，fftfilt似乎为专门的FIR滤波函数并一般应用于音频信号处理中。按照网上的说法，fftfilt只对FIR滤波器有效（从其输入系数也可以看出这一点）。但对于具体的fftfilt实现公式无论是matlab自带的help还是网上的介绍资料都比较少。故在此未能对fftfilt函数做更多讨论。

     综合上述资料，对三个函数的总体感觉是：filter既能进行IIR滤波又能进行FIR滤波，而fftfilt虽然名曰“基于FFT和重叠相加法的FIR滤波”但总体感觉与filter相差无几，只不过fftfilt只能对FIR滤波器有效。而对于filtfilt来说，它是一个使用了filter、比filter略高一个层次的滤波器，它可以实现零相位，这是它的优点，但对于普通滤波器来说，我们要知道的正是滤波器的幅度和相位特点，实现零相位会掩盖滤波器的一些原有内涵，故filtfilt没有filter更广的适用性，它应该说是为特意达到零相位滤波而设计的。通过以上的学习，我们知道，已经接触到的这三个滤波函数在功能各有千秋，虽没有完全找到相关资料，但基本上都从原理上认识了三个滤波函数，这为以后的使用提供了很重要的参考。

 

 

做完了相关滤波函数的探讨这个插曲后，现在让我们把内容继续回归到本次实验的后续分析部分上来：

下图是为了更准确地分析滤波前后信号特点而专门增加的频谱分析部分。图中示出了滤波前后的频谱分布，两个图均以实际频率为横坐标。

从第一个图可以看出原信号由于取样频率的影响存在许多噪声纹波，但由于抽样率仍符合奈奎斯特准则，故原信号中的几个频率不会受太多影响，也就是说，虽然图中有这么多的噪声纹波但信号仍是能有效表达原信号的。

在此我们重点分析第二个图形（滤波后频谱）。从第二个图可以看出在0到10Hz，频谱比较平缓而在10Hz到20Hz却出现大量波纹。这点是否符合理论呢？首先，让我们从频响图（本题第一个图）上分析。从频响图可以看出，该带滤波器大致以0.5pi rad/s（归一化）为中心，左右成对称特点。而下面第一个图（原信号）也大致以10Hz呈一定的对称性，而10Hz恰好位于带通滤波器的通带内。由公式：Y(jw)=X(jw)H(jw)可知，滤波后频谱实际上是频率响应函数与原信号频谱函数的乘积，那么这样一来，从理论上来说，本题中信号滤波后频谱图也应该有一定的对称性。也就是说，在滤波后频谱图中，0到10Hz范围内也应该像10到20Hz波形相似有大量的波纹，而实际却是平缓的，这是为什么呢？

这里首先引用上面刚才对filter、fftfilt等三个函数进行分析的相关内容。从上面的分析可知使用这些函数进行滤波时，由于算法上的原因，滤波后信号在最前面的部分信号是平坦的、趋于直线的失真信号，那么当对该段信号进行频谱分析时，自然得到的频率就是比较稳定的“低频”成份了，频谱图中的低频部分就主要是这些因素造成的。

另外这里还有一种理解：由于题中滤波设计的阻带位于[0.35,0.65]之间，转换为以50Hz为取样率的实际频率就是：[8.75,16.25]Hz之间，可见，对于原信号中包含的10Hz频率，其分布靠近8.75Hz截止频率，也就是说，10Hz以下的频率衰减很大，而10到16.25Hz的噪声基本没被衰减，这也就造成滤波后10Hz以下频谱平滑而10Hz以上波纹较多这一“不对称”现象。个人觉得这种理解更有道理些。

当然，从本图也可以更清楚地看到原信号1Hz与20Hz频率成份基本被滤除而10Hz频率成分被保留的这一结果。

 

 

 

 作者：华南理工大学08电信微电2班---飞鸟