数字信号处理研讨

1、数字信号处理课程研究性学习报告数字滤波器设计专题研讨【目的】(1) 掌握IIR和FIR数字滤波器的设计方法及各自的特点。(2) 掌握各种窗函数的时频特性及对滤波器设计的影响。(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。【研讨题目】 基本题 1分析矩形窗、汉纳窗、哈明窗、布莱克曼窗、凯泽窗的频域特性，并进行比较。【题目分析】分析不同的窗函数的频率特性，可以看出其主瓣宽度与旁瓣宽度的差异，过渡带的宽度的差异。【仿真结果】【结果分析】各种窗特点：各种窗函数都采用了相同的长度，计算fft的长度均为512个点。从结果可以看出矩形窗的主瓣幅度最大，但其的宽度最小，并且其旁瓣幅度也

2、较大。其他的几种窗函数恰恰与矩形窗相反，blackman窗与kaiser窗的旁瓣幅度很小，几乎为0。【自主学习内容】1、 各种窗函数的调用；2、 利用fftshift函数处理FFT的结果，便于观察。【阅读文献】 数字信号处理教材及MATLAB教程【发现问题】 (专题研讨或相关知识点学习中发现的问题)：探讨Kaiser窗时bate值的改变可以改变窗函数的形状，从而达到不同的阻带衰减。bate值为0时，Kaiser窗就是矩形窗，并且随着bate值增加，Kaiser窗在两端的衰减逐渐加大。【问题探究】在谱分析中如何选择窗函数，在滤波器设计中如何选择窗函数？答：根据设计要求的p、s和Ap、As确定滤波

3、器的c和窗函数的类型及其长度N,再确定窗函数的幅度函数和相位函数，计算窗函数的IDTFT得到hdk,加窗截断而得到hk。【仿真程序】L=512;N=20;figure(1)w1=zeros(1,100);w2=ones(1,N);w3=zeros(1,100);wh1=w1 w2 w3;WH1=fftshift(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

4、);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(WH4);title('blackman的频域图');figure(5)bate=20;wh5=kaiser(N,bate);WH5=fftshift(fft(

5、wh5,L);w=(0:L-1)-L/2;plot(w,abs(WH5);title('kaiser窗的频域图');【研讨题目】 基本题 2（M5-5）在用窗口法设计FIR滤波器时，由于理想滤波器的频幅响应在截频处发生突变，使得设计出的滤波器的频幅响应发生振荡，这个现象被称为Gibbs现象。解决这个问题的一个方案是本书中介绍的用逐步衰减的窗函数。另一个方案是使理想滤波器过渡带为渐变的，如下图所示具有线性过渡带的理想低通滤波器的频率响应，试用窗口法设计逼近该频率响应的FIR滤波器。 题2图【设计步骤】（1） 根据所需设计的滤波器，确定线性相位滤波器的类型。（2） 确定理想滤波器的

6、幅度函数Ad（）。（3） 确定理想滤波器的相位d（）=-0.5M+。（4） 计算hdk的IDTFT。（5） 用窗函数截断hdk而得到hk。【单位脉冲响应证明】试证该滤波器的单位脉冲响应为其中：，【仿真结果】M1=8M1=16M1=64M1=128【结果分析】FIR滤波器增加采样点，即提高分辨率，对于阻带波动没有效果。而滤波器渐变设计提高分辨率能够减少波动。【自主学习内容】【阅读文献】数字信号处理教材【发现问题】 【问题探究】【仿真程序】wp=0.7*pi;ws=0.3*pi;Ap=1;As=30;N=ceil(7*pi/(wp-ws);N=mod(N+1,2)+NM=N-1;w=hamming

7、(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(omega/pi,magdb,'b');grid;w=ws-wp;M1=1;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,

8、omega2);magdb2=abs(mag2);hold on;plot(omega2/pi,magdb2,'r');legend('逐步','渐变');grid on;【研讨题目】中等题 3Dhexian.wav是对频率为293.66, 369.99, 440Hz的D大调和弦以8000Hz抽样所得的数字音乐信号，试设计一数字滤波器从和弦中分离出440Hz的音符。要求：（1）设计IIR数字高通滤波器，通过实验研究，的选择对滤波效果及滤波器阶数的影响，给出滤波器指标选择的基本原则，确定你认为最合适的滤波器指标。（2）能否用IIR数字带通滤波器从和

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

10、数等方面进行比较。【温磬提示】在IIR数字滤波器的设计中，不管是用双线性变换法还是冲激响应不变法，其中的参数T的取值对设计结果没有影响。但若所设计的数字滤波器要取代指定的模拟滤波器时，则抽样频率（或抽样间隔T）将对设计结果有影响。【设计步骤】 1、首先应先画出原信号的时域和频域图 2、IIR滤波器设计的主要方法是先设计低通模拟滤波器，然后转换为高通、带通或带阻数字滤波器。对于其他如高通，带通，则通过频率变换转换为设计相应的高通，带通等。在设计的全过程的各个步骤，matlab都提供相应的工具箱函数，使得IIR数字滤波器设计变得非常简单。总的来说，我的设计思路主要有以下两种：思路一：从归一化模拟低

11、通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S域变换到Z域，而得到所需类型的数字滤波器。归一化模拟低通原型数字高,带通或带阻模拟高,带通或带阻 模拟域 冲激响应不变法频率变换 双线性变换法图2-1 先频率变换再离散思路二：先进行双线性变换，将模拟低通原型滤波器变换成数字低通滤波器；然后在Z域内经数字频率变换为所需类型的数字滤波器。 归一化模拟低通原型数字高,带通或带阻数字原型低通 数字域 双线性变换法 频率变换图2-2 先离散再频率变换【仿真结果】矩形窗hann窗Hamming窗Blackman窗Kasier窗（4）ds=0.001ds=0.01ds=0.

12、1【结果分析】1、通过仿真可得到，若想取得较好的实验结果，必须选取合适的fp.fs,ap,as值，才能是滤波效果较好。知可选择Wp=0.85pi,Ws=0.78pi,Ap=1,As=50指标进行滤波。 2、知道通过仿真能够进行滤波，对于带通滤波器，选取合适的参数可以完成与高通滤波器同样的滤波效果3、在相同截取长度下，衰减逐渐增加，同时相应的近似过渡带也逐渐加宽。可以看出在增加衰减的同时会增加过渡带的宽度，所以在设计滤波器时要合理的选择窗函数。而凯泽窗则是完全与实际的衰减幅度类似。可以自主调节所选窗函数的长度。4、可知，Parks-McClellan算法设计出的滤波器阻带的衰减是等波纹的。且通过

13、改变波动的幅度可以控制衰减的大小。相位响应比较：通过相位响应我们可以看出BW行IIR滤波器的线性性比较好，而其余的IIR滤波器的相位响应都有非线性失真，即相位非线性。而FIR滤波器在通带区间的相位响应都是直线，即线性的。幅度响应比较：通过IIR滤波器和FIR滤波器的幅度响应可以看出，不同方法设计的IIR滤波器的通阻带衰减波动不一样，而FIR滤波器都有较大的波动。两者在衰减幅度方面都能达到滤波器设计的要求。都比较符合设计。阶数方面：通过各个图形的N值可以看出，设计相同效果的滤波器，FIR滤波器的阶数更少。【自主学习内容】【阅读文献】【发现问题】 (专题研讨或相关知识点学习中发现的问题)：【问题探

14、究】【仿真程序】（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*tan(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*

15、log10(abs(h);axis(0 1 -50 0);grid;xlabel('Normalized frequency');ylabel('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

16、.*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;（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=

17、mod(N+1,2)+N;M=N-1;w=ones(1,N);fprintf('N=%.0fn',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); grid;xlabel('Normalized frequency');ylabel('Gain,

18、dB');N=23hann窗wp=0.84*pi;ws=0.76*pi;Ap=1;As=50;N=ceil(6.4*pi/(wp-ws);N=mod(N+1,2)+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');g

19、rid;w=ws-wp;N=79Hamming窗wp=0.84*pi;ws=0.76*pi;Ap=1;As=50;N=ceil(6.4*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);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,'

20、b');grid;w=ws-wp;N=89 Blackman窗N=143wp=0.84*pi;ws=0.76*pi;Ap=1;As=50;N=ceil(6.4*pi/(wp-ws);N=mod(N+1,2)+N;M=N-1;w=balckman(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(ome

21、ga/pi,magdb,'b');grid;w=ws-wp;Kasier窗N=74wp=0.84*pi;ws=0.76*pi;Ap=1;As=50;N=ceil(6.4*pi/(wp-ws);N=mod(N+1,2)+N;M=N-1;beta=0.1102*(As-8.7);w=kaiser(N,beta);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)

22、;magdb=20*log10(abs(mag);plot(omega/pi,magdb,'b');grid;w=ws-wp;（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);w=linspace(0,pi,1000);mag=freqz(h,1,w);plot(w/pi,20*log10(abs(mag);xlabel('Normalized frequency');ylabel('G

23、ain, dB');grid;相位谱程序Wp=0.84*pi; Ws=0.76*pi; Ap=1; As=50;T=2;Fs=1/T; wp=2*tan(Wp/2)/T;ws=2*tan(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);

24、 plot(w/pi, angle(h); grid;xlabel('Normalized frequency');ylabel('Gain,dB');title(' BWÐÍÏàλÏìÓ¦');Wp=0.84*pi;Ws=0.76*pi;Ap=1;As=50;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=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi);hd(0.5*M+1)=hd(0.5*M+1)+1;h=hd.*