数字滤波器的设计实验

1、实验二 IIR数字滤波器的设计实验内容及步骤:数字滤波器的性能指标：通带临界频率fp、阻带临界频率fr；通带内的最大衰减Ap；阻带内的最小衰减Ar；采样周期T；(1)、fp=0.3KHz,Ap=0.8dB, fr=0.2KHz,Ar=20dB,T=1ms；设计一Chebyshev高通滤波器；观察其通带损耗和阻带衰减是否满足要求。程序如下：fp=300; fr=200; Ap=0.8; Ar=20; T=0.001;fs=1/T; wp=2*pi*fp*T;wr=2*pi*fr*T; Wp=2/T*tan(wp/2);Wr=2/T*tan(wr/2);N,Wn=cheb1ord(Wp,Wr,Ap

2、,Ar,'s');B,A = cheby1(N,Ap,Wn,'high','s'); num,den=bilinear(B,A,1/T); h,w=freqz(num,den);plot(w*fs/(2*pi),20*log10(abs(h); %衰减及频率都用归一化的1为单位显示axis(0,500,-30,0);title('Chebyshev高通滤波器');xlabel('频率');ylabel('衰减');grid on;根据下图知道通带损耗与阻带衰减满足要求(2)、fp=0.2KHz,Ap

3、=1dB, fr=0.3KHz,Ar=25dB,T=1ms；分别用脉冲响应不变法及双线性变换法设计一Butterworth数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。程序如下：fp=200; fr=300; Ap=1;Ar=25; T=0.001;fs=1/T; wp=2*pi*fp*T;wr=2*pi*fr*T; Wp=2/T*tan(wp/2);Wr=2/T*tan(wr/2);N,Wn=buttord(Wp,Wr,Ap,Ar,'s');B,A = butter(N,Wn,'s'); num

4、1,den1=impinvar(B,A,1/T); %脉冲响应不变法得出设计的传递函数 num2,den2=bilinear(B,A,1/T); %双线性变换法得出设计的传递函数h1,w=freqz(num1,den1); plot(w*fs/(2*pi),20*log10(abs(h2),w*fs/(2*pi),20*log10(abs(h1),'r.');grid on; %衰减及频率都用归一化的1为单位显示axis(0,500,-30,0);title('Butterworth低通滤波器（红线脉冲响应不变法 蓝线双线性变换法）');xlabel('

5、ƵÂÊ');ylabel('Ë¥¼õ');grid on; 优缺点：采用脉冲响应不变法优点： 1.h(n)完全模仿模拟滤波器的单位抽样响应时域逼近良好 2线性相位模拟滤波器转变为线性相位数字滤波器 缺点： 1.对时域的采样会造成频域的“混叠效应”，故有可能使所设计数字滤波器的频率响应与原来模拟滤波器的频率响应相差很大2不能用来设计高通和带阻滤波器。只适用于限带的低通、带通滤波器采用双线性变换法优点：1避免了频率响应的混迭失真现象 2在特定数字滤波器和特定模拟滤波器处，频率响应是严格相等的，

6、它可以较准确地控制截止频率的位置。 3它是一种简单的代数关系，设计十分方便。 缺点： 1除了零频率附近，w与W之间严重非线性，即线性相位模拟滤波器变为非线性相位数字滤波器 2要求模拟滤波器的幅频响应为分段常数型，不然会产生畸变 3对于分段常数型AF滤波器，经双线性变换后，仍得到幅频特性为分段常数的 DF.但在各个分段边缘的临界频率点产生畸变，这种频率的畸变，可通过频率预畸变加以校正(3)、利用双线性变换法分别设计满足下列指标的Butterworth型和Chebyshev型数字低通滤波器，并作图验证设计结果。fp=1.2kHz, Ap0.5dB, fr=2KHz, Ar40dB, fs=8KHz

7、程序如下：fp=1200;fr=2000; Ap=0.5; Ar=40; fs=8000; T=1/fs; wp=2*pi*fp*T; wr=2*pi*fr*T; Wp=2/T*tan(abs(wp/2); Wr=2/T*tan(abs(wr/2);N,Wn=buttord(Wp,Wr,Ap,Ar,'s');B,A = butter(N,Wn,'s'); num2,den2=bilinear(B,A,1/T); %双线性变换法得到设计的传递函数h2,w=freqz(num2,den2); figure(1);subplot(2,1,1)plot(w/pi*fs/

8、2,10*log10(abs(h2).2); axis(0 fs/2 -50 0); title('Butterworth低通滤波器');xlabel('频率 Hz'); ylabel('衰减 dB'); grid on; subplot(2,1,2)N,Wn=cheb1ord(Wp,Wr,Ap,Ar,'s');B,A = cheby1(N,Ap,Wn,'s'); num,den=bilinear(B,A,1/T); h1,w=freqz(num,den);plot(w/pi*fs/2,10*log10(abs(h

9、1).2); axis(0 fs/2 -50 0); title('chebyshev低通滤波器');xlabel('频率 Hz'); ylabel('衰减 dB');grid on;    (4)、利用双线性变换法设计一Butterworth型数字带通滤波器，已知fs=30KHz，其等效的模拟滤波器指标为Ap3dB, 2KHzf3KHz, Ar5dB, f6KHz, Ar20dB, f1.5KHz  程序如下：f1=2000; f3=3000;fsl=1500; fsh=6000;rp=3; rs=20;

10、Fs=30000;wp1=2*pi*f1/Fs;wp3=2*pi*f3/Fs;wsl=2*pi*fsl/Fs;wsh=2*pi*fsh/Fs;wp=wp1 wp3;ws=wsl wsh;wap=2*Fs*tan(wp./2);was=2*Fs*tan(ws./2);n,wn=buttord(wap,was,rp,rs,'s');z,p,k=buttap(n);bp,ap=zp2tf(z,p,k);bw=wap(2)-wap(1);w0=sqrt(wap(2)*was(1);bs,as=lp2bp(bp,ap,w0,bw);bz1,az1=bilinear(bs,as,Fs);h

11、,w=freqz(bz1,az1,256,Fs);plot(w,20*log10(abs(h);axis(0 8000 -30 0); title('Butterworth带通滤波器');xlabel('频率 Hz');ylabel('衰减 dB');grid on; 1 双线性变换法中和之间的关系是非线性的，在实验中你注意到这种非线性关系了吗？从哪几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？   答：在双线性变换法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟器经过双线性变换后得到的数字滤波

12、器不再保持原有的线性相位了。如以上实验过程中，采用双线性变化法设计的butter和cheby1数字滤波器，从图中可以看到这种非线性关系。  2能否利用公式完成脉冲响应不变法的数字滤波器设计？为什么？  答：不能，这样会使得H(z)中的z以对数形式出现，使H(z)的分子分母不再是z的有理多项式，这样会给系统的分析和实现带来很大的困难。实验三、FIR数字滤波器的设计实验内容及步骤(1) N=15, 。用Hanning窗设计一线性相位带通滤波器，观察它的实际3dB和20dB带宽。N=45，重复这一设计，观察幅频和相位特性的变化，注意长度N变化的影响；程序如下：

13、N1=15;w1=0.3;w2=0.5;wn=w1,w2; b1=fir1(N1-1,wn,hanning(N1); %用Hanning窗作为冲击响应的窗函数figure(1);freqz(b1,1);title('N=15 hanning窗幅频相频特性'); N2=45;b2=fir1(N2-1,wn,hanning(N2);figure(2);freqz(b2,1);title('N=45 hanning窗幅频相频特性');由图可知：随着N值的增大，取样值增大，主瓣的宽度减小，但是幅频特性与相频特性曲线的波动频率会增加。(2) 分别改用矩形窗和Blackma

14、n窗，设计(1)中的带通滤波器，观察并记录窗函数对滤波器幅频特性的影响，比较三种窗的特点；程序如下：N1=15;w1=0.3; w2=0.5; wn=w1,w2; b1=fir1(N1-1,wn,boxcar(N1); %用boxcar作为冲击响应的窗函数figure(1);freqz(b1,1);title('N=15 boxcar窗幅频相频特性'); N2=45;b2=fir1(N2-1,wn,boxcar(N2);figure(2);freqz(b2,1);title('N=45 boxca窗幅频相频特性'); b3=fir1(N1-1,wn,blackm

15、an(N1); %用blackman作为冲击响应的窗函数figure(3);freqz(b3,1);title('N=15 blackman窗幅频相频特性'); N2=45;b4=fir1(N2-1,wn,blackman(N2);figure(4);freqz(b4,1);title('N=45 blackman窗幅频相频特性');分析：由图可知，在三个窗函数中，Blackman窗的衰减性最好，矩形窗的衰减最差。，而矩形窗的过渡带带宽最小，Blackman窗的过渡带带宽最大。Hanning窗的性能比较平均。其实，改善阻带衰减的一种方法是加宽过渡带宽，以牺牲过渡

16、带宽换取阻带衰减的增大。也就是以增加主瓣宽度为代价来降低旁瓣。    (3)用Kaiser窗设计一专用线性相位滤波器，N=40， 如图，当0=4，6，10时，分别设计，比较它们的幅频和相频特性，注意0取不同值时的影响；程序如下：N=40;beta=4,6,10; wn=0.2,0.4,0.6,0.8; s='rgb' for i=1:3 window=kaiser(N,beta(i); b1=fir1(N-1,wn,window); H,w=freqz(b1,1); magH=abs(H);angH=angle(H); subplot(2,1,1)

17、 plot(w/pi, magH,s(i);hold on; grid on xlim(0,1); ylim(0,1); xlabel('w/（2*pi）');ylabel('|H( ejw)|') title('N=40 幅度部分'); legend('=4','=6','=10'); subplot(2,1,2); plot(w/pi,angH,s(i);hold on; grid on; xlim(0,1); xlabel('w/(2*pi)');ylabel('(w)

18、'); title('N=40 相频部分' );end由下图可知，随着的增大，幅频曲线变得越平缓，但其边峰瓣值的衰减越大，相频特性变得越好。(4)、用频率采样法设计(3)中的滤波器，过渡带分别设一个过渡点，令H(k)=0.5。比较两种不同方法的结果；分析：，。而在通带两侧的过渡带需要插入0.5。采用第一种相位特性，因此，可以得到以下式子程序如下：Hr=zeros(1,3),0.5,ones(1,5),0.5,0,0.5,ones(1,5),0.5,zeros(1,5),-0.5,-ones(1,5),-0.5,0,-0.5,-ones(1,5),-0.5,0,0;N =

19、 40;k = 0:N-1; subplot(2,2,1);stem(k,Hr,'.');grid on;xlabel('n'),ylabel('h(n)'),title('插值图'); wm = 2*pi*k./N;Hd = Hr.*exp(-j*(N-1)*wm*0.5);hh = real(ifft(Hd);w1 = linspace(0,pi-0.1,1000);H1 = freqz(hh,1,w1);amgH = abs(H1);dbH = 20*log10(amgH);angH = angle(H1); subplot

20、(2,2,2);plot(w1/pi,dbH),grid on;title('幅度响应');xlabel('w/pi'),ylabel('|H( ejw)|/dB')axis(0 1 -100 20);subplot(2,2,3);plot(w1/pi,angH),grid on;title('相频特性');xlabel('w/pi'),ylabel('')subplot(2,2,4);plot(w1/pi,amgH),grid on;title('幅频特性');xlabel(

21、9;w/pi'),ylabel('|H( ejw)|')（5）、用雷米兹(Remez)交替算法设计(3)中的滤波器，并比较(3)、(4)、(5)三种方法的结果。 程序如下：f=0 0.38 0.42 0.78 0.82 1;A=0 0 1 1 0 0;weigh=1 8 1; b=remez(40,f,A,weigh);%h,w=freqz(b,1,512,'whole');magh=abs(h);dbh=20*log10(magh);phah=angle(h); subplot(2,1,1)plot(w/(2*pi),dbh);grid on

22、;xlabel('w/(2*pi)'),ylabel('|H( ejw)|/dB'),title('幅频特性');axis(0 1 -70 5) subplot(2,1,2)plot(w/(2*pi),phah);grid on;xlabel('w/（2*pi）'),ylabel('(w)'),title('相频特性'); 四、实验思考1. 定性地说明用本实验程序设计的FIR滤波器的3dB截止频率在什么位置？它等于理想频率响应Hd(ej)的截止频率吗？ 答：在之间，大致等于理想频率响应Hd

23、(ej)的截止频率。2. 如果没有给定h(n)的长度N，而是给定了通带边缘截止频率c和阻带临界频率p，以及相应的衰减，你能根据这些条件用窗函数法设计线性相位FIR低通滤波器吗？答：能，可以根据通带边缘截止频率c和阻带临界频率p，以及相应的衰减与所选用的窗函数来大致判断其长度N，从而设计线性相位FIR低通滤波器。实验四、用MATLAB信号处理工具箱进行滤波器设计训练针对一个含有5Hz、50Hz和100Hz的混和正弦波信号，设计一个FIR带通滤波器并对其进行滤波。利用MATLAB实现的三种方法：程序设计法、 FDATool设计法和SPTool设计法。参数要求：采样频率fs=100Hz，通带下限截止

24、频率fc1=10 Hz，通带上限截止频率fc2=20 Hz，过渡带宽6 Hz，通阻带波动0.01，采用凯塞窗设计。方法1 ：程序设计法fsamp = 100;fcuts = 4 10 20 26;mags = 0 1 0;devs = 0.01 0.01 0.01;n,wn,beta,ftype = kaiserord(fcuts,mags,devs,fsamp);n = n + rem(n,2);hh = fir1(n,wn,ftype,kaiser(n+1,beta),'noscale');H,f = freqz(hh,1,1024,fsamp); figure(1)plo

25、t(f,abs(H),grid on;xlabel('f/Hz'),ylabel('|H( ejw)|'),title('kaiser数字带通滤波器');用一个序列测试滤波器性能figure(2)N = 256;n = 0:1:N-1;t =n/fsamp;xn = sin(2*pi*5*t)+sin(2*pi*15*t)+sin(2*pi*35*t); subplot(2,2,1)stem(n,xn,'.');grid on;xlabel('n');ylabel('x(n)');title(

26、9;序列x(n)');axis(0 256 -2 2); X = fft(xn,N);magX = abs(X);f = n*fsamp/N; subplot(2,2,2);plot(f(1:N/2),magX(1:N/2);grid on;%axis(0 fs 0 1);xlabel('f/Hz');ylabel('|X(ejw)|');title('序列x(n)幅频特性'); x = xn;y = filter(hh,1,x);Y = fft(y,N);magY = abs(Y); subplot(2,2,3);stem(n,y,&#

27、39;.');grid on;xlabel('n');ylabel('y(n)');title('滤波后的序列y(n)');axis(0 256 -2 2); subplot(2,2,4);plot(f(1:N/2),magY(1:N/2);grid on;xlabel('f/Hz');ylabel('|Y(ejw)|');title('序列y(n)幅频特性');方法2：利用FDATool设计法在MATLAB命令窗口输入FDATool后回车就会弹出FDATool界面。  带通滤波器设计:   已知滤波器的阶数n=38，beta=3.4。(1) 首先在Filter Type中选择Bandpass；(2) 在Design Method选项中选择FIR Window，(3) 接着在Window选项中选取Kaiser，Beta值为3.4；(4) 指定Filter Order项中的Specify order为38；采样频率Fs=100Hz,截止频率Fc1=10Hz,Fc2=20Hz。(5) 设置完以后点击窗口下方的Design Filter，在窗口上方就会看