实验四 IIR数字滤波器设计

实验四 IIR数字滤波器设计一、实验目的1、熟悉用双线性变换法设计IIR数字滤波器的原理与方法；2、熟悉用脉冲响应不变法设计IIR数字滤波器的原理与方法。

二、实验内容1、已知低通滤波器的指标为：通带边缘频率：0.4π，Ap=0.5dB；阻带边缘频率：0.6π，As=50dB；a、采用脉冲响应不变法设计巴特沃斯，T=1.画出幅度响应和数字滤波器的脉冲响应h（n）；b、b、采用脉冲响应不变法设计巴特沃斯，T=1.画出幅度响应和数字滤波器的脉冲响应h（n）。

2、用双线性变换法设计低通滤波器，满足技术指标wp=0.2π，Ap=0.25dB；ws=0.4π，As=50dB，并对方波信号进行滤波，画出滤波前后的波形图并进行简要分析3、设计一个数字高通滤波器H（z），它用在结构xa（t）A/D H（z) D/A ya（t）中，满足下列要求：a、采样速率为10kHZ；b、阻带边缘频率为1.5kHZ，衰减为40dB；c、通带边缘频率为2kHZ，衰减为3dB；d、单调的通带和阻带。

4、设计一个带阻滤波器，要求通带上下截止频率为0.8π，0.2π，通带内衰减不大于1dB，阻带起始频率为0.7π，0.4π，阻带内衰减不小于30dB。

设计巴特沃斯带阻滤波器并画出该数字高通滤波器的幅度响应和脉冲响应。

三、实验程序及解释和实验分析及图形1a、clear;close all; %清屏wp=0.4*pi;ws=0.6*pi;Ap=0.5As=50;T=1;Fs=1/T;OmegaP=wp/T;OmegaS=ws/T;[cs,ds]=afd_butt(OmegaP,OmegaS,Ap,As);[b,a]=impinvar(cs,ds,Fs);[h,w]=freqz(b,a);subplot(2,2,1);plot(w/pi,abs(h));title('幅度响应');grid;subplot(2,2,2);plot(w/pi,angle(h));title('相位响应');grid;subplot(2,2,2);plot(w/pi,angle(h));title('相位响应');grid;subplot(2,2,3);plot(w/pi,20*log(abs(h)));title('幅度响应dB');grid;n=[0:1:59];imp=[1;zeros(59,1)];y=filter(b,a,imp);subplot(2,2,4);plot(n,y);title('脉冲响应');grid;分析：因为w=Ωt,w与Ω呈线性关系，所以其相位响应图是线性的，如相位响应图所示；因为其设计的是低通滤波器，所以会把其高频部分滤掉，留下低频，从而不是特别陡峭，如脉冲响应图所示。

1b、clear;close all;wp=0.4*pi; %通带边缘频率设为0.4*piws=0.6*pi; %阻带边缘频率设为0.6*piAp=0.5; %通带最大衰减为0.5dBAs=50; %阻带最小衰减为50dBT=1; %把周期T设为1Fs=1/T;OmegaP=(2/T)*tan(wp/2); %计算通带截止频率原型OmegaS=(2/T)*tan(ws/2); %计算阻带截止频率原型[cs,ds]=afd_butt(OmegaP,OmegaS,Ap,As);%巴特沃斯原型滤波器计算[b,a]=bilinear(cs,ds,Fs); %双线性不变法变换[h,w]=freqz(b,a); %计算数字滤波器的Z变换频率响应subplot(2,2,1);plot(w/pi,abs(h)); %画幅频响应曲线图title('幅度响应');grid;subplot(2,2,2);plot(w/pi,angle(h)); %画出相频响应曲线图title('相位响应');grid;subplot(2,2,3);plot(w/pi,20*log(abs(h)));%画出取其对数的幅频响应曲线title('幅度响应dB');axis([0,1,-80,5]);grid; %设置坐标轴，加网格线n=[0:1:59]; %设n值imp=[1;zeros(59,1)]; %输入矩形波y=filter(b,a,imp); %对imp进行滤波subplot(2,2,4);plot(n,y); %画出y的图形title('脉冲响应');grid;分析：用双线性变换法实现低通巴特沃斯滤波器时，由Ω=(2/T)tan(w/2),可以知道w与Ω是呈非线性关系的，所以如“相位响应”图所示，呈现的是不是直线，而是曲线。

由于其设计的是低通滤波器，所以把其高频的部分滤掉了，留下低频的部分，从而在变化的时候表现的平缓些，而不是特别的陡峭，如“脉冲响应”图所示，从而实现低通滤波的功能。

