有限长单位脉冲响应滤波器设计

1、实验五 有限长单位脉冲响应滤波器设计1、 实验目的1、 掌握用窗函数法、频率采样法以及优化设计法设计FIR滤波器的原理及方法，熟悉相应的MATLAB编程。2、 熟悉线性相位FIR滤波器的幅频特性和相频特性。3、 了解各种不同窗函数对滤波器性能的影响。 2、 实验原理window=ones(1, N): 产生N点矩形窗，行向量。window=hann(N): 产生N点汉宁窗，列向量。window=hanning(N): 产生N点非零汉宁窗，列向量。等价于去除hann(N+2)的第一个零元素和最后一个零元素，得到的N点非零窗函数。window=hamming(N): 产生N点海明窗，列向量。win

2、dow=blackman(N): 产生N点布莱克曼窗，列向量。window=kaiser(N, beta): 产生参数为beta的N点凯塞窗，列向量。M, Wd, beta, ftype=kaiserord(f, a, dev, fs): 凯塞窗参数估计。f为一组边界频率，最高频率为fs/2。a为f中各个频带的幅度值，通带取1，阻带取0。如果f中有2个元素，则形成3个频带，其中第1个和第3个是通带或阻带，第2个是过渡带，a中也有2个元素，指明第1个和第3个频带是通带还是阻带；如果f中有4个元素，则形成5个频带，其中1，3和5是通带或阻带，2和4是过渡带，a中有3个元素，指明1，3和5是通带还是

3、阻带。dev的维数与a相同，指明每个频带上的波动值。fs为采样频率。M为FIR滤波器的阶数，M=N-1。Wd为归一化边界频率，等于数字边界角频率除以，或者边界频率除以fs/2。beta就是凯塞窗的参数。ftype为滤波器的类型。b = fir1(M, Wd, 'ftype', window): 用窗函数法求FIR滤波器的系数b(单位脉冲响应)。M为滤波器的阶数，M=N-1。Wd为一组归一化边界频率，通带和阻带间隔分布，无过渡带；只有一个元素，表示低通或高通滤波器；有两个元素表示带通和带阻滤波器；有三个及以上元素，表示多带滤波器。'ftype'表示滤波器类型，&#

4、39;high'表示高通滤波器，'stop'表示带阻滤波器，'DC-0'表示多带滤波器的第一个频带为阻带，'DC-1'表示多带滤波器的第一个频带为通带。window为窗口类型，缺省为海明窗。b = fir2(M, f, m, window): 用频率采样法求FIR滤波器的系数b。M为滤波器的阶数，M=N-1。f为一组归一化频率，第一个元素必须为0，最后一个元素必须为1(对应奈奎斯特频率，即采样频率的一半)，中间的元素按升序排列。m的维数与f相同，指明f中每个频率上的理想幅度。window为窗口类型，缺省为海明窗。Fir2可以实现任意幅度特

5、性的滤波器。三、实验内容1、用窗函数法设计一个线性相位FIR高通滤波器，通带边界频率为0.7，阻带边界频率为0.5，要求阻带衰减不小于50dB，窗函数从矩形窗、汉宁窗、海明窗和布莱克曼窗中选取，且要求滤波器的阶数最小。给出该滤波器的单位脉冲响应h(n)的解析式，并用MATLAB绘出时域波形和幅频特性。2、用窗函数法设计一个线性相位FIR带通滤波器，通带边界频率为3kHz和7kHz，阻带边界频率为2kHz和9kHz，采样频率为20kHz，要求阻带衰减不小于40dB，窗函数从矩形窗、汉宁窗、海明窗和布莱克曼窗中选取，且要求滤波器的阶数最小。给出该滤波器的单位脉冲响应h(n)的解析式，并用MATLA

