北交大DSP研究性学习报告数字滤波器设计2014

1、数字信号处理课程研究性学习报告姓名 学号 同组成员 指导教师 时间 2014年5月18日星期日 数字滤波器设计专题研讨【目的】(1) 掌握IIR和FIR数字滤波器的设计方法及各自的特点。(2) 掌握各种窗函数的时频特性及对滤波器设计的影响。(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】 基本题 分析矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗的频域特性，并进行比较。【题目分析】本题分析不同的窗函数的频域特性，预计可以看出不同的窗有不同的过渡带大小和不同的旁瓣宽度，可以满足对不同的设计设计要求。【仿真结果】【结果分析】各种窗有何特点?计算过程中，不同的窗都采用了

2、相同的长度（均为10），进行fft计算的长度均为512点，从以上结果可以看出，矩形窗的幅度最大，主瓣最窄，旁瓣幅度大；其余几个窗函数的主瓣幅度均小于矩形窗的主瓣幅度，主瓣宽度约为矩形窗主瓣宽度的2倍多；其中汉纳窗和哈明窗的旁瓣有较小幅度，而布莱克曼窗和凯泽窗旁瓣基本为0。【阅读文献】1数字信号处理.陈后金【发现问题】 (专题研讨或相关知识点学习中发现的问题)：matlab提供的函数产生的窗为N行1列的列向量，不是1行N列的一行数据，计算时必须将其进行转置。使用fft计算的结果是原序列周期化后求出的频谱的主值区间，必须使用fftshift将其整理，然后通过向左平移L/2才能求出原序列的频谱。【问

3、题探究】在谱分析中如何选择窗函数，在滤波器设计中如何选择窗函数？在谱分析中窗函数的选择需要考虑泄露现象，当需要分辨的频率为f时，若使用矩形窗，其宽度必须大于fsamf。当要求分辨较小的频率而计算量较小可以选择其他窗函数来改善。在FIR滤波器的设计中，使用矩形窗会由于其锐截止而产生Gibbs现象，会有9%的过冲，为了改善效果，可以在满足过渡带宽度和增益要求的前提下适当选择其他种类的窗函数。【仿真程序】%clear all;L=512;%L2N=10;%figure(1)w1=zeros(1,50);w2=ones(1,N);w3=zeros(1,50);wh1=w1 w2 w3;WH1=ffts

4、hift(fft(wh1,L);w=(0:L-1)-L/2;plot(w,abs(WH1);title();%figure(2)wh2=hann(N);WH2=fftshift(fft(wh2,L);w=(0:L-1)-L/2;plot(w,abs(WH2);title();%figure(3)wh3=hamming(N);WH3=fftshift(fft(wh3,L);w=(0:L-1)-L/2;plot(w,abs(WH3);title();%figure(4)wh4=blackman(N);WH4=fftshift(fft(wh4,L);w=(0:L-1)-L/2;plot(w,abs(

5、WH4);title();%figure(5)beta=6;wh5=kaiser(N,beta);WH5=fftshift(fft(wh5,L);w=(0:L-1)-L/2;plot(w,abs(WH5);title();【研讨题目】中等题 Dhexian.wav是对频率为293.66, 369.99, 440Hz的D大调和弦以8000Hz抽样所得的数字音乐信号，试设计数字滤波器从和弦中分离出369.99Hz的音符。要求：（1）设计IIR数字带通滤波器，通过实验研究不同、过渡带、对滤波器设计的影响，确定本题最合适的滤波器指标。（2）通过将IIR数字低通滤波器和数字高通滤波器级联也可实现369.

6、99Hz音符的分离。分别确定 IIR数字低通滤波器和数字高通滤波器的指标，并进行滤波器设计。将结果与带通滤波器比较，给出你的结论。（3）用窗函数法设计FIR数字带通滤波器，分别利用矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗截断。讨论用窗函数法设计FIR数字带通滤波器时如何确定滤波器的指标，比较相同过渡带时用矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗设计滤波器的阶数。（4）采用Parks-McClellan算法，设计FIR数字带通滤波器。试参照（1）确定的最合适的带通滤波器指标，给出FIR数字带通滤波器的指标。将设计结果与（1）中的IIR数字滤波器，从幅度响应、相位响应、滤波器阶数等方面进行比较

7、。【温磬提示】在IIR数字滤波器的设计中，不管是用双线性变换法还是冲激响应不变法，其中的参数T的取值对设计结果没有影响。但若所设计的数字滤波器要取代指定的模拟滤波器时，则抽样频率（或抽样间隔T）将对设计结果有影响。【设计步骤】1.讲所给模拟频率转换为数字指标：通过公式=2ffsam可以计算出三个音符对应的数字频率1=0.073 2=0.092 3=0.11 2.采用双线性变换法，通过公式=2Ttan(2)将数字频率指标转换成模拟指标3.根据模拟指标设计模拟滤波器4.利用公式Hz=Hs|s=2T1-z-11+z-1将所得模拟滤波器转换为数字滤波器【仿真结果】(1) 确定数字带通滤波器的参数：首先

8、讨论通阻带截频参数，考虑到一般滤波器对通阻带截频的要求，我们可以先让As=1,Ap=20(此两个指标暂不作考虑）【1】取p1=0.089 ,p2=0.095 ,s1=0.073, s2=0.11 从滤波后频谱可以看出，仍留存有高低频分量，因此应调整阻带频率【2】取p1=0.089 ,p2=0.095 ,s1=0.080, s2=0.10 滤波效果较好，经多次调整后，最终确定的滤波器参数为p1=0.085 ,p2=0.099 ,s1=0.078, s2=0.105。 下面讨论Ap, As两个指标为变量研究其对低通滤波器的影响一、研究As对滤波器的影响 As10dB As20dB As32dB A

9、s40dB由上面四幅图可知，一味增大阻带衰减幅度可能是结果更糟，只需选取适当阻带衰减幅度即可，故令As20dB 二、研究Ap对滤波器的影响 Ap0.5dB Ap1dBAp2dB由图可知Ap0.5dB时阻带衰减幅度最大，但Ap1dB时也可以满足滤波的要求，考虑到系统成本，选取Ap1dB综上所述，带通滤波器的指标确定为：p1=0.085 ,p2=0.099 ,s1=0.078, s2=0.105,Ap1dB, As20dB (2) 指标： Wp2=0.095p, Ws2=0.1p 指标：Wp1=0.0875p, , Ws1=0.085p, (3)有表格可以看出，选取不同的窗函数，要通过确定其过渡带

10、宽度来确定其阶数，因此，在相同过渡带宽度的前提下，我们通过不同的窗函数设计了不同的带通滤波器，其幅度相应及阶数如下图所示。矩形窗 N=301汉纳窗 N=1034哈明窗 N=1167布莱克曼窗 N=1901凯泽窗 N=237(beta=3.395)(4)(1)中带通滤波器的指标如下：Wp1=0.085p, Wp2=0.099p, Ws1=0.078p, Ws2=0.105p Ap1dB As20dB 转换为FIR数字带通滤波器指标Wp1=0.085p, Wp2=0.099p, Ws1=0.078p, Ws2=0.105p dp = 0.109 ds = 0.1 等纹滤波器幅度谱和相位谱 N=54

11、5IIR带通滤波器 N=6【结果分析】(1)(2)采用高低通滤波器级联的方法也能达到滤波的效果，同时成本较带通而言基本一致，但不易将杂波完全滤去，仍会保留一些高频、低频分量。从频谱来看，滤波效果不如带通滤波器；从声音效果来看，两者几乎一样。高低通级联滤波效果不如用带通滤波器好，但在物理上更易实现；在误差允许的范围内，高低通级联也可以实现带通的效果。因此，工程上可以考虑使用高低通级联代替带通滤波器。(3)用窗函数发设计FIR数字带通滤波器时，首先确定过渡带宽度，然后要选择窗函数，根据窗函数近似过渡带的宽度确定滤波器的阶数N。由实验结果可知，按照矩形窗、汉纳窗、哈明窗、布莱克曼窗的顺序，阻带衰减逐

12、渐增大，阶数逐渐增高。而选择凯泽窗时，既能保证较高的阻带衰减，又有较小的阶数，与其他窗函数相比，更为灵活。(4)比较采用Parks-McClellan算法和IIR数字滤波器。 幅度响应来看，前者在通带和阻带都会波动，而后者波纹明显较小，但过渡带较前者宽； 从相位响应来看，前者在0到pi间波动，且较均匀，而后者也在0到pi间波动，但是在低频附近波动更厉害，随着频率升高，波动减弱。 从阶数来看，前者为545阶，后者为6阶，IIR数字滤波器的阶数明显低于采用Parks-McClellan算法设计的滤波器的阶数。【发现问题】 (专题研讨或相关知识点学习中发现的问题)：设计滤波器的过程中使用BW型有时频

13、谱会有振荡产生,并且当阻带增益As越大时，震荡越明显，甚至会出现严重的失真。【问题探究】过渡带越窄，通带最大衰减越小，阻带最小衰减越大，则滤波器设计难度越大；反之，过渡带越宽，通带最大衰减越大，阻带最小衰减越小，滤波器更易设计。我们发现，当滤波器设计难度到达一定程度时，matlab已经不能设计出相应滤波器或是出现上图震荡的情况。经查阅资料，发现应该是滤波器指标不合适导致达到计算机计算的极限，结果会有相当大的误差。原本理论可实现的系统计算机无法算出。但这也说明这样的系统物理上也很难实现。所以当发现这种情况时，我们应该改变设计指标，使滤波器更易实现。【仿真程序】(1)Fsam=8000;x1 =

14、audioread(Dhexian.wav);y1,f1=ctft1(x1,Fsam,8001);t=0:8000;figure(1)plot(f1,abs(y1)axis(200,500,0,0.18)sound(x1,Fsam);pause(2); Ws1=0.078*pi; Wp1=0.085*pi; Wp2=0.099*pi ; Ws2=0.105*pi; Ap=1;As=20;%ws1=tan(Ws1/2); wp1=tan(Wp1/2); wp2=tan(Wp2/2); ws2=tan(Ws2/2);%B=wp2-wp1; w0=sqrt(wp1*wp2); wp=1; ws=mi

15、n(abs(ws12-wp1*wp2)/(B*ws1),abs(ws22-wp1*wp2)/(B*ws2);%N,Wc=buttord(wp,ws,Ap,As,s);num,den = butter(N,Wc,s);numt,dent = lp2bp(num,den,w0,B);numd,dend=bilinear(numt,dent,0.5);w=linspace(0,pi,1024);h=freqz(numd,dend,w);figure(2)plot(w/pi,20*log10(abs(h),r);grid;xlabel(Normalized frequency);ylabel(Gain

16、,dB);axis(0,0.2,-100,0)%hold on; figure(3)x2=filter(numd,dend,x1);y2,f2=ctft1(x2,Fsam,8001);plot(f2,abs(y2);axis(-500,500,0,0.18)player1=audioplayer(x2,fs);play(player1);axis(200,500,0,0.18) (2)fs=8000;nbits=16;L=8192;x,fs,nbits=wavread(E:DSP2Dhexian.wav);wavplay(x,fs);f=(0:L-1)*fs/L-0.5*fs;X=fftshi

17、ft(fft(x,L);figure(1);plot(f,abs(X);title();%hpWp1=365*2/fs;Ws1=330*2/fs;Wp2=3560*2/fs;Ws2=3605*2/fsWp=Wp1,Wp2;Ws=Ws1,Ws2;Ap=1;As=12;N,wc=buttord(Wp,Ws,Ap,As);num,den=butter(N,wc);w=linspace(0,pi,1000);h1=freqz(num,den,w);figure(2);plot(w/pi,20*log10(abs(h1);grid on;title();axis(0 0.8 -180 20);y1=fi

18、lter(num,den,x);%?wavplay(y1,fs);%?Y1=fftshift(fft(y1,L);figure(3);f=(0:L-1)*fs/L-0.5*fs;plot(f,abs(Y1);title();%lp Wp4=370*2/fs;Ws4=400*2/fsAp=1;As=7;N,wc=buttord(Wp4,Ws4,Ap,As);num,den=butter(N,wc);w=linspace(0,pi,1000);h2=freqz(num,den,w);figure(4);plot(w/pi,20*log10(abs(h2);grid on;title();y2=fi

19、lter(num,den,y1);%?wavplay(y2,fs);%?Y2=fftshift(fft(y2,L);figure(5);f=(0:L-1)*fs/L-0.5*fs;plot(f,abs(Y2);title();h=h1.*h2;%plot(w/pi,20*log10(abs(h);grid on;(3)Ws1=0.078*pi; Wp1=0.085*pi; Wp2=0.099*pi ; Ws2=0.105*pi; Ap=1;As=20;Wc1=(Wp1+Ws1)/2;Wc2=(Wp2+Ws2)/2;% N1=ceil(1.8*pi/(Wp1-Ws1);% N2=ceil(1.8

20、*pi/(Ws2-Wp2);% N=max(N1 N2);% numd = fir1(N,Wc1 Wc2,boxcar(N+1);% N1=ceil(6.2*pi/(Wp1-Ws1);% N2=ceil(6.2*pi/(Ws2-Wp2);% N=max(N1 N2);% numd = fir1(N,Wc1 Wc2,hann(N+1);% N1=ceil(7*pi/(Wp1-Ws1);% N2=ceil(7*pi/(Ws2-Wp2);% N=max(N1 N2);% numd = fir1(N,Wc1 Wc2);% N1=ceil(11.4*pi/(Wp1-Ws1);% N2=ceil(11.4

21、*pi/(Ws2-Wp2);% N=max(N1 N2);% numd = fir1(N,Wc1 Wc2,blackman(N+1);%f=Ws1 Wp1 Wp2 Ws2;a=0,1,0;Rs=0.01;dev=Rs*ones(1,length(a);N,Wc,beta,ftype = kaiserord(f,a,dev); numd = fir1(N,Wc,ftype,kaiser(N+1,beta);w=linspace(0,pi,1024);h=freqz(numd,1,w);plot(w/(pi*pi),20*log10(abs(h),r);grid;xlabel(Normalized frequency);ylabel(Gain,dB);axis(0,0.2,-100,0)(4)clear all;fs=8000;Ap=0.5;As=50;ds=0.1;dp=0.109;fs1