北交大数字信号研讨实验2

1、数字信号处理课程研究性学习报告姓名 陶威东 学号 14212129 同组成员 张欣悦 14212138 指导教师 刘留 时间 2015 年6月1日 数字滤波器设计专题研讨【目的】(1) 掌握IIR和FIR数字滤波器的设计方法及各自的特点。(2) 掌握各种窗函数的时频特性及对滤波器设计的影响。(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】 基本题 1分析矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗的频域特性，并进行比较。【题目分析】不同窗函数的主瓣宽度，主瓣高度以及过渡带长度和阻带波动都有明显差别。下面我们将通过对时域长度相同的窗函数进行分析来比较不同窗函数的频

2、域特性。取时域长度为50【仿真结果】【结果分析】各种窗有何特点?通带衰减由小到大的排列为：，；阻带衰减由小到大的排列为：，；过渡带由小到大的排列为：，；在设计滤波器的时候，我们希望得到最大的通带衰减，最小的阻带衰减及最小的过渡带。但在实践中，这三者是不可能同时满足的。所以，我们要根据实际需求，选择合适的窗函数，以达到最佳的滤波效果。最后我们还可以看出，主瓣宽度与主瓣幅值的乘积是定值，即主瓣宽度的增加必将导致主瓣幅值的降低。最后我们会发现凯泽窗幅度响应曲线和矩形窗的完全吻合。这是因为凯泽窗是是一个可以变幻形状的窗函数，当系统没有给值时，默认的为矩形窗，所以是完全符合的。【自主学习内容】几种不同窗

3、函数的设计【阅读文献】1 陈后金.数字信号处理M.北京：高等教育出版社.2008.11【发现问题】 (专题研讨或相关知识点学习中发现的问题)：对于不同窗函数的选择上，我们需要根据需要在主瓣宽度和过渡带长短之间做出选择【问题探究】在谱分析中如何选择窗函数，在滤波器设计中如何选择窗函数？提高FIR滤波器阻带衰减是以增加过渡带宽度为代价的。实际中，在满足阻带衰减的前提下，尽可能地选择主瓣宽度较小的窗函数。凯泽窗是可调的，它通过改变窗函数的形状来控制窗函数旁瓣的大小，设计中可根据滤波器的衰减指标来确定窗函数的形状。【仿真程序】M=50;L=512; fsam=1000;wsam=2*pi*fsam;

4、w=(-wsam/2+(0:L-1)*wsam/L)/fsam;w1=ones(1,M);X1=fftshift(fft(w1,L);plot(w,abs(X1),'g'); hold on;xlabel('Normalized frequency'); legend('square'); title('·幅度响应¦');grid on;w2=hanning(M)' X2=fftshift(fft(w2,L); subplot(2,2,1),plot(w,abs(X2),'r'); ho

5、ld on;xlabel('Normalized frequency'); legend('hanning'); title('·幅度响应');grid on;w3=hamming(M)' X3=fftshift(fft(w3,L); subplot(2,2,2),plot(w,abs(X3),'y'); hold on;xlabel('Normalized frequency'); legend('hamming'); title('幅度响应¦');gr

6、id on;w4=blackman(M)' X4=fftshift(fft(w4,L); subplot(2,2,3),plot(w,abs(X4),'b'); hold on;xlabel('Normalized frequency'); legend('blackman'); title('幅度响应');grid on;w5=kaiser(M)' X5=fftshift(fft(w5,L); subplot(2,2,4),plot(w,abs(X5),'-');xlabel('Norma

7、lized frequency'); legend('kaiser'); title('幅度响应');grid on;【研讨题目】 基本题 2（M5-5）在用窗口法设计FIR滤波器时，由于理想滤波器的频幅响应在截频处发生突变，使得设计出的滤波器的频幅响应发生振荡，这个现象被称为Gibbs现象。解决这个问题的一个方案是本书中介绍的用逐步衰减的窗函数。另一个方案是使理想滤波器过渡带为渐变的，如下图所示具有线性过渡带的理想低通滤波器的频率响应，试用窗口法设计逼近该频率响应的FIR滤波器。 题2图【设计步骤】渐变的窗选为hamming窗。取wp=0.55pi,w

8、s=0.45pi,As=25db,Ap=1db。设hamming窗的长度为M,矩形窗的长度为M1。M1取不同的值作图，其中M=7【单位脉冲响应证明】试证该滤波器的单位脉冲响应为其中：，证明如图【仿真结果】M=8M=16M=32M=64【结果分析】1. 通过逐步衰减的窗函数法和渐变过渡带两种方法设计的滤波器都能较好的设计较小Gibbs现象的滤波器。通过理论分析我们可以知道，FIR滤波器的波动是由于窗函数的主瓣面积和旁瓣面积决定的，所以增加采样点并不能改变通阻带的波动。而对于滤波器的渐变法，则可以通过改变点数很好的减少波动的范围。【自主学习内容】逐步减变法的设计【阅读文献】1 陈后金.数字信号处理

9、M.北京：高等教育出版社.2008.11【发现问题】 (专题研讨或相关知识点学习中发现的问题)：【问题探究】通过实验讨论如何控制滤波器的阻带衰减【仿真程序】wp=0.55*pi;ws=0.45*pi;Ap=1;As=25;N=ceil(7*pi/(wp-ws);N=mod(N+1,2)+N;M=N-1;w=hamming(N);wc=(wp+ws)/2;k=0:M;hd=(wc/pi)*sinc(wc*(k-0.5*M)/pi);h=hd'.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=abs(mag);plot(omeg

10、a/pi,magdb,'b');grid;w=ws-wp;M1=8;k2=-M1:M1;wc=(wp+ws)/2;hd=sinc(w*k2/2).*(sin(wc*k2)./(k2.*pi);hd(M1+1)=wc/pi;omega2=linspace(0,pi,512);mag2=freqz(hd,1,omega2);magdb2=abs(mag2);hold on;plot(omega2/pi,magdb2,'r');legend('逐步衰减','过渡带渐变 ')grid on 【研讨题目】 基本题 3试用频率取样法设计频率响

11、应逼近的FIR数字微分器。分别采用III型IV型线性相位滤波器，并将所得结果进行比较。【题目分析】可利用III型滤波器的频率响应函数和频率取样法设计FIR滤波器【设计步骤】【仿真结果】III型IV型【自主学习内容】【阅读文献】1 陈后金.数字信号处理M.北京：高等教育出版社.2008.11【发现问题】 (专题研讨或相关知识点学习中发现的问题)：在0点处，所设计的滤波器效果并不好，有一定的衰减【问题探究】未发现好的解决办法【仿真程序】III型M=64;Wp1=-pi;Wp2=pi;m=0:M/2;Wm=2*pi*m./(M+1);Ad=double(Wm>=Wp1)&(Wm<

12、=Wp2).*Wm;Hd=j*Ad.*exp(-j*0.5*M*Wm);Hd=Hd conj(fliplr(Hd(2:M/2+1);h=real(ifft(Hd);w=linspace(-pi,pi,1000);H=freqz(h,1,w);plot(w/pi,20*log10(abs(H);axis(-1 1 -50 20);grid;IV型M=63;Wp1=-pi;Wp2=pi;m=0:M/2;Wm=2*pi*m./(M+1);Ad=double(Wm>=Wp1)&(Wm<=Wp2).*Wm;Hd=j*Ad.*exp(-j*0.5*M*Wm);Hd=Hd conj(fl

13、iplr(Hd(2:(M+1)/2);h=real(ifft(Hd);w=linspace(-pi,pi,1000);H=freqz(h,1,w);plot(w/pi,20*log10(abs(H);grid;4（1）（M5-6）利用频率取样法设计III型线性相位带通FIR滤波器，其通带截频分别为Wp1=0.3p rad, Wp2=0.5p rad（2）（M5-7）在通带和阻带间设置1个过渡点，探讨过渡点幅度的确定方法，研究过渡点幅度对设计结果的影响。【设计步骤】（1） 学习利用频率取样法设计某I型线性相位带通FIR滤波器的方法；（2） 确定理想滤波器的幅度函数；（3） 确定理想滤波器的相位；

14、（4） 用实验法确定满足要求的滤波器的阶数；（5） 用实验法确定过度点。【仿真结果】衰减大约为-16.87db(2)取过渡点t1=0.29,t2=0.3,此时可得实验最大衰减最大衰减大约为-52.9【结果分析】滤波器的阻带衰减，滤波器的阶数与设计结果的关系。知滤波器阶数越大，阻带衰减越大，但滤波器阶数越大，实际设计难度越大。【自主学习内容】III型FIR滤波器的设计方法【阅读文献】1 陈后金.数字信号处理M.北京：高等教育出版社.2008.11【发现问题】 (专题研讨或相关知识点学习中发现的问题)：应用实验法寻找过渡点费时费力又不一定得出最好的过渡点【问题探究】如何确定过渡点的最佳幅度。由大及

15、小，不断试验确定【仿真程序】（1）M=64;Wp1=0.3*pi;Wp2=0.5*pi;m=0:M/2;Wm=2*pi*m./(M+1);Ad=double(Wm>=Wp1)&(Wm<=Wp2);Hd=Ad.*exp(-j*0.5*M*Wm);Hd=Hd conj(fliplr(Hd(2:M/2+1);h=real(ifft(Hd);w=linspace(0.1,pi,1000);H=freqz(h,1,w);plot(w/pi,20*log10(abs(H);grid;（2）M=64;Wp1=0.3*pi;Wp2=0.5*pi; Wp=(Wp1+Wp2)/2;m=0:M/

16、2;Wm=2*pi*m./(M+1);mtr1=ceil(Wp1*(M+1)/(2*pi);mtr2=floor(Wp2*(M+1)/(2*pi)+2; Ad=double(Wm>=Wp1)&(Wm<=Wp2);Ad(mtr1)=0.29;Ad(mtr2)=0.30;Hd=Ad.*exp(-j*0.5*M*Wm);Hd=Hd conj(fliplr(Hd(2:M/2+1);h=real(ifft(Hd);w=linspace(0.1,pi,1000);H=freqz(h,1,w);plot(w/pi,20*log10(abs(H);grid;【研讨题目】中等题 5Dhexi

17、an.wav是对频率为293.66, 369.99, 440Hz的D大调和弦以8000Hz抽样所得的数字音乐信号，试设计一数字滤波器从和弦中分离出440Hz的音符。要求：（1）设计IIR数字高通滤波器，通过实验研究，的选择对滤波效果及滤波器阶数的影响，给出滤波器指标选择的基本原则，确定你认为最合适的滤波器指标。（2）能否用IIR数字带通滤波器从和弦中分离出440Hz的音符？利用（1）确定的基本原则，给出数字带通滤波器的指标。设计IIR数字带通滤波器，并将结果与高通滤波器比较，给出你的结论。（3）用窗函数法设计FIR数字高通滤波器，分别利用矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗截断。讨论用窗

18、函数法设计FIR数字高通滤波器时如何确定滤波器的指标，比较相同过渡带时用矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗设计滤波器的阶数。（4）采用Parks-McClellan算法，设计FIR数字高通滤波器。试参照（1）确定的最合适的高通滤波器指标，给出FIR数字高通滤波器的指标。将设计结果与（1）中的IIR数字滤波器，从幅度响应、相位响应、滤波器阶数等方面进行比较。【温磬提示】在IIR数字滤波器的设计中，不管是用双线性变换法还是冲激响应不变法，其中的参数T的取值对设计结果没有影响。但若所设计的数字滤波器要取代指定的模拟滤波器时，则抽样频率（或抽样间隔T）将对设计结果有影响。【设计步骤】（1）首先

19、应先画出原信号的时域和频域图（2）IIR滤波器设计的主要方法是先设计低通模拟滤波器，然后转换为高通、带通或带阻数字滤波器。对于其他如高通，带通，则通过频率变换转换为设计相应的高通，带通等。在设计的全过程的各个步骤，matlab都提供相应的工具箱函数，使得IIR数字滤波器设计变得非常简单。总的来说，我的设计思路主要有以下两种：思路一：从归一化模拟低通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型的数字滤波器。归一化模拟低通原型数字高,带通或带阻模拟高,带通或带阻 模拟域 冲激响应不变法频率变换 双线性变换法图2-1 先频率变换再离散

20、思路二：先进行双线性变换，将模拟低通原型滤波器变换成数字低通滤波器；然后在Z域内经数字频率变换为所需类型的数字滤波器。 归一化模拟低通原型数字高,带通或带阻数字原型低通 数字域 双线性变换法 频率变换图2-2 先离散再频率变换【仿真结果】（1）Wp=0.85*pi; Ws=0.78*pi; Ap=1; As=50;可知此时的滤波器满足设计指标能够进行滤波当改变wp,ws,ap,as时，滤波器的截频会发生一些变化，此时我们可根据实际情况选择大致合适的截频即可Wp=0.82pi,ws=0.79piWp=0.83pi.ws=0.76piWp=0.85*pi; Ws=0.78*pi; Ap=1; As

21、=30;(2)知此时带通滤波器满足设计指标（3）矩形窗N=23HammingN=89 HannN=79BlackmanN=143KasierN=74（4）ds=0.001ds=0.01ds=0.1（4）IIR滤波器相位谱矩形窗函数法相位谱【结果分析】（1）通过仿真可得到，若想取得较好的实验结果，必须选取合适的fp.fs,ap,as值，才能是滤波效果较好。知可选择Wp=0.85pi,Ws=0.78pi,Ap=1,As=50指标进行滤波。（2）知道通过仿真能够进行滤波，对于带通滤波器，选取合适的参数可以完成与高通滤波器同样的滤波效果（3）在相同截取长度下，衰减逐渐增加，同时相应的近似过渡带也逐渐加

22、宽。可以看出在增加衰减的同时会增加过渡带的宽度，所以在设计滤波器时要合理的选择窗函数。而凯泽窗则是完全与实际的衰减幅度类似。可以自主调节所选窗函数的长度。（4）可知，Parks-McClellan算法设计出的滤波器阻带的衰减是等波纹的。且通过改变波动的幅度可以控制衰减的大小。相位响应比较：通过相位响应我们可以看出BW行IIR滤波器的线性性比较好，而其余的IIR滤波器的相位响应都有非线性失真，即相位非线性。而FIR滤波器在通带区间的相位响应都是直线，即线性的。幅度响应比较：通过IIR滤波器和FIR滤波器的幅度响应可以看出，不同方法设计的IIR滤波器的通阻带衰减波动不一样，而FIR滤波器都有较大的

23、波动。两者在衰减幅度方面都能达到滤波器设计的要求。都比较符合设计。阶数方面：通过各个图形的N值可以看出，设计相同效果的滤波器，FIR滤波器的阶数更少。【自主学习内容】【阅读文献】1 陈后金.数字信号处理M.北京：高等教育出版社.2008.11【发现问题】 (专题研讨或相关知识点学习中发现的问题)：对于IIR滤波器的指标选择没有什么比较好的办法，本题中，我才用了逐个实验的方法，可找出一个满足滤波设计指标的参数。【问题探究】对于kasier【仿真程序】（1）Wp=0.85*pi; Ws=0.78*pi; Ap=1; As=50;T=2;Fs=1/T; wp=2*tan(Wp/2)/T;ws=2*t

24、an(Ws/2)/T;wp1=1/wp;ws1=1/ws;N,wc=buttord(wp1,ws1,Ap,As,'s');num,den=butter(N,wc,'s');numa,dena=lp2hp(num,den,1);numd,dend=bilinear(numa,dena,Fs);w=linspace(0,pi,1024);h=freqz(numd,dend,w);plot(w/pi,20*log10(abs(h);axis(0 1 -50 0);grid;xlabel('Normalized frequency');ylabel(&#

25、39;Gain,dB');（2）M=64;Wp1=0.8*pi;Wp2=0.9*pi; Wp=(Wp1+Wp2)/2;m=0:M/2;Wm=2*pi*m./(M+1);mtr1=ceil(Wp1*(M+1)/(2*pi);mtr2=floor(Wp2*(M+1)/(2*pi)+2; Ad=double(Wm>=Wp1)&(Wm<=Wp2);Ad(mtr1)=0.29;Ad(mtr2)=0.30;Hd=Ad.*exp(-j*0.5*M*Wm);Hd=Hd conj(fliplr(Hd(2:M/2+1);h=real(ifft(Hd);w=linspace(0.1,pi

26、,1000);H=freqz(h,1,w);plot(w/pi,20*log10(abs(H);grid;（3）矩形窗clear allFsamp = 10e3;fp = 4.2e3;fs = 3.8e3;Ap = 1;As = 50;Wp = 2*pi*fp/Fsamp; %¹éÒ»»¯Ws = 2*pi*fs/Fsamp;N=ceil(1.8*pi/(Wp-Ws);N=mod(N+1,2)+N;M=N-1;w=ones(1,N);fprintf('N=%.0fn',N);Wc=(Wp+Ws)/2;k=0:M;hd=

27、-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);hd(0.5*M+1)=hd(0.5*M+1)+1;h=hd.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag);plot(omega/pi,magdb); grid;xlabel('Normalized frequency');ylabel('Gain,dB');hann窗wp=0.84*pi;ws=0.76*pi;Ap=1;As=50;N=ceil(6.4*pi/(wp-ws);N=mod(N+1,2)+

28、N;M=N-1;w=hann(N);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)+1;h=hd'.*w;omega=linspace(0,pi,512);mag=freqz(h,1,omega);magdb=20*log10(abs(mag);plot(omega/pi,magdb,'b');grid;w=ws-wp;N对于其他窗函数只需改变窗函数写法即可得到（4）Fp=0.84;Fs=0.76;ds=0.001;dp=ds; f=Fs Fp;a=0 1;dev=ds dp;M,fo,ao,w = remezord(f,a,dev);h = remez(M,fo,ao,w)