6、B绘出时域波形和幅频特性。3、用窗函数法设计一线性相位FIR带阻滤波器，通带边界频率为600Hz和1400Hz，阻带边界频率为700Hz和1200Hz，采样频率为4000Hz，要求阻带衰减不小于50dB，窗函数从矩形窗、汉宁窗、海明窗和布莱克曼窗中选取，且要求滤波器的阶数最小。给出该滤波器的单位脉冲响应h(n)的解析式，并用MATLAB绘出时域波形和幅频特性。4、分别用矩形窗、汉宁窗、海明窗和布莱克曼窗设计FIR低通滤波器，理想滤波器的边界频率为0.8，N=51，用Fir1函数设计。根据所设计滤波器的阻带最小衰减值确定实际阻带边界频率。5、用频率采样法设计一个线性相位低通滤波器，N=40，通带

7、边界频率为/4，过渡带设置一个采样点|H(k)|=0.5，给出单位脉冲响应h(n)的序列值，并作出幅频响应。6、用频率采样法设计一个线性相位高通滤波器，通带边界频率为2/3，过渡带设置一个采样点|H(k)|=0.5，分别求N=31和N=32时的单位脉冲响应h(n)的序列值，并分别作出幅频响应。7、用频率采样法设计一个多带线性相位滤波器，理想幅频特性如下图所示，N=50，过渡带设置一个采样点，过渡点的值设为0.5，窗函数选矩形窗。给出该滤波器的单位脉冲响应h(n)的解析式，并用MATLAB绘出时域波形和幅频特性。4、 实验结果与分析 1. 解：求通带和阻带的数字边界频率:求理想高通滤波器的边界频

