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

1、数字信号处理课程研究性学习报告姓名 初步草稿 学号 同组成员 指导教师 时间 数字滤波器设计专题研讨【目的】(1) 掌握IIR和FIR数字滤波器的设计方法及各自的特点。(2) 掌握各种窗函数的时频特性及对滤波器设计的影响。(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】中等题 Dhexian.wav是对频率为293.66, 369.99, 440Hz的D大调和弦以8000Hz抽样所得的数字音乐信号，试设计一数字滤波器从和弦中分离出440Hz的音符。要求：（1）设计IIR数字高通滤波器，通过实验研究不同、过渡带、对滤波器设计的影响，确定本题最合适的滤波器指标。（

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

3、的IIR数字滤波器，从幅度响应、相位响应、滤波器阶数等方面进行比较。【温磬提示】在IIR数字滤波器的设计中，不管是用双线性变换法还是冲激响应不变法，其中的参数T的取值对设计结果没有影响。但若所设计的数字滤波器要取代指定的模拟滤波器时，则抽样频率（或抽样间隔T）将对设计结果有影响。【设计步骤】:首先以8000Hz抽样时，293.66, 369.99, 440Hz对应的数字频率分别为0.07341*pi,0.0924975*pi,0.11*pi。此处在设计低通滤波器时，Wp显然要小于0.11*pi，而Ws要大于0.0924975*pi。另外此处设计的IIR滤波器均为BW型滤波器，采用双线性变换法。

4、【仿真结果】：考虑到一般滤波器对通阻带截频的要求，我们可以先让As=0.5,Ap=15(此两个指标暂不作考虑）【结果分析】 1.Wp=0.11*pi,Ws=0.092497*pi,Ap=0.5,As=15 音乐滤波后的频谱关键位置放大看评价：听声音感觉滤波前后有一些变化，嘈杂、浑浊的声音变得较为清脆，显然是频谱成分减少了。仔细观察滤波前后音乐的频谱，440Hz的信号保留了，但300400Hz间也有部分信号没有被滤掉，尽管其幅度较小。猜测是阻带截频太小导致此情况，所以增大阻带截频。 改变Ws2. Wp= 0.11*pi Ws=0.0944975 Ap=0.5,As=15 音乐滤波后的频谱关键位置

5、放大看评价：听声音感觉和上次一样，滤波前后有一些变化，嘈杂、浑浊的声音变得较为清脆，显然是频谱成分减少了。仔细观察滤波前后音乐的频谱，440Hz的信号保留了，但300400Hz间也有部分信号没有被滤掉，但其幅度较上次有所减小。3. Wp= 0.11*pi Ws=0.0954975 Ap=0.5,As=15音乐滤波后的频谱关键位置放大看评价：听声音感觉滤波前后有一些变化，嘈杂、浑浊的声音变得较为清脆，显然是频谱成分减少了。仔细观察滤波前后音乐的频谱，440Hz的信号保留了，但300400Hz间的信号基本被滤掉，剩下的部分幅度很小，结果较为满意。4. Wp= 0.11*pi Ws=0.096497

6、5 Ap=0.5,As=15音乐滤波后的频谱关键位置放大看评价：听声音感觉不如前几次效果好，无论从滤波器的幅度响应还是滤波后信号的频谱看，效果均较差。综上，确定Ws=0.0954975pi 改变Wp5. Wp= 0.1099*pi Ws=0.0954975 *pi Ap=0.5,As=15 音乐滤波后的频谱关键位置放大看评价：听声音感觉不如前几次效果好，无论从滤波器的幅度响应还是滤波后信号的频谱看，效果均较差。Wp减小一点就会造成较大误差。 再次总结：取Wp=0.11*pi，Ws=0.0954975*pi 探究Ap、As的影响6. Wp= 0.11*pi Ws=0.0954975*pi Ap=

7、0.5,As=13音乐滤波后的频谱关键位置放大看 关键位置放大看7. Wp= 0.11*pi Ws=0.0954975*pi Ap=0.5,As=14 音乐滤波后的频谱关键位置放大看8. Wp= 0.11*pi Ws=0.0954975*pi Ap=0.5,As=15 音乐滤波后的频谱关键位置放大看 9.Wp= 0.11*pi Ws=0.0954975*pi Ap=0.5,As=16 关键位置放大看再次总结：从声音的变化、滤波器的幅度响应和滤波后的频谱除440Hz其他成分的幅度三方面来看，确定As=1510.Wp= 0.11*pi Ws=0.0954975*pi Ap=0.3,As=1511.

8、 Wp= 0.11*pi Ws=0.0954975*pi Ap=1.0,As=15 音乐滤波后的频谱关键位置放大看12. Wp= 0.11*pi Ws=0.0954975*pi Ap=1.5,As=15 音乐滤波后的频谱关键位置放大看综上：确定高通滤波器的指标为Wp= 0.11*pi（440Hz） Ws=0.0954975*pi （381.99Hz） Ap=0.5,As=15(2) 从理论上分析，只要440Hz在数字带通滤波器的通带范围内，293.66Hz,369.99Hz处在其阻带内，就能从和弦中分离出440Hz的音符。根据以上设计的高通滤波器的指标先确定带通滤波器的以下指标，W1.11*p

9、i,Ws1=0.0954975*pi,As=15,Ap=0.5,至于Wp2和Ws2的确定如下：1. Wp2=0.190*pi,Ws2=0.999*pi实验程序：fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%听滤波前的音乐f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音乐频谱');%设计数字滤波器Ws=0.0954975;W

10、p=0.11;Ap=0.5;As=15;N,Wc=buttord(Wp,Ws,Ap,As)num,den=butter(N,Wc,'high');%设计数字滤波器w=linspace(0,pi,1024);h=freqz(num,den,w);%求滤波器的增益响应figure(2);plot(w/pi,20*log10(abs(h);gridtitle('高通数字滤波器的增益响应');y=filter(num,den,x);%滤波操作wavplay(y,fs);%听滤波后的音乐Y=fftshift(fft(y,L);figure(3);f=(0:L-1)*fs/

11、L-0.5*fs;plot(f,abs(Y);title('音乐滤波后的频谱'); 音乐滤波后的频谱关键位置放大看2. Wp2= 0.19*pi,Ws2=0.8*pi音乐滤波后的频谱关键位置放大看3. Wp2=0.19 *pi,Ws2=0.5*pi 音乐滤波后的频谱关键位置放大看4. Wp2= 0.4*pi,Ws2=0.5*pi5.Wp2=0.4 *pi,Ws2=0.7*pi6. Wp2=0.89 *pi,Ws2=0.9075025*pi总结分析：由以上实验结果可知不同的Wp2,Ws2指标对实验结果的影响很小，究其原因是所给的音乐的频谱较为单一，它只在293.66Hz,369.

12、99Hz和440Hz及其附近幅度较大，在大于440Hz处幅度基本为零，所以带通滤波器的Wp2,Ws2指标对实验结果的影响非常小。带通滤波器确实能从和弦中分离出440Hz的音符，在此处从滤波器指标的对称性考虑，我们选择带通滤波器的指标为Ws1=0.0954975*pi(381.99Hz);Wp1=0.11*pi(440Hz);Wp2=0.89*pi(3560Hz);Ws2=0.9045025*pi(3619Hz);As=15;Ap=15。实验程序：fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian

13、.wav');wavplay(x,fs);%听滤波前的音乐f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音乐频谱');%设计数字滤波器Ws1=0.0954975;Wp1=0.11;Wp2=0.89;Ws2=0.9075025;Wp=Wp1,Wp2;Ws=Ws1,Ws2;Ap=0.5;As=15;N,Wc=buttord(Wp,Ws,Ap,As)num,den=butter(N,Wc);%设计数字滤波器w=linspace(0,pi,1024);h=freqz(num,

14、den,w);%求滤波器的增益响应figure(2);plot(w/pi,20*log10(abs(h);gridtitle('高通数字滤波器的增益响应');y=filter(num,den,x);%滤波操作wavplay(y,fs);%听滤波后的音乐Y=fftshift(fft(y,L);figure(3);f=(0:L-1)*fs/L-0.5*fs;plot(f,abs(Y);title('音乐滤波后的频谱');（3）1. 用矩形窗（阶数M=124) 音乐滤波后的频谱关键位置放大看%用矩形窗fs=8000;nbits=16;L=8192;x,fs,nbits

15、=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%听滤波前的音乐f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音乐频谱');N=ceil(1.8*pi/(Wp-Ws);N=mod(N+1,2)+N;M=N-1w=linspace(1,N,1);Wc=(Wp+Ws)/2;k=0:M;hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);hd(0.5*M+1)=hd(0.5*M+1)

16、+1;h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag);figure(2);plot(omega/pi,magdb);title('用矩形窗设计的高通滤波器的幅度响应');y=filter(h,1,x);wavplay(y,fs);%听滤波后音乐的频谱f=(0:L-1)*fs/L-0.5*fs;Y=fftshift(fft(y,L);figure(3);plot(f,abs(Y);title('滤波后音乐的频谱');2. 用Hanning窗(阶数M=428）

17、音乐滤波后的频谱关键位置放大看汉宁窗法的程序：%用hanning窗fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%听滤波前的音乐f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音乐频谱');Ws=0.0954972*pi;Wp=0.11*pi;Wc=(Wp+Ws)/(2*pi);%使滤波器为I型N=ceil(6.2*pi/(W

18、p-Ws);N=mod(N+1,2)+N;M=N-1h = fir1(M,Wc,'high',hanning(N);omega=linspace(0,pi,512);mag=freqz(h,1,omega);figure(2);plot(omega/pi,20*log10(abs(mag);title('用hanning窗设计的高通滤波器的幅度响应');y=filter(h,1,x);wavplay(y,fs);%听滤波后音乐的频谱f=(0:L-1)*fs/L-0.5*fs;Y=fftshift(fft(y,L);figure(3);plot(f,abs(Y);

19、title('滤波后音乐的频谱');3. 用哈明窗（阶数M=482) 音乐滤波后的频谱关键位置放大看哈明窗法的程序：%用Hamming窗fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%听滤波前的音乐f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音乐频谱');Ws=0.0954972*pi;Wp=0.11*pi

20、;Wc=(Wp+Ws)/(2*pi);%使滤波器为I型N=ceil(7*pi/(Wp-Ws);N=mod(N+1,2)+N;M=N-1h = fir1(M,Wc,'high',hamming(N);omega=linspace(0,pi,512);mag=freqz(h,1,omega);figure(2);plot(omega/pi,20*log10(abs(mag);title('用Hamming窗设计的高通滤波器的幅度响应');y=filter(h,1,x);wavplay(y,fs);%听滤波后音乐的频谱f=(0:L-1)*fs/L-0.5*fs;Y=f

21、ftshift(fft(y,L);figure(3);plot(f,abs(Y);title('滤波后音乐的频谱');4. 布莱克曼窗(阶数M=786） 音乐滤波后的频谱关键位置放大看Blackman窗法的程序：%用blackman窗fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%听滤波前的音乐f=(0:L-1)*fs/L-0.5*fs;X=fftshift(fft(x,L);figure(1);plot(f,abs(X);t

22、itle('原始音乐频谱');Ws=0.0954972*pi;Wp=0.11*pi;Wc=(Wp+Ws)/(2*pi);%使滤波器为I型N=ceil(11.4*pi/(Wp-Ws);N=mod(N+1,2)+N;M=N-1h = fir1(M,Wc,'high',blackman(N);omega=linspace(0,pi,512);mag=freqz(h,1,omega);figure(2);plot(omega/pi,20*log10(abs(mag);title('用blackman窗设计的高通滤波器的幅度响应');y=filter(h,

23、1,x);wavplay(y,fs);%听滤波后音乐的频谱f=(0:L-1)*fs/L-0.5*fs;Y=fftshift(fft(y,L);figure(3);plot(f,abs(Y);title('滤波后音乐的频谱');5. 凯泽窗（阶数M=166）音乐滤波后的频谱关键位置放大看采用凯泽窗的程序：%用kaiser窗fs=8000;nbits=16;L=8192;x,fs,nbits=wavread('C:UserslenovoDesktopDhexian.wav');wavplay(x,fs);%听滤波前的音乐f=(0:L-1)*fs/L-0.5*fs;X

24、=fftshift(fft(x,L);figure(1);plot(f,abs(X);title('原始音乐频谱');%Program:利用Kaiser窗设计FIR高通滤波器Ap=0.5;As=15;Rp=1-10.(-0.05*Ap);Rs=10.(-0.05*As);f=0.0954975,0.11;a=0,1;dev=Rp,Rs;M,Wc,beta,ftype = kaiserord(f,a,dev);%使滤波器为I型M=mod(M,2)+Mh = fir1(M,Wc,ftype,kaiser(M+1,beta);omega=linspace(0,pi,512);mag=

25、freqz(h,1,omega);figure(2);plot(omega/pi,20*log10(abs(mag);title('用kaiser窗设计的高通滤波器的幅度响应');y=filter(h,1,x);wavplay(y,fs);%听滤波后音乐的频谱f=(0:L-1)*fs/L-0.5*fs;Y=fftshift(fft(y,L);figure(3);plot(f,abs(Y);title('滤波后音乐的频谱');【结果分析】：一、用窗函数法设计FIR滤波器的方法及步骤：1.根据Wp,Ws，Ap,As确定逼近的理想滤波器及窗函数类型及窗的长度。由Wp,

26、Ws确定理想滤波器Wc及窗长度N,由As确定窗类型，2. 根据要求确定线性相位FIR滤波器的类型3. 确定理想滤波器的幅度函数Ad()和相位d()，d()= -0.5M+ 4.计算IDTFT得hdk 5.截短hdk, hk= hdkwk， 0<=k<= M二、用MATLAB实现的方法： h = fir1(M,Wc,'ftype',window)M表示滤波器的阶数；Wc表示理想FIR滤波器的B个频带； ftype表示滤波器的类型，缺省值为空；window是一长度为N+1的向量, Hamming窗3、 本题中相同过渡带时，矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗的阶数

27、分别是 矩形窗（124）<凯泽窗（166）<汉纳窗（428）<哈明窗（482）<布莱克曼窗（786）(4) 采用Parks-McClellan算法 滤波后音乐频谱关键位置放大看【结果分析】：比较采用Parks-McClellan算法和IIR数字滤波器。从幅度响应来看，前者在通带和阻带都会波动，而后者只在阻带波动，通带不波动从相位响应来看，前者在0到pi间波动，且较均匀，而后者也在0到pi间波动，但是在低频附近波动更厉害，随着频率升高，波动减弱。从阶数来看，前者为102阶，后者为20阶，IIR数字滤波器的阶数明显低于采用Parks-McClellan算法设计的滤波器的阶数。%采用Parks-McClellan算法fs=8000;nbits=16;L=8192;x,fs,nbits=wavread(&#