2、wp=0.2*pi; %通带边缘频率设为0.2*piws=0.4*pi; %阻带边缘频率设为0.4*piAp=0.25; %通带最大衰减为0.25dBAs=50; %阻带最小衰减为50dB%原型指标的频率逆映射T=1; %把周期T设为1Fs=1/T;OmegaP=(2/T)*tan(wp/2); %计算通带截止频率原型OmegaS=(2/T)*tan(ws/2); %计算阻带截止频率原型N=ceil((log10((10^(Ap/10)-1)/(10^(As/10)-1)))/(2*log10(OmegaP/OmegaS))); %确定N值fprintf('\n***Butterworth Filter Order =%2.0f\n',N);OmegaC=OmegaP/((10^(Ap/10)-1)^(1/(2*N)));wn=2*atan((OmegaC*T)/2); %计算截止频率wn=wn/pi;[b,a]=butter(N,wn); %计算截止频率为wn的N阶低通巴特沃斯滤波器[h,w]=freqz(b,a); %计算数字滤波器的Z变换频率响应subplot(2,2,1);plot(w/pi,abs(h)); %画幅频响应曲线图title('幅度响应');grid ;subplot(2,2,3);plot(w/pi,20*log10(abs(h))); %画出取其对数的幅频响应曲线title('幅度响应dB');axis([0,1,-80,5]);grid ; %设置坐标轴，加网格线n=[0:1:59];imp=[1;zeros(59,1)];y=filter(b,a,imp); %对imp进行滤波subplot(2,2,4);plot(n,y); %画出y的图形title('脉冲响应');grid ;n=[0:100];t=0.001*n; %fs=1/t=1000HZx=2*sin(2*pi*400*t)+5*sin(2*pi*50*t);%输入含有两种频率的波形x1=ones(1,101); %输入单位抽样序列x1y1=filter(b,a,x1); %对x1进行滤波figure(2);subplot(2,1,1);plot(n,x1,'p-'); %绘制一条由实线串起来的五角星形线x1subplot(2,1,2);plot(n,y1,'p-'); %绘制y1曲线3、f=10;fp=2;fs=1.5;wp=2*pi*fp*(1/f); %设置通带边缘频率ws=2*pi*fs*(1/f); %设置阻带边缘频率Ap=3;As=40;[N,wn]=buttord(wp/pi,ws/pi,Ap,As); %确定N值和截止频率wn[b,a]=butter(N,wn,'high'); %设计截止频率为wn的高通滤波器[h,w]=freqz(b,a); %计算数字滤波器的Z变换频率响应subplot(2,2,1);plot(w/pi,abs(h)); %画幅频响应曲线图title('幅度响应');grid ;subplot(2,2,3);plot(w/pi,20*log10(abs(h)));%画幅频的对数的幅频曲线title('幅度响应dB');grid ;n=[0:1:59];imp=[1;zeros(59,1)];y=filter(b,a,imp); %对imp进行滤波subplot(2,2,4);plot(n,y); %画出y的图形title('脉冲响应');grid ;n=[0:1000];t=0.0001*n; %fs=1/t=10000HZx=2*sin(2*pi*4000*t)+5*sin(2*pi*50*t);%输入混有两种频率的波形y1=filter(b,a,x); %对x进行滤波figure(2);subplot(2,1,1);plot(n,x,'p-'); %绘制一条由实线串起来的五角星形xsubplot(2,1,2);plot(n,y1,'p-'); %绘制一条由实线串起来的五角星形线y14、clear;close all;wp=[0.2*pi,0.8*pi];ws=[0.4*pi,0.7*pi];Ap=1;As=30;[N,wn]=buttord(wp/pi,ws/pi,Ap,As); %确定N和截止频率wn的值[b,a]=butter(N,wn,'stop'); %设计截止频率为wn的带阻滤器[h,w]=freqz(b,a); %计算数字滤波器的Z变换频率响应subplot(2,2,1);plot(w/pi,abs(h)); %画幅频响应曲线图title('幅度响应'); %加图名grid;subplot(2,2,3);plot(w/pi,20*log(abs(h))); %画幅频的对数的幅频曲线title('幅度响应dB'); %加图名axis([0,1,-80,5]); %设置坐标grid; %加网格线n=[0:1:59];imp=[1;zeros(59,1)];y=filter(b,a,imp); %对imp进行滤波subplot(2,2,4);plot(n,y); %画出滤波后的曲线ytitle('脉冲响应'); %加图名grid;n=[0:100];t=0.001*n; %fs=1/t=1000HZ x=2*sin(2*pi*400*t)+5*sin(2*pi*50*t);x1=ones(1,101);y=filter(b,a,x); %对x进行滤波，得到yfigure(2);subplot(2,1,1);plot(n,x,'p-'); %绘制一条由实线串起来的五角星形线xsubplot(2,1,2);plot(n,y,'p-'); %绘制一条由实线串起来的五角星形线yfunction [b,a] =afd_butt(Omegp,Omegs,Ap,As);if Omegp<= 0error('Passband edge must be larger than 0')endif Omegs<=Omegperror('Stopband edge must be larger than Passband edge')endif (Ap <= 0) | (As < 0)error('PB ripple and/or SB attenuation must be larger than 0')endN=ceil((log10((10^(Ap/10)-1)/(10^(As/10)-1)))/(2*log10(Omegp/Omegs)) );fprintf('\n*** Butterworth Filter Order = %2.0f \n',N)OmegaC = Omegp/((10^(Ap/10)-1)^(1/(2*N)));[b,a]=u_buttap(N,OmegaC）；Endfunction [b,a] = u_buttap(N,Omegac);[z,p,k] = buttap(N);p = p*Omegac;k = k*Omegac^N;B = real(poly(z));b = k*B;b0= k;a = real(poly(p));end实验分析：2用双线性法设计的低通滤波器，fs=1000HZ，当输入为方波时，把其高频部分滤去，留下低频部分，从而实现低通滤波功能；3设计的为高通数字滤波器，fs=10000HZ，而输入信号为x=x1+x2，x1,x2是由两种不同频率的sin 函数混合而成，可以通过高频，滤去低频，从而实现高频滤波器功能；4设计的为带阻滤波器，fs=1000HZ，输入信号为x=x1+x2，x1,x2由两种不同的频率混合而成，带阻滤波器的幅度响应是关于π对称的，从图中可以看出，没经过滤波前的波形与滤波后波形是一样的，即两种波形没有被阻隔。