实验四___数字滤波器设计及应用综合实验

实验四数字滤波器设计及应用综合实验班级：电子121班姓名：李金恩学号：122114一、实验目的1．熟悉IIR数字滤波器的设计原理及方法。2．熟悉FIR数字滤波器的设计原理及方法。3.掌握利用Matlab实现数字滤波器的方法4.掌握利用数字滤波器进行信号处理的方法。5.了解基于Simulink的动态仿真实现信号滤波的基本方法。二、实验内容及要求综合运用数字滤波器设计的相关知识，根据给定设计方法要求，用脉冲响应不变法和双线性变换法设计IIR数字滤波器；利用窗函数设计法设计FIR数字滤波器。根据实际信号的频谱特性，分析、确定滤波器设计技术指标，实现对信号的滤波。1．IIR数字滤波器设计（1）用脉冲响应不变法设计巴特沃斯数字滤波器。（2）用双线性变换法设计切比雪夫数字滤波器。（3）用双线性变换法设计巴特沃斯数字滤波器。并将直接型结构转换成级联型结构。实验所需程序及函数%butterworth低通滤波器原型设计函数，要求Ws﹥Wp﹥0，As﹥Rp﹥0。function[b,a]=afd_butt(Wp,Ws,Rp,As)；N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws)))；%上条语句为求滤波器阶数，函数ceil朝正无穷大方向取整；fprintf('\nButterworthFilterOrder=%2.0f\n',N)；OmegaC=Wp/((10^(Rp/10)-1)^(1/(2*N)))；%求对应于N的3db截止频率；[b,a]=u_buttap(N,OmegaC)；%非归一化Butterworth模拟低通滤波器原形设计函数%得到的b,a分别为传输函数分子、分母多项式系数；function[b,a]=u_buttap(N,Omegac)；[z,p,k]=buttap(N)；%归一化巴特沃思模拟低通滤波器原形%传输函数用极点形式表示()((1))((2))...(())kHsspspspnp=p*Omegac;%将css/代入上式，相当于分子乘以Nc，极点乘以ck=k*Omegac^N;B=real(poly(z));%poly为构造具有指定根的多项式real为求实部b=k*B;a=real(poly(p));%利用脉冲响应不变法从模拟到数字滤波器变换函数function[b,a]=imp_invr(c,d,T)[R,p,k]=residue(c,d);%部分分式展开p=exp(p*T);%从模拟到数字极点对应关系sTez,部分分式系数相同[b,a]=residuez(R,p,k);%将部分分式的形式变换成多项式之比的形式b=real(b');%求出数字滤波器系数a=real(a');%非归一化切比雪夫I型模拟低通滤波器原型设计function[b,a]=u_chb1ap(N,Rp,Omegac)[z,p,k]=cheb1ap(N,Rp);%归一化切比雪夫1型模拟低通滤波器原形a=real(poly(p));%以下步骤实际上与求巴特沃思滤波器的原理aNn=a(N+1);%一样，只是所用方法稍有不同。p=p*Omegac;a=real(poly(p));aNu=a(N+1);k=k*aNu/aNn;B=real(poly(z));b=k*B;%频率响应函数freqz的修正，此函数可获得滤波器的幅值响应、相位响应及群延迟响应function[db,mag,pha,w]=freqz_m(b,a);[H,w]=freqz(b,a,1000,'whole');%在0-2*pi之间选取N个点计算频率响应H=(H(1:501))'；%频率响应w=(w(1:501))';%频率mag=abs(H);%响应幅度db=20*log10((mag+eps)/max(mag));%增益pha=angle(H);%相位%变直接形式为级联形式22,11,22,11,01111011...1...)(zAzAzBzBbzazazbzbbzHkkkkkNNNNfunction[b0,B,A]=dir2cas(b,a)b0=b(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0;%以上步骤求出系数0bM=length(b);N=length(a);ifNMb=[bzeros(1,N-M)];elseifMNa=[azeros(1,M-N)];elseNM=0;endK=floor(N/2);B=zeros(K,3);A=zeros(K,3);ifK*2==Nb=[b0];a=[a0];endbroots=cplxpair(roots(b));%以下程序将每两个极点和两个零点组合成二阶因子aroots=cplxpair(roots(a));%roots：求多项式的根fori=1:2:2*KBrow=broots(i:1:i+1,:);Brow=real(poly(Brow));B(fix(i+1)/2,:)=Brow;Arow=aroots(i:1:i+1,:);Arow=real(poly(Arow));A(fix(i+1)/2,:)=Arow;End2.FIR数字滤波器设计用窗函数法设计FIR数字低通滤波器（窗函数选择两种）。要求求出数字滤波器的阶数，并画出频率特性。FIR数字滤波器设计指标为wp=0.2*pi；Rp=1db；ws=0.3*pi；As=20db。①该程序由学生自己设计。②产生理想低通滤波器的冲激响应函数functionhd=ideal_lp(wc,M)，见本实验“具体实验内容说明→有关函数说明”。3.利用数字滤波器实现DTMF信号的提取双音多频（DualToneMultiFrequency）信号是音频电话中的拨号信号，每一个数字（0-9）由两个不同频率单音组成（每个单音用正弦信号表示），所用频率分为高频带和低频带两组，每个数字由高、低频带中各一个频率组成，例如数字9使用852Hz和1477Hz两个频率。数字与符号对应频率关系见表3-1所示。实验涉及到的程序：functionhd=ideal_lp(wc,M)alpha=(M-1)/2;%21Mn=[0:(M-1)];m=n-alpha+eps;%nmeps是一个非常小的数，防止m为零hd=sin(wc*m)./(pi*m);%mpimhdc*)*sin(表3-1双频拨号的频率分配1209Hz1336Hz1477Hz1633Hz697Hz123A770Hz456B852Hz789C941Hz*0#D设计要求：（1）通过查阅资料，了解电话中DTMF信号的产生与检测方法。（2）自己选择一个数字（0-9），通过数字方法产生该数字的双频信号；设采样频率为8000Hz。（3）根据该信号的频率特性，确定滤波器类型及技术指标，将两个单音分别提取出来。（4）画出数字滤波器输入、输出信号的波形及频谱。（5）要求自拟实验方案实现。三、实验步骤IIR数字滤波器设计：已知：IIR数字滤波器设计指标为wp=0.2*pi；ws=0.3*pi；Rp=1db；As=15db。（1）用脉冲响应不变法设计巴特沃斯数字滤波器实验程序：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;OmegaP=wp/T;OmegaS=ws/T;[cs,ds]=afd_butt(OmegaP,OmegaS,Rp,As);[b,a]=imp_invr(cs,ds,T)[db,mag,pha,w]=freqz_m(b,a);subplot(2,1,1);plot(w/pi,mag);title('digitalfilterMagnitudeResponse')axis([0,1,0,1.1])subplot(2,1,2);plot(w/pi,db);title('digitalfilterMagnitudeinDB')axis([0,1,-40,5]);实验结果：N=6ButterworthFilterOrder=6OmegaC=0.7032b=0.1209a=1.00002.71703.69103.17881.82520.66440.1209b=0.00000.00060.01010.01610.00410.0001a=1.0000-3.36355.0684-4.27592.1066-0.57060.0661（2）用双线性变换法设计切比雪夫数字滤波器。实验程序：wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2);ep=sqrt(10^(Rp/10)-1);Ripple=sqrt(1/(1+ep*ep));Attn=1/(10^(As/20));A1=1/Attn;a1=sqrt(A1*A1-1)/ep;a2=OmegaS/OmegaP;N=ceil(logm(a1+sqrt(a1*a1-1))/logm(a2+sqrt(a2*a2-1)));fprintf('\nChebyshevFilterOrder=%2.0f\n',N)[cs,ds]=u_chb1ap(N,Rp,OmegaP);[b,a]=bilinear(cs,ds,1/T)[db,mag,pha,w]=freqz_m(b,a);subplot(2,1,1);plot(w/pi,mag);title('digitalfilterMagnitudeResponse');axis([0,1,0,1.1])subplot(2,1,2);plot(w/pi,db);title('digitalfilterMagnitudeinDB');axis([0,1,-40,5]);实验结果：ChebyshevFilterOrder=4b=0.00180.00730.01100.00730.0018a=1.0000-3.05433.8290-2.29250.5507（3）用双线性变换法设计巴特沃斯数字滤波器。并将直接型结构转换成级联型结构。wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2);[cs,ds]=afd_butt(OmegaP,OmegaS,Rp,As);[b,a]=bilinear(cs,ds,T)[db,mag,pha,w]=freqz_m(b,a);subplot(2,1,1);plot(w/pi,mag);title('digitalfilterMagnitudeResponse');axis([0,1,0,1.1])subplot(2,1,2);plot(w/pi,db);title('digitalfilterMagnitudeinDB');axis([0,1,-40,5]);[b0,B,A]=dir2cas(b,a)实验结果：N=6ButterworthFilterOrder=6OmegaC=0.7273b=0.1480a=1.00002.81003.94823.51682.08840.78620.1480b=0.00060.00350.00870.01160.00870.00350.0006a=1.0000-3.31434.9501-4.14332.0275-0.54580.0628b0=5.7969e-004B=1.00002.03351.03381.00001.99961.00001.00001.96690.9673A=1.0000-0.94590.23421.0000-1.05410.37531.0000-1.31430.71492.FIR数字滤波器设计：本实验所需M文件如下：functionhd=ideal_lp(wc,M);alpha=(M-1)/2;n=[0:(M-1)];m=n-alpha+eps;hd=sin(wc*m)./(pi*m);本实验程序如下：f