8、率。理想高通滤波器的边界频率是实际滤波器幅度函数正负肩峰频率的中心，但这里不知道正负肩峰的频率，因此用通带和阻带边界频率的中心近似:选择窗函数和窗口长度。阻带衰减不小于50dB，因此选择海明窗(其阻带衰减为53dB，满足衰减要求)。滤波器的过渡带宽为，因此窗口长度为:线性相位延迟常数为:根据理想边界频率和线性相位延迟常数，求理想单位脉冲响应:窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的单位脉冲响应:时域波形：幅频特性：程序：N=33;n=0:1:N-1;wn=0.6*pi;a=(N-1)/2;hd=(sin(pi*(n-a)-sin(n-a)*wn)./(pi*(n-a);hd

9、(a+1)=0.4;win=hamming(N);h=win'.*hd;figure;stem(n,h);xlabel('n');ylabel('h(n)');grid;title('线性相位高通滤波器FIR的单位脉冲响应h(n)');H,w=freqz(h,1);H=20*log10(abs(H);figure;plot(w/pi,H);axis(0 1 -100 0);xlabel('频率/Hz');ylabel('幅度/dB');title('线性相位高通滤波器FIR，海明窗，N=33'

10、;);2.解：求通带和阻带的数字边界频率:求理想带滤波器的边界频率。理想带通滤波器的边界频率是实际滤波器幅度函数正负肩峰频率的中心，但这里不知道正负肩峰的频率，因此用通带和阻带边界频率的中心近似:选择窗函数和窗口长度。阻带衰减不小于40dB，且阶数最小，因此选择汉宁窗(其阻带衰减为44dB，满足衰减要求)。滤波器的过渡带宽为,选择较小过渡带宽，因此窗口长度为:线性相位延迟常数为:根据理想边界频率和线性相位延迟常数，求理想单位脉冲响应:窗函数与理想单位脉冲响应相乘，得到线性相位FIR带通滤波器的单位脉冲响应:时域波形：幅频特性：程序：N=31;n=0:1:N-1;w1=0.25*pi;w2=0.

11、8*pi;a=(N-1)/2;hd=(0.5-0.5.*cos(pi/15*n).*(sin(w2*(n-a)-sin(w1*(n-a)./(pi*(n-a);hd(a+1)=0.55;win=hanning(N);h=win'.*hd;figure;stem(n,h);xlabel('n');ylabel('h(n)');grid;title('线性相位带通滤波器FIR的单位脉冲响应h(n)');H,w=freqz(h,1);H=20*log10(abs(H);figure;plot(w/pi,H);axis(0 1 -100 0);x

12、label('频率/Hz');ylabel('幅度/dB');title('线性相位带通滤波器FIR，汉宁窗，N=31');3.解：求通带和阻带的数字边界频率:求理想带滤波器的边界频率。理想带阻滤波器的边界频率是实际滤波器幅度函数正负肩峰频率的中心，但这里不知道正负肩峰的频率，因此用通带和阻带边界频率的中心近似:选择窗函数和窗口长度。阻带衰减不小于50dB，且阶数最小，因此选择海明窗(其阻带衰减为53dB，满足衰减要求)。滤波器的过渡带宽为,选择较大过渡带宽，由于选择的是带阻滤波器，N必须为奇数，因此窗口长度为:线性相位延迟常数为:根据理想边界频

13、率和线性相位延迟常数，求理想单位脉冲响应:窗函数与理想单位脉冲响应相乘，得到线性相位FIR带通滤波器的单位脉冲响应:时域波形：频域特性：程序：N=67;n=0:1:N-1;w1=0.35*pi;w2=0.65*pi;a=(N-1)/2;hd=(0.54-0.46.*cos(pi/33*n).*(sin(pi*(n-a)+sin(w1*(n-a)-sin(w2*(n-a)./(pi*(n-a);hd(a+1)=0.7;win=hanning(N);h=win'.*hd;figure;stem(n,h);xlabel('n');ylabel('h(n)');

14、grid;title('线性相位带阻滤波器FIR的单位脉冲响应h(n)');H,w=freqz(h,1);H=20*log10(abs(H);figure;plot(w/pi,H);axis(0 1 -90 0);xlabel('频率/Hz');ylabel('幅度/dB');title('线性相位带阻滤波器FIR，海明窗，N=33');4.实验结果：程序：N=51;M=N-1;n=0:1:M;Wd=0.8;w1=ones(N,1);w2=hanning(N);w3=hamming(N);w4=blackman(N);b1=fir

15、1(M,Wd,'low',w1);b2=fir1(M,Wd,'low',w2);b3=fir1(M,Wd,'low',w3);b4=fir1(M,Wd,'low',w4);figure;stem(n,b1);xlabel('n');ylabel('h(n)');grid;title('FIR低通滤波器的单位脉冲响应h(n)，矩形窗');H,w=freqz(b1,1);H=20*log10(abs(H);figure;plot(w/pi,H);axis(0 1 -100 0);xlab

16、el('频率/Hz');ylabel('幅度/dB');title('FIR低通滤波器，矩形窗');figure;stem(n,b2);xlabel('n');ylabel('h(n)');grid;title('FIR低通滤波器的单位脉冲响应h(n)，汉宁窗');H,w=freqz(b2,1);H=20*log10(abs(H);figure;plot(w/pi,H);axis(0 1 -100 0);xlabel('频率/Hz');ylabel('幅度/dB');

17、title('FIR低通滤波器，汉宁窗');figure;stem(n,b3);xlabel('n');ylabel('h(n)');grid;title('FIR低通滤波器的单位脉冲响应h(n)，海明窗');H,w=freqz(b3,1);H=20*log10(abs(H);figure;plot(w/pi,H);axis(0 1 -100 0);xlabel('频率/Hz');ylabel('幅度/dB');title('FIR低通滤波器，海明窗');figure;stem(n,

18、b4);xlabel('n');ylabel('h(n)');grid;title('FIR低通滤波器的单位脉冲响应h(n)，布莱克曼窗');H,w=freqz(b4,1);H=20*log10(abs(H);figure;plot(w/pi,H);axis(0 1 -100 0);xlabel('频率/Hz');ylabel('幅度/dB');title('FIR低通滤波器，布莱克曼窗');5.解：由于(2/40)×5<0.25，因此下边界在4和5之间，上边界在6和35之间。又因为

19、N为偶数，属于第二类线性相位滤波器，幅度函数在02上以为中心呈奇对称分布。因此频率采样值的幅度为:频响采样的相位为:频率采样值H(k)为：单位脉冲响应h(n)为：实验结果：程序：clear;N=40;k=0:N-1;Hk=ones(1,5) 0.5 zeros(1,29) -0.5 -ones(1,4);h1=real(ifft(Hk.*exp(-j*pi*(N-1)*k/N);h2=hanning(N)'.*h1;H1,w=freqz(h1,1);H2,w=freqz(h2,1);figure;plot(w/pi,20*log10(abs(H1),':r',w/pi,

20、20*log10(abs(H2),'-b');axis(0 1 -70 10); grid; xlabel('omega/pi'); ylabel('幅度/dB');legend('矩形窗','汉宁窗'); title('频率采样法，N=40，过渡带设一个采样点');6.解：（1）由于(2/31)×10<2/3，因此下边界在9和10之间，上边界在11和20之间。又因为N为奇数，高通，属于第一类线性相位滤波器，幅度函数在02上以为中心呈偶对称分布。因此频率采样值的幅度为:频响采样的相位

21、为:频率采样值H(k)为：单位脉冲响应h(n)为：（2）由于(2/32)×10<2/3，因此下边界在9和10之间，上边界在11和20之间。又因为N为偶数，高通，属于第四类线性相位滤波器，幅度函数在02上以为中心呈偶对称分布。因此频率采样值的幅度为:频响采样的相位为:频率采样值H(k)为：单位脉冲响应h(n)为：实验结果：程序：N=31;n=0:N-1;k=0:N-1;Hk=zeros(1,10) 0.5 ones(1,10) 0.5 zeros(1,9);h1=real(ifft(Hk.*exp(-j*pi*(N-1)*k/N);h2=hamming(N)'.*h1;f

22、igure;stem(n,h1);xlabel('n');ylabel('h(n)');grid;title('线性相位高通滤波器的单位脉冲响应h(n)，N=31');H1,w=freqz(h1,1);H2,w=freqz(h2,1);figure;plot(w/pi,20*log10(abs(H1),':r',w/pi,20*log10(abs(H2),'-b');axis(0 1 -70 10); grid; xlabel('omega/pi'); ylabel('幅度/dB')

23、;legend('矩形窗','海明窗'); title('频率采样法，N=31，过渡带设一个采样点');clear;N=32;n=0:N-1;k=0:N-1;Hk=zeros(1,10) 0.5 ones(1,11) 0.5 zeros(1,9);h1=ifft(Hk.*exp(-j*pi*(N-1)*k/N);h2=hamming(N)'.*h1;figure;stem(n,h1);xlabel('n');ylabel('h(n)');grid;title('线性相位高通滤波器的单位脉冲响应h(n

24、)，N=32');H1,w=freqz(h1,1);H2,w=freqz(h2,1);figure;plot(w/pi,20*log10(abs(H1),':r',w/pi,20*log10(abs(H2),'-b');axis(0 1 -120 10); grid; xlabel('omega/pi'); ylabel('幅度/dB');legend('矩形窗','海明窗'); title('频率采样法，N=32，过渡带设一个采样点');7.实验结果：程序：clear;N=50; n=0:N-1;k=0:N-1;Hk=zeros(1,5) 0.5 ones(1,4) 0.5 zeros(1,4) 0.5 ones(1,4) 0.5 zeros(1,9) 0.5 ones(1,4) 0.5 zeros(1,4) 0.5 ones(1,4) 0.5 zeros(1,4);h1=real(ifft(exp(-j*pi*(N-1)*k/N-j*pi/2).*Hk);figure;stem(n,h1);xlabel('n');ylabe