matlab数字滤波器设计程序
使用MATLAB进行数字滤波器设计的步骤与方法
使用MATLAB进行数字滤波器设计的步骤与方法数字滤波器是用于信号处理的重要工具,它可以对信号进行去噪、频率调整等操作。
而MATLAB作为一种强大的数学计算软件,提供了丰富的数字信号处理工具箱,可以方便地进行数字滤波器的设计与仿真。
本文将介绍使用MATLAB进行数字滤波器设计的步骤与方法。
1. 了解数字滤波器的基本原理在进行数字滤波器设计之前,首先需要了解数字滤波器的基本原理。
数字滤波器根据其频率响应特性可以分为低通、高通、带通和带阻滤波器等。
此外,数字滤波器的设计还需要考虑滤波器的阶数、截止频率以及滤波器类型等因素。
在设计中,我们可以选择滤波器的类型和相应的参考模型,然后利用MATLAB工具箱提供的函数进行设计。
2. 导入MATLAB中的数字信号处理工具箱使用MATLAB进行数字滤波器设计需要先导入数字信号处理工具箱。
通过在MATLAB命令窗口输入`>> toolbox`即可打开工具箱窗口,并可以选择数字信号处理工具箱进行加载。
加载完成后,就可以调用其中的函数进行数字滤波器设计。
3. 设计数字滤波器在MATLAB中,常用的数字滤波器设计函数有`fir1`、`fir2`、`iirnotch`等。
这些函数可以根据系统特性需求设计相应的数字滤波器。
以FIR滤波器为例,可以使用`fir1`函数进行设计。
该函数需要输入滤波器的阶数和截止频率等参数,输出设计好的滤波器系数。
4. 评估滤波器性能设计好数字滤波器后,需要进行性能评估。
可以使用MATLAB提供的`fvtool`函数绘制滤波器的幅频响应、相频响应和群延迟等。
通过观察滤波器在频域的性能表现,可以判断设计的滤波器是否满足要求。
5. 对滤波器进行仿真在对滤波器性能进行评估之后,还可以使用MATLAB进行滤波器的仿真。
通过将需要滤波的信号输入设计好的滤波器中,观察输出信号的变化,可以验证滤波器的去噪效果和频率调整能力。
MATLAB提供了函数`filter`用于对信号进行滤波处理。
用MATLAB设计FIR数字滤波器
实验八 用MATLAB 设计FIR 数字滤波器(二)一、实验目旳:1、加深对窗函数法设计FIR 数字滤波器旳基本原理旳理解。
2、学习用MATLAB 语言旳窗函数法编写设计FIR 数字滤波器旳程序。
3、理解MATLAB 语言有关窗函数法设计FIR 数字滤波器旳常用函数用法。
二、实验原理:1、用窗函数法设计FIR 数字滤波器 FIR 数字滤波器旳系统函数为N-1-n n=0H(z)=h(n)z ∑这个公式也可以当作是离散LSI 系统旳系统函数M-m -1-2-mmm=0012m N -1-2-k-k12k k k=1bz b +b z +b z ++b z Y(z)b(z)H(z)====X(z)a(z)1+a z +a z ++a z1+a z ∑∑ 分母a 0为1,其他a k 全都为0时旳一种特例。
由于极点所有集中在零点,稳定和线性相位特性是FIR 滤波器旳突出长处,因此在实际中广泛使用。
FIR 滤波器旳设计任务是选择有限长度旳h(n),使传播函数H(e j ω)满足技术规定。
重要设计措施有窗函数法、频率采样法和切比雪夫等波纹逼近法等。
本实验重要简介窗函数法。
用窗函数法设计FIR 数字滤波器旳基本环节如下:(1)根据过渡带和阻带衰减指标选择窗函数旳类型,估算滤波器旳阶数N 。
(2)由数字滤波器旳抱负频率响应H(e j ω)求出其单位脉冲响应h d (n)。
可用自定义函数ideal_lp实现抱负数字低通滤波器单位脉冲响应旳求解。
程序清单如下:function hd=ideal_lp(wc,N) %点0到N-1之间旳抱负脉冲响应%wc=截止频率(弧度)%N=抱负滤波器旳长度tao=(N-1)/2;n=[0:(N-1)];m=n-tao+eps; %加一种小数以避免0作除数hd=sin(wc*m)./(pi*m);其他选频滤波器可以由低通频响特性合成。
如一种通带在ωc1~ωc2之间旳带通滤波器在给定N值旳条件下,可以用下列程序实现:Hd=ideal_lp(wc2,N)-ideal_lp(wc1,N)(3)计算数字滤波器旳单位冲激响应h(n)=w(n)h d(n)。
数字滤波器matlab的程序
数字滤波器matlab的源代码function lvbo(Ua,Ub,choise)%参考指令:lvbo(2*pi,10*pi,1/0/-1)U1=min(Ua,Ub);U2=max(Ua,Ub);Us=16*U2;T=2*pi/Us;T_sum=4*max(2*pi/Ua,2*pi/Ub);sum=T_sum/T;t=T:T:T_sum;x=sin(U1*t)+0.8*sin(U2*t);X=DFT(x);figure(1); subplot(221)U=Us/sum:Us/sum:Us;stem(U,abs(X));grid onaxis([Us/sum,Us/2,0,1.2*max(abs(X))])title('原模拟信号采样频谱图')Ucd=U1+(U2-U1)*1/5;Usd=U2-(U2-U1)*1/5;switch choisecase 1Hz_ejw=IIR_DF_BW(Ucd,1,Usd,30,T,sum);case -1Hz_ejw=IIR_DF_CF(Ucd,1,Usd,30,T,sum);case 0Hz_ejw=FIR_DF_HM(U1,U2,T,sum);otherwiseHz_ejw=IIR_DF_BW(Ucd,1,Usd,30,T,sum);endY=X.*Hz_ejw;y=1/sum*conj(DFT(conj(Y)));figure(1); subplot(224)plot(t,real(y)); title('模拟信号滤波后');grid on axis([0,T_sum,-max(real(y))*1.5,max(real(y))*1.5]) subplot(222);plot(t,x); hold onaxis([0,T_sum,-max(x)*1.2,max(x)*1.2])x=sin(U1*t);plot(t,x,':r');grid ontitle('模拟信号滤波前')function Hz_ejw=IIR_DF_BW(Ucd,Ap,Usd,As,t,sum)% 巴特沃思滤波器E=(10^(0.1*Ap)-1)^0.5;V=(10^(0.1*As)-1)^0.5;Wc=Ucd*t; Ws=Usd*t;Ucd=Wc/t; Usd=Ws/t;Uca=(2/t)*tan(Ucd*t/2); Usa=(2/t)*tan(Usd*t/2);N=ceil(log10(V/E)/log10(Usa/Uca));k=[1:2*N];Spk=exp(j*(pi/2+(2*k-1)/(2*N)*pi));i=find(real(Spk)<0);Sk(1:N)=Spk(i);den=real(poly(Sk'));k0=polyval(den,0);disp('模拟巴特沃思滤波器的归一化统函数 Ha(s) 为')tf(k0,den)syms s z T;den_jU=1;s=s/Uca;for i=1:Nden_jU=s^(N-i+1)*den(i)+den_jU;endHa_s=simple(1/den_jU);H_z=subs(Ha_s,'s',(2/T)*((1-1/z)/(1+1/z)));k=1:sum;w=(2*pi/sum)*k;ejw=exp(j*w);Hz_ejw=subs(H_z,{z,T},{ejw,t*ones(1,length(ejw))}); figure(1); subplot(223)plot(w,abs(Hz_ejw)); grid ontitle('巴特沃思低通滤波器')axis([2*pi/sum,pi,-0.2,1.2*max(abs(Hz_ejw))]) function Hz_ejw=IIR_DF_CF(Ucd,Ap,Usd,As,t,sum)% 切比雪夫低通滤波器E=(10^(0.1*Ap)-1)^0.5;V=(10^(0.1*As)-1)^0.5;Wc=Ucd*t; Ws=Usd*t;Ucd=Wc/t; Usd=Ws/t;Uca=(2/t)*tan(Ucd*t/2); Usa=(2/t)*tan(Usd*t/2);N=ceil(acosh(V/E)/acosh(Usa/Uca));;A=1/E+(1/E^2+1)^0.5;a=1/2*(A^(1/N)-A^(-1/N));b=1/2*(A^(1/N)+A^(-1/N));k=1:2*N;Spk=-a*sin((2*k-1)/(2*N)*pi)+j*b*...cos((2*k-1)/(2*N)*pi);i=find(real(Spk)<0);Sk(1:N)=Spk(i);den=real(poly(Sk'));k0=1;disp('模拟切比雪夫低通滤波器的归一化统函数 Ha(s) 为') tf(k0,den)if (rem(N,2)==1)for i=1:Nk0=k0*(-Sk(i));endelseif ((rem(N,2))==0)k0=1;for i=1:Nk0=k0*(-Sk(i));endendif (rem(N,2)==0)k0=10^(-0.05*Ap)*k0;endk0=real(k0);syms s z T;den_jU=1;s=s/Uca;for i=1:Nden_jU=s^(N-i+1)*den(i)+den_jU;endHa_s=simple(1/den_jU);H_z=subs(Ha_s,'s',(2/T)*((1-1/z)/(1+1/z)));k=1:sum;w=(2*pi/sum)*k;ejw=exp(j*w);Hz_ejw=subs(H_z,{z,T},{ejw,t*ones(1,length(ejw))}); figure(1); subplot(223)plot(w,abs(Hz_ejw));grid ontitle('切比雪夫低通滤波器')axis([2*pi/sum,pi,-0.5,max(abs(Hz_ejw))])function Hz_ejw=FIR_DF_HM(U1,U2,T,sum)wp=U1*T;ws=U2*T;kuan=ws-wp;M=sum;n=[0:1:M-1];wc=(ws+wp)/2;hd=H_D(wc,M);window=hamming_m(M);h_z=hd.*window;Hz_ejw=DFT(h_z);k=1:sum;w=(2*pi/sum)*k;figure(1); subplot(223)plot(w,abs(Hz_ejw));grid onaxis([2*pi/sum,pi,-0.2,1.2*max(abs(Hz_ejw))]);title('海明窗函数低通滤波器')function hd=H_D(wc,N)M=(N-1)/2;for k=-M:Mif k==0hd(k+M+1)=wc/pi;elsehd(k+M+1)=sin(wc*k)/(pi*k);endendfunction wn=hamming_m(M)n=0:M-1;wn(n+1)=0.54-0.46*cos((2*pi*n)/(M-1));function Xk=DFT(xn)% 离散傅立叶变换,xn为原序列,Xk为DFT变换后的序列N=length(xn);n=0:N-1;k=0:N-1;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk;Xk=xn*WNnk;。
matlabfir滤波器设计
matlabfir滤波器设计MATLAB是一个高级编程语言和交互式环境,被广泛应用于各种科学和工程问题的数值分析、数据可视化和编程开发等领域。
FIR滤波器是数字信号处理中经常使用的一种滤波器,它是基于有限长冲激响应的滤波器。
在MATLAB平台上,我们可以使用fir1函数来设计FIR滤波器。
一、FIR滤波器设计基础1.1 什么是FIR滤波器FIR滤波器是有限长冲激响应滤波器,由于其具有线性相位特性和可控阶数等优点,在数字信号处理中得到了广泛的应用。
一般来说,FIR滤波器的频率响应特性由滤波器的系数函数确定。
FIR滤波器的设计一般采用窗函数法、最小二乘法、频率抽取法等方法。
窗函数法是最常见的一种方法,大部分情况下选择的是矩形窗、汉宁窗、布莱克曼窗等。
1.3 fir1函数介绍fir1函数是MATLAB中用于FIR滤波器设计的函数,用法为:h = fir1(N, Wn, type)N为滤波器的阶数,Wn是用于指定滤波器截止频率的参数,type指定滤波器类型,可以是低通、高通、带通、带阻等。
二、使用fir1函数设计FIR滤波器2.1 设计要求采样率为300Hz;滤波器阶数为50;截止频率为50Hz。
2.2 实现步骤(1)计算规范化截止频率规范化截止频率是指在数字滤波器设计中使用的无单位量,通常范围为0到1。
在本例中,我们需要将50Hz的截止频率转化为规范化截止频率。
Wn = 2*50/300 = 1/3根据计算出的规范化截止频率和滤波器阶数,我们可以使用fir1函数来进行滤波器设计。
此处滤波器的阶数为50,规范化截止频率为1/3,类型为低通。
(3)绘制滤波器的幅频响应图为了验证设计的低通FIR滤波器是否符合要求,我们需要绘制其幅频响应图。
freqz(h,1,1024,300)经过上述步骤后,我们就得到了一张低通FIR滤波器的幅频响应图,如下图所示:图1.低通FIR滤波器的幅频响应图三、总结通过上述例子,我们可以看出在MATLAB中与fir1函数可以非常方便的进行FIR滤波器的设计。
基于MATLAB的数字滤波器设计
编号淮安信息职业技术学院毕业论文学生姓名。
学号。
系部。
专业。
班级。
指导教师。
顾问教师。
摘要本论文介绍了FIR数字滤波器的设计方法,即窗函数法。
在此基础上,用MATLAB实现IIR数字滤波器。
介绍了IIR数字滤波器的传统设计思想与步骤,及其计算机辅助设计方法。
以一数字带通滤波器为例,着重说明了基于MATLAB的三种滤波器的实现手段:模拟低通原型、合适模拟带通及直接原型,为数字滤波器设计带来全新的实现手段。
关键词:滤波 IIR滤波器 FIR滤波器MATLAB淮安信息职业技术学院目录第一章前言1.1 MATLAB 软件简介 (4)1.2数字滤波器技术的发展状况 (5)第二章数字滤波器的基本概念2.1数字滤波器的概况 (6)2.2 FIR 数字滤波器的基本概念 (6)2.2.1 FIR 数字滤波器的窗函数设计法 (7)2.2.2 窗函数设计法的步骤 (7)2.3 MATLAB环境下的实例 (9)2.3.1高通滤波器的设计 (9)2.3.2低通滤波器的设计 (10)第三章 IIR数字滤波器的设计过程及方法3.1 IIR滤波器的基本特点 (13)3.2 IIR滤波器的设计思路与步骤 (14)3.3 IIR 滤波器的设计 (14)3.4 IIR滤波器设计方法MATLAB的实现 (15)3.4.1 基于模拟低通原型的MATLAB实现 (15)3.4.2基于合适类型模拟滤波器的MATLAB实现 (16)3.4.3 基于直接原型变换法的MATLAB实现 (18)总结 (19)参考文献 (20)第一章前言1.1. MATLAB简介MATLAB (Matrix Laboratory)为美国Mathworks公司1983年首次推出的一套高性能的数值分析和计算软件,其功能不断扩充,版本不断升级,1992年推出划时代的4.0版,1993年推出了可以配合Microsoft Windous使用的微机版,95年4.2版,97年5.0版,99年5.3版,5.X版无论是界面还是内容都有长足的进展,其帮助信息采用超文本格式和PDF格式,可以方便的浏览。
matlab频率采样法设计滤波器
一、介绍频率采样法设计滤波器的背景和意义1.1 频率采样法设计滤波器的概念及其在数字信号处理中的作用 1.2 频率采样法设计滤波器与其他设计方法的比较1.3 频率采样法设计滤波器的优势和适用范围二、频率采样法设计滤波器的原理和方法2.1 频率采样法设计滤波器的基本原理2.2 频率采样法设计滤波器的设计步骤2.3 频率采样法设计滤波器的常用工具和软件三、matlab频率采样法设计滤波器的实现步骤3.1 设定滤波器的规格和要求3.2 使用matlab进行频域设计3.3 使用matlab进行时域设计3.4 验证设计的滤波器性能四、matlab频率采样法设计滤波器的案例分析4.1 案例一:低通滤波器设计4.1.1 滤波器规格要求4.1.2 频率采样法设计滤波器的实现步骤4.1.3 设计参数及性能分析4.2 案例二:带通滤波器设计4.2.1 滤波器规格要求4.2.2 频率采样法设计滤波器的实现步骤4.2.3 设计参数及性能分析五、matlab频率采样法设计滤波器的应用前景和挑战5.1 应用前景分析5.2 技术发展趋势5.3 面临的挑战和解决方案六、总结与展望6.1 频率采样法设计滤波器的优势和不足6.2 matlab工具在频率采样法设计滤波器中的应用6.3 未来发展方向和趋势在数字信号处理中,滤波器设计是一项重要的工作。
频率采样法设计滤波器是其中一种常用的设计方法,在matlab软件中进行频率采样法设计滤波器具有高效、便捷的特点。
本文将介绍频率采样法设计滤波器的原理、方法以及在matlab中的实现步骤,通过案例分析和应用前景展望来全面解析这一设计方法的优势和发展趋势。
在数字信号处理领域,滤波器设计是至关重要的一环。
而频率采样法设计滤波器作为一种常用的设计方法,在matlab软件中具有高效、便捷的特点。
接下来,我们将深入探讨频率采样法设计滤波器的原理、方法以及在matlab中的实现步骤,并通过案例分析和应用前景展望来全面解析这一设计方法的优势和发展趋势。
数字滤波器的设计及其MATLAB实现
设计低通数字滤波器,要求在通带内频率低于0.2pirad时,允许幅度误差在1dB以内,在频率0.3pi rad~pi rad之间的阻带衰减大于15dB,用脉冲响应不变法设计数字滤波器,T=1: 切比雪夫I型模拟滤波器的设计子程序:function [b,a]=afd_chb1(Omegap,Omegar,Ar)if Omegap<=0error('通带边缘必须大于0')endif(Dt<=0)|(Ar<0)error('通带波动或阻带衰减必须大于0');endep=sqrt(10^(Dt/10)-1);A=10^(Ar/20);OmegaC=Omegap;OmegaR=Omegar/Omegap;g=sqrt(A*A-1)/ep;N=ceil(log10(g+sqrt(g*g-1))/log10(OmegaR+sqrt(OmegaR*OmegaR-1)));fprintf('\n***切比雪夫I型模拟低通滤波器阶数=%2.0f\n',N);[b,a]=u_chblap(N,Dt,OmegaC);设计非归一化切比雪夫I型模拟低通滤波器原型程序:function [b,a]=u_chblap(N,Dt,OmegaC)[z,p,k]=cheb1ap(N,Dt);a=real(poly(p));aNn=a(N+1);p=p*OmegaC;a=real(poly(p));aNu=a(N+1);k=k*aNu/aNn;b0=k;B=real(poly(z));b=k*B;直接形式转换成级联形式子程序:function [C,B,A]=sdir2cas(b,a)Na=length(a)-1;Nb=length(b)-1;b0=b(1);b=b/b0;a0=a(1);a=a/a0;C=b0/a0;p=cplxpair(roots(a));K=floor(Na/2);if K*2==NaA=zeros(K,3);for n=1:2:NaArow=p(n:1:n+1,:);Arow=poly(Arow);A((fix(n+1)/2),:)=real(Arow);elseif Na==1A=[0 real(poly(p))];elseA=zeros(K+1,3);for n=1:2:2*KArow=p(n:1:n+1,:);Arow=poly(Arow);A((fix(n+1)/2),:)=real(Arow);endA(K+1,:)=[0 real(poly(p(Na)))];endz=cplxpair(roots(b));K=floor(Nb/2);if Nb==0B=[0 0 poly(z)];elseif K*2==NbB=zeros(K,3);for n=1:2:NbBrow=z(n:1:n+1,:);Brow=poly(Brow);B((fix(n+1)/2),:)=real(Brow);endelseif Nb==1B=[0 real(poly(z))];elseB=zeros(K+1,3);for n=1:2:2*KBrow=z(n:1:n+1,:);Brow=poly(Brow);B((fix(n+1)/2),:)=real(Brow);endB(K+1,:)=[0 real(poly(z(Nb)))];End计算系统函数的幅度响应和相位响应子程序:function [db,mag,pha,w]=freqs_m(b,a,wmax)w1=0:500;w=w1*wmax/500;h=freqs(b,a,w);mag=abs(h);db=20*log10((mag+eps)/max(mag));pha=angle(h);脉冲响应不变法程序:function [b,a]=imp_invr(c,d,T)[R,p,k]=residue(c,d);p=exp(p*T);[b,a]=residuez(R,p,k);b=real(b).*T;数字滤波器响应子程序:function [db,mag,pha,grd,w]=freqz_m(b,a);[H,w]=freqz(b,a,1000,'whole');H=(H(1:501))';w=(w(1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);grd=grpdelay(b,a,w);直接转换成并联型子程序:function [C,B,A]=dir2par(b,a)M=length(b);N=length(a);[r1,p1,C]=residuez(b,a);p=cplxpair(p1,10000000*eps);x=cplxcomp(p1,p);r=r1(x);K=floor(N/2);B=zeros(K,2);A=zeros(K,3);if K*2==Nfor i=1:2:N-2br=r(i:1:i+1,:);ar=p(i:1:i+1,:);[br,ar]=residuez(br,ar,[]);B((fix(i+1)/2),:)real(br');A((fix(i+1)/2),:)real(ar');end[br,ar]=residuez(r(N-1),p(N-1),[]);B(K,:)=[real(br') 0];A(K,:)=[real(ar') 0];elsefor i=1:2:N-1br=r(i:1:i+1,:);ar=p(i:1:i+1,:);[br,ar]=residuez(br,ar,[]);B((fix(i+1)/2),:)real(br);A((fix(i+1)/2),:)real(ar);endEnd比较两个含同样标量元素但(可能)有不同下标的复数对及其相位留数向量子程序:function I=cplxcomp(p1,p2)I=[];for i=1:length(p2)for j=1:length(p1)if(abs(p1(j)-p2(i))<0.0001)I=[I,j];endendendI=I';双线性变换巴特沃斯低通滤波器设计:巴特沃思模拟滤波器的设计子程序:function [b,a]=afd_butt(wp,ws,Rp,rs)if wp<=0error('通带边缘必须大于0');endif ws<=wperror('阻带边缘必须大于通带边缘');endif(Rp<=0)|(Rs<0)error('通带波动或阻带衰减必须大于0');endN=ceil((log10((10^(Rp/10)-1)/(10^(Rs/10)-1)))/(2*log10(wp/ws))); fprintf('\n***Butterworth Filter Order=%2.0f\n',N);OmegaC=wp/((10^(Rp/10)-1)^(1/(2*N)));[b,a]=u_buttap(N,OmegaC)设计非归一化巴特沃思模拟低通滤波器原型子程序:function [b,a]=u_buttap(N,OmegaC)[z,p,k]=buttap(N);p=p*OmegaC;k=k*OmegaC^N;B=real(poly(z));b0=k;b=k*B;a=real(poly(p));直接型到级联型形式的转换:function [b0,B,A]=dir2cas(b,a)b0=b(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0;M=length(b);N=length(a);if N>Mb=[b,zeros(1,N-M)];a=[a,zeros(1,M-N)];elseNM=0;endk=floor(N/2);B=zeros(k,3);A=zeros(k,3);if k*2==Nb=[b,0];a=[a,0];endbroots=cplxpair(roots(b));aroots=cplxpair(roots(a));for i=1:2:2*kbr=broots(i:1:i+1,:);br=real(polt(br));B((fix(i+1)/2),:)=br;ar=aroots(i:1:i+1,:);ar=real(polt(ar));A((fix(i+1)/2),:)=ar;Endfunction [db,mag,pha,grd,w]=freqz_m(b,a)[h,w]=freqz(b,a,1000,'whole');h=(h(1:501))';w=(w(1:501))';mag=abs(h);db=20*log10((mag+eps)/max(mag));pha=angle(h);grd=grdelay(b,a,w);设计一个巴特沃思高通滤波器,要求通带截止频率为0.6pi,通带内衰减不大于1dB,阻带·起始频率为0.4pi,阻带内衰减不小于15dB,T=1:>> wp=0.6*pi;ws=0.4*pi;>> Rp=1;Rs=15;T=1;>> [N,wn]=buttord(wp/pi,ws/pi,Rp,Rs) 计算巴特沃思滤波器阶数和截止频率N =4wn =>> [b,a]=butter(N,wn,'high'); 频率变换法计算巴特沃思高通滤波器>> [C,B,A]=dir2cas(b,a)C =0.0751B =1.0000 -2.0000 1.00001.0000 -2.0000 1.0000A =1.0000 0.1562 0.44881.0000 0.1124 0.0425>> [db,mag,pha,grd,w]=freqz_m(b,a);>> subplot(2,1,1);plot(w/pi,mag);>> subplot(2,1,2);plot(w/pi,db);椭圆带通滤波器的设计--ellip函数的应用:>> ws=[0.3*pi 0.75*pi]; 数字阻带边缘频率>> wp=[0.4*pi 0.6*pi]; 数字通带边缘频率>> Rp=1;Rs=40;>> Ripple=10^(-Rp/20); 通带波动>> Attn=10^(-Rs/20); 阻带衰减>> [N,wn]=ellipord(wp/pi,ws/pi,Rp,Rs) 计算椭圆滤波器参数N =4wn =0.4000 0.6000>> [b,a]=ellip(N,Rp,Rs,wn); 数字椭圆滤波器的设计>> [b0,B,A]=dir2cas(b,a) 级联形式实现b0 =0.0197B =1.0000 1.5066 1.00001.0000 0.9268 1.00001.0000 -0.9268 1.00001.0000 -1.5066 1.0000A =1.0000 0.5963 0.93991.0000 0.2774 0.79291.0000 -0.2774 0.79291.0000 -0.5963 0.9399>> figure(1);>> [db,mag,pha,grd,w]=freqz_m(b,a);>> subplot(2,2,1);plot(w/pi,mag);>> grid on;>> subplot(2,2,3);plot(w/pi,db);grid on;>> subplot(2,2,2);plot(w/pi,pha/pi);grid on;>> subplot(2,2,4);plot(w/pi,grd);设计一个巴特沃思带阻滤波器,要求通带上下截止频率为0.8pi、0.2pi,通带内衰减不大于1dB,阻带上起始频率为0.7pi、0.4pi,阻带内衰减不小于30dB:>> wp=[0.2*pi 0.8*pi];>> ws=[0.4*pi 0.7*pi];>> Rp=1;Rs=30;>> [N,wn]=buttord(wp/pi,ws/pi,Rp,Rs);>> [b,a]=butter(N,wn,'stop');>> [C,B,A]=dir2cas(b,a)C =0.0394B =1.0000 0.3559 0.99941.0000 0.3547 1.00401.0000 0.3522 0.99541.0000 0.3499 1.00461.0000 0.3475 0.99601.0000 0.3463 1.0006A =1.0000 1.3568 0.79281.0000 1.0330 0.46331.0000 0.6180 0.17751.0000 -0.2493 0.11131.0000 -0.6617 0.37551.0000 -0.9782 0.7446>> [db,mag,pha,grd,w]=freqz_m(b,a); >> subplot(2,1,1);plot(w/pi,mag);>> subplot(2,1,2);plot(w/pi);数字低通---数字带阻:function [bz,az]=zmapping(bZ,aZ,Nz,Dz) bzord=(length(bZ)-1)*(length(Nz)-1); azord=(length(aZ)-1)*(length(Dz)-1);bz=zeros(1,bzord+1);for k=0:bzordpln=[1];for i=0:k-1pln=conv(pln,Nz);endpld=[1];for i=0:bzord-k-1pld=conv(pld,Dz);endbz=bz+bZ(k+1)*conv(pln,pld); endfor k=0:azordpln=[1];for i=0:k-1pln=conv(pln,Nz);endpld=[1];for i=0:azord-k-1pld=conv(pld,Dz);endaz=az+aZ(k+1)*conv(pln,pld); endall=az(1);az=az/az1;bz=bz/az1;线性相位FIR滤波器的幅度特性:function pzkplot(num,den)hold on;axis('square');x=-1:0.01:1;y=(1-x.^2).^0.5;y1=-(1-x.^2).^0.5;plot(x,y,'b',x,y1,'b');num1=length(num);den1=length(den);if(num1>1)z=roots(num);elsez=0;endif(den1>1)p=roots(den);elsep=0;endif(num>1&den1>1)r_max_z=max(abs(real(z)));i_max_z=max(abs(imag(z)));a_max_z=max(r_max_z,i_max_z);r_max_p=max(abs(real(p)));i_max_p=max(abs(imag(p)));a_max_p=max(r_max_p,i_max_p);a_max=max(a_max_z,a_max_p);elseif (num1>1)r_max_z=max(abs(real(z)));i_max_z=max(abs(imag(z)));a_max=max(r_max_z,i_max_z);elser_max_p=max(abs(real(p)));i_max_p=max(abs(imag(p)));a_max=max(r_max_p,i_max_p);endaxis([-a_max a_max -a_max a_max]);plot([-a_max a_max],[0 0],'b');plot([0 0],[-a_max a_max],'b');plot([-a_max a_max],[a_max a_max],'b');plot([a_max a_max],[-a_max a_max],'b');Lz=length(z);for i=1:Lz;plot(real(z(i)),imag(z(i)),'bo');endLp=length(p);for j=1:Lpplot(real(p(j)),imag(p(j)),'bx');endtitle('The zeros-pole plot');xlabel('虚部');ylabel('实部');function [Hr,w,a,L]=Hr_Type1(h)M=length(h);L=(M-1)/2;a=[h(L+1) 2*h(L:-1:1)];n=[0:1:L];w=[0:1:500]'*pi/500;Hr=cos(w*n)*a';设计I型线性相位FIR滤波器:>> h=[-4 1 -1 -2 5 6 5 -2 -1 1 -4];>> M=length(h);n=0:M-1;>> [Hr,w,a,L]=Hr_Type1(h);>> amax=max(a)+1;>> amin=min(a)-1;>> subplot(2,2,1);stem(n,h);>> axis([-1 2*L+1 amin amax]);text(2*L+1.5,amin,'n'); >> xlabel('n');ylabel('h(n)');title('脉冲响应');>> subplot(2,2,3);stem(0:L,a);>> axis([-1 2*L+1 amin amax]);>> xlabel('n');ylabel('a(n)');title('a(n) 系数');>> subplot(2,2,2);plot(w/pi,Hr);>> grid on;text(1.05,-20,'频率pi');>> xlabel('频率');ylabel('Hr');title('I 型振幅响应');>> subplot(2,2,4);pzkplot(h,1);>> title('零极点分布');function [hr,w,b,L]=Hr_Type2(h)M=length(h);L=M/2;b=2*h(L:-1:1);n=[1:1:L];n=n-0.5;w=[0:1:500]'*pi/500;hr=cos(w*n)*b';II型线性相位FIR滤波器:>> h=[-4 1 -1 -2 5 6 5 -2 -1 1 -4];>> M=length(h);n=0:M-1;>> [Hr,w,b,L]=Hr_Type2(h);Warning: Integer operands are required for colon operator when used as index. > In Hr_Type2 at 2>> bmax=max(b)+1;bmin=min(b)-1;>> subplot(2,2,1);stem(n,h);axis([-1 2*L+1 bmin bmax]);text(2*L+1.5,bmin,'n');xlabel('n');ylabel('h(n)');title('脉冲响应');>> subplot(2,2,3);stem(1:L,b);axis([-1 2*L+1 bmin bmax]);xlabel('n');ylabel('b(n)');title('b(n) 系数');>> subplot(2,2,2);plot(w/pi,Hr);grid on;text(1.05,-20,'频率pi');xlabel('频率');ylabel('Hr');title('II 型振幅响应');>> subplot(2,2,4);pzkplot(h,1);title('零极点分布');function [hr,w,c,L]=Hr_Type3(h)M=length(h);L=(M-1)/2;b=2*h(L+1:-1:1);n=[1:1:L];w=[0:1:500]'*pi/500;hr=cos(w*n)*c';用MA TLAB编程绘制各种窗函数的形状。
利用MATLAB实现数字低通滤波器的设计
西南石油大学实验报告一实验目的:1学习用Matlab直接设计模拟滤波器和数字滤波器。
2学习用冲激响应不变法和双线性变换法的Matlab的实现。
二实验内容:设计满足下列指标的数字低通滤波器:Wp=0.2*pi, Rp=1db Ws=0.5*pi Rs=20db Fs=1khz1.利用B、C1型设计出模拟低通滤波器,采用冲激响应不变法、双线性发转换成数字低通滤波器。
2.直接设计出B、C1型数字低通滤波器。
三实验步骤:程序1Wp=2*pi*0.1*1000;Ws=2*pi*0.25*1000;Rp=1;Rs=20;[N,Wn]=buttord(Wp,Ws,Rp,Rs,'s');[z,p,k]=buttap(N);[B,A]=butter(N,Wn,'s');freq1=linspace(0,Wp,5);freq2=linspace(Wp,Ws,15);freq3=linspace(Ws,10*pi*2,25);h1=20*log10(abs(freqs(B,A,freq1)));h2=20*log10(abs(freqs(B,A,freq2)));h3=20*log10(abs(freqs(B,A,freq3)));plot([freq1 freq2 freq3]/(2*pi),[h1,h2,h3]);grid;Xlabel('Frequency in Hz');Ylabel('gain in DB');图一程序2wp=0.2*pi;ws=0.5*pi;rp=1;rs=20;fs=1000;omegap=wp*fs;omegas=ws*fs;[N,Wn]=buttord(omegap,omegas,rp,rs,'s');[B A]=butter(N,Wn,'s');[b,a]=impinvar(B,A,fs);[h,w]=freqz(b,a,256);h=20*log10(abs(h));plot(w/pi,h);图二程序3wp=0.2*pi;ws=0.5*pi;rp=1;rs=20;fs=1000;omegap=2*fs*tan(wp/2);omegas=2*fs*tan(ws/2);[N,Wn]=cheb1ord(omegap,omegas,rp,rs,'s');[B A]=cheby1(N,rp,Wn,'s');[b,a]=bilinear(B,A,fs);[h,w]=freqz(b,a,256);h=20*log10(abs(h));plot(w/pi,h);图三程序4wp=0.2*pi;ws=0.5*pi;rp=1;rs=20;[N,Wn]=buttord(wp/pi,ws/pi,rp,rs);[B A]=butter(N,Wn);[h,w]=freqz(B,A,256);h=20*log10(abs(h));plot(w/pi,h);图四程序5Wp=0.2*pi;Ws=0.5*pi;Rp=1;Rs=20;T=0.001;Fs=1000;omegap=(2/T)*tan(Wp/2);omegas=(2/T)*tan(Ws/2);[N,Wn]=cheb1ord(omegap,omegas,Rp,Rs,'s'); [B,A]=cheby1(N,Rp,Wn,'s');[b,a]=bilinear(B,A,Fs);[h,w]=freqz(b,a,256);h1=20*log10(abs(h));plot(w/pi,h1);grid;xlabel('Digital Frequency in pi units'); ylabel('Gain in DB');axis([0 1 -50 10]);图五Wp=0.2;Ws=0.5;Rp=1;Rs=20;disp('ÇбÈÑ©·òIÐÍ')[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs)[B,A]=cheby1(N,Rp,Wn);disp('ÇбÈÑ©·òÐÍ·Ö×Ó¶àÏîʽ');fprintf('%.4e\n',B);disp('ÇбÈÑ©·ò·Öĸ¶àÏîʽ');fprintf('%.4e\n',A);w=linspace(0,0.8*pi,50);h1=20*log10(abs(freqz(B,A,w)));plot(w/pi,h1);grid;xlabel('Normalized frequency');ylabel('Gain in DB ');axis([0 0.8 -50 1]);图六四、实验小结通过本次实验,对MA TLAB软件有了进一步的了解,也在不断的实践中,更多的熟悉了MATLAB的编程,在编程方面一点点的有了进步。
matlab iir低通滤波器设计
I. 简介Matlab是一种非常常用的科学计算软件,它广泛用于信号处理、图像处理、控制系统等领域。
在信号处理中,IIR(Infinite Impulse Response)滤波器是一种常见的数字滤波器,常被用于模拟滤波、数字滤波等应用中。
这篇文章将介绍如何使用Matlab进行IIR低通滤波器的设计。
II. 什么是IIR低通滤波器1. IIR滤波器IIR滤波器是一种数字滤波器,其特点是其单位脉冲响应是无限长的。
它通常具有较为复杂的频率响应特性,且具有较小的阶数,能够更好地逼近某些复杂的频率响应曲线。
IIR滤波器分为低通滤波器、高通滤波器、带通滤波器和带阻滤波器等。
2. 低通滤波器低通滤波器是一种常见的滤波器,其特点是只允许低频信号通过,而抑制高频信号。
在信号处理中,低通滤波器常被用于去除高频噪声、提取低频信号等应用中。
III. Matlab中的IIR低通滤波器设计1. 使用Matlab进行IIR低通滤波器设计Matlab提供了丰富的信号处理工具箱,包括了数字滤波器设计工具。
在Matlab中,可以使用函数butter、cheby1、cheby2、ellip等来设计IIR低通滤波器。
2. 设计步骤设计IIR低通滤波器的一般步骤如下:a. 确定通带和阻带的频率范围b. 选择滤波器的通带和阻带的最大允许衰减c. 选择滤波器的类型(Butterworth、Chebyshev等)以及阶数d. 使用Matlab中相应的函数设计滤波器e. 对设计的滤波器进行频率响应分析IV. 实例分析以下是一个在Matlab中设计IIR低通滤波器的简单实例:设计IIR低通滤波器fs = 1000; 采样频率fpass = 100; 通带截止频率fstop = 200; 阻带截止频率apass = 1; 通带最大允许衰减astop = 80; 阻带最小要求衰减[num, den] = butter(4, fpass/(fs/2), 'low');freqz(num, den, 512, fs); 绘制滤波器频率响应曲线V. 结论使用Matlab进行IIR低通滤波器设计是一种简单而有效的方法。
MATLAB数字滤波器的设计
MATLAB数字滤波器的设计一、数字滤波器设计数字滤波器就是由数字乘法器、加法器和延时单位组成的一种装置。
它的作用是对输入离散信号的数字代码进行运算处理,从而达到改变信号频谱的目的。
数字滤波器可以用计算机软件实现或者大规模集成数字硬件实时实现。
在数字滤波器处理模拟信号的时候,先要对输入的模拟信号进行限带、抽样和模数转换。
数字滤波器输入信号的抽样率要大于信号带宽的两倍,它的频率响应还要具有以抽样频率为间隔的周期重复特性,并且它还要以1/2抽样频率点呈镜像对称。
数字滤波器具有高精度、高可靠性、可控制改变特性、便于集成等优点。
现在数字滤波器的使用范围比较广,例如,在语言信号处理、图像信号处理、医学生物信号处理以及其他应用领域。
数字滤波器有低通、高通、带通、带阻和全通等类型。
函数butter可用于设计Butterworth滤波器,butter函数可以设计低通、高通、带通和带阻的数字和模拟滤波器,用butter函数设计的滤波器通带内的幅度响应最大限度的平滑,但损失了截止频率处的下降斜度。
函数chebyl可用于设计Chebyshev I型滤波器,chebyl函数可以设计低通、高通、带通和带阻的数字和模拟Chebyshev I型滤波器,用chebyl函数设计的Chebyshev I型滤波器的通带内为等波纹,它的阻带内为单调。
Chebyshev I型滤波器的下降斜度比Chebyshev II型的要大,但是Chebyshev II的代价是通带内波纹比较大。
除了典型设计以外,MATLAB信号处理工具箱提供了几个直接设计IIR数字滤波器的函数,直接调用就可以设计滤波器,这为设计通用滤波器提供了方便。
在MATLAB中设计IIR滤波器时要设计出巴特沃斯滤波器可以使用Butterworth函数,要设计契比雪夫I型滤波器时可以使用Cheby1函数,设计契比雪夫II型滤波器要使用Cheby2,而椭圆滤波器则需要ellipord函数。
设计数字低通滤波器(用matlab实现)
DSP 设计滤波器报告姓名:张胜男 班级:07级电信(1)班 学号:078319120一·低通滤波器的设计(一)实验目的:掌握IIR 数字低通滤波器的设计方法。
(二)实验原理:1、滤波器的分类滤波器分两大类:经典滤波器和现代滤波器。
经典滤波器是假定输入信号)(n x 中的有用成分和希望取出的成分各自占有不同的频带。
这样,当)(n x 通过一个线性系统(即滤波器)后可讲欲去除的成分有效的去除。
现代滤波器理论研究的主要内容是从含有噪声的数据记录(又称时间序列)中估计出信号的某些特征或信号本身。
经典滤波器分为低通、高通、带通、带阻滤波器。
每一种又有模拟滤波器(AF )和数字滤波器(DF )。
对数字滤波器,又有IIR 滤波器和FIR 滤波器。
IIR DF 的转移函数是:∑∑=-=-+==Nk k k M r r rz a z b z X z Y z H 101)()()(FIR DF 的转移函数是:∑-=-=10)()(N n nz n h z HFIR 滤波器可以对给定的频率特性直接进行设计,而IIR 滤波器目前最通用的方法是利用已经很成熟的模拟滤波器的设计方法进行设计。
2、滤波器的技术要求低通滤波器:p ω:通带截止频率(又称通带上限频率) s ω:阻带下限截止频率 p α:通带允许的最大衰减s α:阻带允许的最小衰减 (p α,s α的单位dB ) p Ω:通带上限角频率 s Ω:阻带下限角频率 (s p p T ω=Ω,s s s T ω=Ω)即 C p p F ωπ2=Ω C s s F ωπ2=Ω3、IIR 数字滤波器的设计步骤:1)按一定规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标。
2)根据转换后的技术指标设计模拟低通滤波器)(s G ;3)再按一定的规则将)(s G 转换成)(z H 。
4)若是高通、带通或带阻数字滤波器则将它们的技术指标先转化为低通模拟滤波器的技术指标,然后按上述步骤2)设计出低通)(s G ,再将)(s G 转换为所需的)(z H 。
FIR数字滤波器的MATLAB详细程序
%FIR-DF hd=ideal_lp(wc,M);Fs=32000;a=2*pi*6500;b=2*pi*7000;c=2*pi*9000;n1=0:4096;t=n1/Fs;xa=cos(a*t)+cos(b*t)+cos(c*t);%+cos(d*t);figure(1);y=fft(xa);plot(abs(y));title('fft(xa)变换');grid;Rp=0.25;As=52;%ws=pi*6600*2/Fs;wp=pi*6700*2/Fs;%gao tong%ws=pi*6700*2/Fs;wp=pi*6600*2/Fs;%di tongws=[pi*6600*2/Fs,pi*7400*2/Fs];wp=[pi*6800*2/Fs,pi*7200*2/Fs];%dai tong%ws=[pi*6600*2/Fs,pi*7400*2/Fs,pi*8600*2/Fs,pi*9400*2/Fs];wp=[pi*6800*2/Fs,pi*7200*2/Fs, pi*8800*2/Fs,pi*9200*2/Fs];%dai tong%ws=[pi*6800*2/Fs,pi*7200*2/Fs];wp=[pi*6600*2/Fs,pi*7400*2/Fs];%dai zuDB=abs(ws(1)-wp(1));%DB=abs(ws-wp);M0=ceil((As-8)/(2.285*DB));%阶数M=M0+mod(M0,2);%确保阶数为偶数N=M+1;%长度为奇数%0.112*(As-8.7);As>50;%0.5842*(As-21)^0.4+0.07886*(As-21);21<=As<=50;%0;As<21;if As>50beta=0.112*(As-8.7);elseif21<=As&&As<=50beta=0.5842*(As-21)^0.4+0.07886*(As-21);elsebeta=0;endwn=kaiser(N,beta);wc=(wp+ws)/2/pi;%fir1使用wc归一化%hn=fir1(M,wc,'high',wn);%gao tonghn=fir1(M,wc,wn);%dai tong,di tong%hn=fir1(M,wc,'stop',wn);%dai zufigure(2);plot(abs(wn));title('窗函数');grid;figure(3);[hh,w]=freqz(hn,1);freqzplot(hh,w,'linear');title('幅频/相频特性');%Hn=fft(hn);%plot(abs(Hn));title('FFT(Hn)');%freqz(hn,1);title('幅频特性');y=filter(hn,1,xa);figure(4);y1=fft(y);plot(abs(y1));title('高通滤波器滤波之后');grid;%=========验证=====================%验证%w_kai=wn';%hd=hn;%h=hd.*w_kai;%freqz_mN1=1000;[H,W]=freqz(hn,1,N1,'whole');%[H,W]=freqz(B,A,N1,'whole');delta_w=2*pi/N1;db=20*log10(abs(H)/max(abs(H))+eps);%通带波纹db数%rp1=-(min(db(wp/delta_w+1:1:501)))%gao tong yan zheng%as1=-round(max(db(1:1:ws/delta_w+1)))%rp1=-(min(db(1:1:wp/delta_w+1)))%di tong yan zheng%as1=-round(max(db(ws/delta_w+1:1:501)))rp1=-(min(db(wp(1)/delta_w+1:1:wp(2)/delta_w+1)))%dai tong yan zhengas1=round(max(-max(db(1:1:ws(1)/delta_w+1)),-max(db(ws(2)/delta_w+1:1:501))))%rp1=-min(min(db(1:1:wp(1)/delta_w+1)),min(db(wp(2)/delta_w+1:1: 501)))%dai zu yan zheng%as1=-round(max(db(ws(1)/delta_w+1:1:ws(2)/delta_w+1)))%N1=4096;%Hw=fft(h,N1);%wk=2*pi*[0:N1-1]/N1;%hgw=Hw.*exp(j*wk*N/2);%rp=max(20*log10(abs(hgw)))%%hgmin=min(real(hgw));%rs=20*log10(abs(hgmin))%w_kai=wn';%hd=hn;%h=hd.*w_kai;%d_w=2*pi/1000;%[db,w0]=freqz(h,1);%rp1=-(min(db(wp/d_w+1:1:501))) %as1=-round(max(db(1:1:ws/d_w+1)))。
基于MATLAB的数字滤波器的设计
基于MATLAB 的数字滤波器的设计1 引言数字滤波器是指完成信号滤波处理功能的,用有限精度算法实现的离散线性非时变系统,其输入是一组(由模拟信号取样和量化的)数字量,其输出是经过变换或说处理的另一组数字量。
数字滤波器具有稳定性高、精度高、灵活性大等突出优点。
这里所说的数字滤波器是指理想带通,低通等的频率选择数字滤波器。
数字滤波器设计的一个重要步骤是确定一个可实现的传输函数H(z),这个确定传输函数H(z)的过程称为数字滤波器设计。
数字滤波器的一般设计过程为:(1)按照实际需要,确定滤波器的性能要求(通常在频域内给定数字滤波的性能要求)。
(2)寻找一满足预定性能要求的离散时间线性系统。
(3)用有限精度的运算实现所设计的系统。
(4)通过模拟,验证所设计的系统是否符合给定性能要求。
2 数字滤波器的设计滤波器分为两种,分别为模拟滤波器和数字滤波器。
数字滤波器可以理解为是一个计算程序或算法,将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列,并在转化的过程中,使信号按预定的形式变化。
数字滤波器有多种分类,从数字滤波器功能上分可分为低通、高通、带阻、带通滤波器,根据数字滤波器冲激响应的时域特征,可将数字滤波器分为两种,即无限长冲激响应滤波器(IIR )和有限长冲激响应滤波器(FIR )。
数字滤波器指标:一般来说,滤波器的幅频特性是分段常数的,以低通为例,在通带内逼近于1,阻带内逼近与0,实际设计的滤波器并非是锐截止的通带和阻带两个范围,两者之间总有一个过渡带。
在设计滤波器时事先给定幅频特性允许误差,在通带范围内幅度响应以误差逼近于1,在阻带内幅1σ度响应以误差逼近于0。
2σ (1)πσσ≤≤≤≤≤≤-w w e H w w e H r jwc jw ,2|)(|,1|)(|11式中wc 和wr 分别为通带边界频率和阻带边界频率,wr-wc 为过渡带。
在具体的技术指标中往往用通带波动来表示,用最小阻带衰减At 来表示,其具体的对应公式这里就不详述了。
matlab 数据滤波处理
matlab数据滤波处理在MATLAB中,数据滤波可以使用不同的方法和函数来实现。
下面是几种常见的数据滤波处理方法:1.移动平均滤波(Moving Average Filter):移动平均滤波是一种简单的滤波方法,通过计算数据序列中相邻数据点的平均值来平滑数据。
可以使用smoothdata函数实现移动平均滤波。
示例:创建示例数据data=randn(1,100);100个随机数应用移动平均滤波smoothed_data=smoothdata(data,'movmean',5);使用移动窗口大小为5的移动平均2.中值滤波(Median Filter):中值滤波是一种非线性滤波方法,将每个数据点替换为相邻数据点的中值。
可以使用medfilt1函数进行中值滤波处理。
示例:创建示例数据data=randn(1,100);100个随机数应用中值滤波smoothed_data=medfilt1(data,5);使用窗口大小为5的中值滤波3.低通滤波器(Low-pass Filter):低通滤波器可以滤除高频噪声,保留信号的低频成分。
MATLAB中可以使用filter函数设计和应用数字低通滤波器。
示例:创建示例数据data=randn(1,100);100个随机数设计低通滤波器fc=0.1;截止频率fs=1;采样频率[b,a]=butter(5,fc/(fs/2));设计5阶巴特沃斯低通滤波器应用低通滤波器smoothed_data=filter(b,a,data);以上示例中的滤波方法和参数可以根据数据的特性和需求进行调整。
使用不同的滤波方法可能需要更多的参数调整和信号处理知识。
根据具体情况,可以选择合适的滤波方法来平滑或处理数据。
matlab滤波器设计命令
matlab滤波器设计命令Matlab滤波器设计命令滤波器是数字信号处理中常用的工具,用于去除信号中的噪声、频率干扰或其他不需要的成分。
Matlab提供了一系列有用的滤波器设计命令,使用户能够轻松设计并应用各种类型的滤波器。
在本文中,我们将详细介绍Matlab中常用的滤波器设计命令,包括滤波器设计函数、滤波器类型和设计过程。
I. Matlab中常用的滤波器设计函数在Matlab中,有几种函数可用于设计滤波器,其中最常用的函数是`designfilt`函数和`fir1`函数。
1. designfilt函数`designfilt`函数是Matlab中最灵活和功能强大的滤波器设计函数之一,可用于设计各种类型的IIR和FIR滤波器。
它的基本语法如下:`filt = designfilt(FilterType, 'PropertyName', PropertyValue, ...)`其中,`FilterType`代表滤波器类型,包括低通滤波器(Lowpass)、高通滤波器(Highpass)、带通滤波器(Bandpass)、带阻滤波器(Bandstop)等。
`PropertyName`和`PropertyValue`是可选的参数,用于设置滤波器的各种属性,如阶数(Order)、截止频率(CutoffFrequency)、通带和阻带的最大衰减(MaximumAttenuation)等。
下面是一个使用`designfilt`函数设计低通滤波器的例子:Fs = 1000; 采样频率Fpass = 20; 通带截止频率Fstop = 30; 阻带截止频率designfilt('lowpassiir', 'FilterOrder', 4, 'PassbandFrequency', Fpass, 'StopbandFrequency', Fstop, 'SampleRate', Fs)该命令将设计一个4阶的低通IIR滤波器,其通带截止频率为20Hz,阻带截止频率为30Hz,采样频率为1000Hz。
基于matlab的数字滤波器设计与仿真
滤波器仿真实验报告题目:基于Matlab的数字滤波器设计与仿真姓名:学号:指导教师:基于Matlab 的数字滤波器设计与仿真数字信号处理在通信、语音、图像,自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。
在数字信号处理中,滤波器占有非常重要的地位。
数字滤波器在语音和图像处理、HDTV 、模式识别、频谱分析等方面得到广泛应用。
所谓滤波器,是指对输入信号进行滤波的软件或硬件。
所谓数字滤波器,是指输入、输出均为数字信号,通过一定运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分的硬件或软件。
由于数字滤波信号形式与实现滤波方法与模拟滤波方法不同,数字滤波器具有比模拟滤波器精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配以及可以避免模拟滤波器所无法克服的电压漂移和噪声问题。
数字滤波器的设计,其实质是数学逼近理论的应用,通过计算使物理可实现的实际滤波器频率特性逼近理想的或给定的频率特性,以达到去除干扰提取有用信号的目的。
1 实验与仿真原理1.1 数字滤波器的工作原理在数字滤波中,我们主要讨论离散时间序列。
如图1所示。
设输入序列为()n x ,离散或数字滤波器对单位抽样序列()n δ的响应为()n h 。
因()n δ在时域离散信号和系统中所起的作用相当于单位冲激函数在时域连续信号和系统中所起的作用。
图1 数字滤波器原理数字滤波器的序列()n y将是这两个序列的离散卷积,即()()()∑∞∞=-=kk n x k h n y同样,两个序列卷积的z 变换等于个自z 变换的乘积,即()()()z X z H z Y =用T j e z ω=代入上式,其中T 为抽样周期,则得到()()()T j T j T j e X e H e Y ωωω=式中()Tj e Xω和()T j e Y ω分别为数字滤波器输入序列和输出序列的频谱,而()Tj e H ω为单位抽样序列响应()n h 的频谱。
MATLAB 用频率采样法设计一带通数字滤波器(程序)
以MATLAB 用频率采样法设计一带通数字滤波器。
低通带边缘:w1p=0.2*pi,低阻带边缘:w1s=0.35*pi,高通带边缘w2p=0.65*pi,高阻带边缘w2s=0.8*pi;设计过渡带中的频率样本值为T1=0.109021,T2=0.59417456cheb1% wp1=0.45*pi;wp2=0.65*pi;ws1=0.3*pi;ws2=0.75*pi;Rp=1;Rs=40% =============双线型变换法========================================= wp1=0.45*pi; wp2=0.65*pi;ws1=0.3*pi; ws2=0.75*pi;Rp=1; Rs=40;Wp1=tan(wp1/2); Wp2=tan(wp2/2);Ws1=tan(ws1/2); Ws2=tan(ws2/2);BW=Wp2-Wp1; W0=Wp1*Wp2; W00=sqrt(W0);WP=1; WS=WP*(W0^2-Ws1^2)/(Ws1*BW);[N,Wn]=cheb1ord(WP,WS,Rp,Rs,'s');[B,A]=cheby1(N,Rp,Wn,'s');[BT,AT]=lp2bp(B,A,W00,BW);[num,den]=bilinear(BT,AT,0.5);[h,omega]=freqz(num,den,64);subplot(2,2,1);stem(omega/pi,abs(h));xlabel('\omega/\pi');ylabel('|H(z)|');subplot(2,2,2);stem(omega/pi,20*log10(abs(h)));xlabel('\omega/\pi');ylabel('增益.dB');% =============直接法=================================wp1=0.45*pi; wp2=0.65*pi;ws1=0.3*pi; ws2=0.75*pi;Rp=1; Rs=40;Wp=[wp1/pi,wp2/pi]; Ws=[ws1/pi,ws2/pi];[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs);[B,A]=cheby1(N,Rp,Wn);[h,omega]=freqz(B,A,64);subplot(2,2,3);stem(omega/pi,abs(h));xlabel('\omega/\pi');ylabel('|H(z)|');subplot(2,2,4);stem(omega/pi,20*log10(abs(h)));xlabel('\omega/\pi');ylabel('增益.dB');%cheby2%% wp1=0.45*pi;wp2=0.65*pi;ws1=0.3*pi;ws2=0.75*pi;Rp=1;Rs=40% =============双线型变换法=========================================wp1=0.45*pi; wp2=0.65*pi;ws1=0.3*pi; ws2=0.75*pi;Rp=1; Rs=40;Wp1=tan(wp1/2); Wp2=tan(wp2/2);Ws1=tan(ws1/2); Ws2=tan(ws2/2);BW=Wp2-Wp1; W0=Wp1*Wp2; W00=sqrt(W0);WP=1; WS=WP*(W0^2-Ws1^2)/(Ws1*BW);[N,Wn]=cheb2ord(WP,WS,Rp,Rs,'s');[B,A]=cheby2(N,Rs,Wn,'s');[BT,AT]=lp2bp(B,A,W00,BW);[num,den]=bilinear(BT,AT,0.5);[h,omega]=freqz(num,den,64);subplot(2,2,1);stem(omega/pi,abs(h));xlabel('\omega/\pi');ylabel('|H(z)|');subplot(2,2,2);stem(omega/pi,20*log10(abs(h)));axis([0 1 -100 0]);xlabel('\omega/\pi');ylabel('增益.dB');% =============直接法=================================wp1=0.45*pi; wp2=0.65*pi;ws1=0.3*pi; ws2=0.75*pi;Rp=1; Rs=40;Wp=[wp1/pi,wp2/pi]; Ws=[ws1/pi,ws2/pi];[N,Wn]=cheb2ord(Wp,Ws,Rp,Rs);[B,A]=cheby2(N,Rs,Wn);[h,omega]=freqz(B,A,64);subplot(2,2,3);stem(omega/pi,abs(h));xlabel('\omega/\pi');ylabel('|H(z)|');subplot(2,2,4);stem(omega/pi,20*log10(abs(h)));axis([0 1 -100 0]);xlabel('\omega/\pi');ylabel('增益.dB');%椭圆%% wp1=0.45*pi;wp2=0.65*pi;ws1=0.3*pi;ws2=0.75*pi;Rp=1;Rs=40% =============双线型变换法========================================= wp1=0.45*pi; wp2=0.65*pi;ws1=0.3*pi; ws2=0.75*pi;Rp=1; Rs=40;Wp1=tan(wp1/2); Wp2=tan(wp2/2);Ws1=tan(ws1/2); Ws2=tan(ws2/2);BW=Wp2-Wp1; W0=Wp1*Wp2; W00=sqrt(W0);WP=1; WS=WP*(W0^2-Ws1^2)/(Ws1*BW);[N,Wn]=ellipord(WP,WS,Rp,Rs,'s');[B,A]=ellip(N,Rp,Rs,Wn,'s');[BT,AT]=lp2bp(B,A,W00,BW);[num,den]=bilinear(BT,AT,0.5);[h,omega]=freqz(num,den,64);subplot(2,2,1);stem(omega/pi,abs(h));grid;xlabel('\omega/\pi');ylabel('|H(z)|');subplot(2,2,2);stem(omega/pi,20*log10(abs(h)));grid;xlabel('\omega/\pi');ylabel('增益.dB');% =============直接法================================= wp1=0.45*pi; wp2=0.65*pi;ws1=0.3*pi; ws2=0.75*pi;Rp=1; Rs=40;Wp=[wp1/pi,wp2/pi]; Ws=[ws1/pi,ws2/pi];[N,Wn]=ellipord(Wp,Ws,Rp,Rs);[B,A]=ellip(N,Rp,Rs,Wn);[h,omega]=freqz(B,A,64);subplot(2,2,3);stem(omega/pi,abs(h));grid;xlabel('\omega/\pi');ylabel('|H(z)|');subplot(2,2,4);stem(omega/pi,20*log10(abs(h)));grid;xlabel('\omega/\pi');ylabel('增益.dB');一个例程,只需要改一改参数就行了clear all;f=[0 0.19 0.2 0.3 0.31 0.59 0.6 0.8 0.81 1];% 给定频率轴分点;m=[0 0 1 1 0 0 1 1 0 0];% 给定在这些频率分点上理想的幅频响应N1=30;N2=90;% 取两种不同的滤波器长度;b1=fir2(N1,f,m);b2=fir2(N2,f,m);% 得到两个滤波器;subplot(311);stem(b1,'.');grid;subplot(312);stem(b2,'.');grid;M=128;[h1,w]=freqz(b1,1,M,1);[h2,w]=freqz(b2,1,M,1);subplot(313);plot(w,abs(h1),'b-',w,abs(h2),'g-');grid;其中,f是归依化以后的频率通过数字滤波器的采样频率算出来,根据通带和阻带算好f和m就行了看一看help,这个函数应该有窗函数的选择默认情况下是汉明窗设计一个模拟的butterworth低通滤波器,设阶数是N,OmegaC=50*pi*2;[z,p,k]=buttap(N);p=p*OmegaC;k=k*OmegaC^N;B=real(poly(z));b0=k;b=k*B;a=real(poly(p));这样一个最简单的低通滤波器就设计好了,b是传递函数的分子系数,a是分母,都是按s的降幂排列的。
基于MATLAB的数字插值滤波器设计
EDA 课程设计报告滤波器设计参数:根据要求,要设计一个输入8位,输出8位的17阶线性相位FIR 滤波器,所以采用图2(a)的方式,其中输入信号范围为:[±99,0,0,0, ±70,0,0,0, ±99,0,0,0, ±70,…],此滤波器 Fs 为44kHz,Fc 为10.4kHz 。
(一)FIR 数字滤波器理论简述有限冲激响应(FIR )数字滤波器和无限冲激响应(IIR )数字滤波器广泛应用于数字信号处理系统中。
IIR 数字滤波器方便简单,但它相位的非线性,要求采用全通网络进行相位校正,且稳定性难以保障。
FIR 滤波器具有很好的线性相位特性,使得它越来越受到广泛的重视。
有限冲击响应(FIR )滤波器的特点:1 既具有严格的线性相位,又具有任意的幅度;2 FIR 滤波器的单位抽样响应是有限长的,因而滤波器性能稳定;3只要经过一定的延时,任何非因果有限长序列都能变成因果的有限长序列,因而能用因果系统来实现;4 FIR 滤波器由于单位冲击响应是有限长的,因而可用快速傅里叶变换(FFT)算法来实现过滤信号,可大大提高运算效率。
5 FIR 也有利于对数字信号的处理,便于编程,用于计算的时延也小,这对实时的信号处理很重要。
6 FIR 滤波器比较大的缺点就是阶次相对于IIR 滤波器来说要大很多。
FIR 数字滤波器是一个线性时不变系统(LTI ),N 阶因果有限冲激响应滤波器可以用传输函数H (z )来描述,()()Nk k H z h k z -==∑(0.1)在时域中,上述有限冲激响应滤波器的输入输出关系如下:[][][][][]Nk y n x n h n x k h n k ==*=-∑(0.2)其中,x [n ]和y [n ]分别是输入和输出序列。
N 阶有限冲激响应滤波器要用N +1个系数描述,通常要用N+1个乘法器和N 个两输入加法器来实现。
乘法器的系数正好是传递函数的系数,因此这种结构称为直接型结构,可通过式(1.2)来实现,如图1。
用MATLAB设计数字FIR低通滤波器(哈明窗)
用MATLAB设计数字FIR低通滤波器(哈明窗)设计数字FIR低通滤波器,技术指标如下:Wp=0.2pi.rp=0.25db;Ws=0.3pi,as=50db.求解,如下图:程序:wp=0.2*pi;ws=0.3*pi;trwidth=ws-wp;M=ceil(6.6*pi/trwidth)+1;%朝正无穷方向取整n=[0:1:M-1];wc=(ws+wp)/2%求的截止频率(弧度)hd=ideal_lp(wc,M) %调用计算理想低通滤波器的单位取样响应wham=(hamming(M))';%哈明窗函数h=hd.*wham;[H,w]=freqz(h,1,1000,'whole'); %求Z变换频率响应函数H=(H(1:1:501))';w=(w(1:1:501))'mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);grd=grpdelay(h,1,w);deltaw=2*pi/1000;Rp=-(min(db(1:1:wp/deltaw+1)))As=-round(max(db(ws/deltaw+1:1:501)))subplot(1,1,1)subplot(2,2,1),stem(n,hd);title('理想脉冲响应')axis([0 M-1 -0.1 0.3]);xlabel('N');ylabel('hd(n)');subplot(2,2,2),stem(n,wham);title('Hamming窗')axis([0 M-1 -0 1.1]);xlabel('N');ylabel('w(n)');subplot(2,2,3),stem(n,h);title('实际脉冲响应')axis([0 M-1 -0.1 0.3]);xlabel('N');ylabel('h(n)');subplot(2,2,4),plot(w/pi,db);title('相对标尺的幅度响应');grid on;axis([0 1 -100 10]);xlabel('frequency in pi units');ylabel('decibels');set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1])%设置或修改LTI对象的属性值set(gca,'XTickMode','manual','XTick',[-50,0])set(gca,'XTickLabelMode','manual','XTickLabels',['50';'0'])子程序:function hd=ideal_lp(wc,N);a=(N-1)/2;n=[0:1:(N-1)];m=n-a+eps;hd=sin(wc*m)./(pi*m);。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
%要求设计一butterworth低通数字滤波器,wp=30hz,ws=40hz,rp=0.5,rs=40,fs=100hz。
>> wp=30;ws=40;rp=0.5;rs=40;fs=100;>> wp=30*2*pi;ws=40*2*pi;>> [n,wn]=buttord(wp,ws,rp,rs,'s');>> [z,p,k]=buttap(n);>> [num,den]=zp2tf(z,p,k);>> [num1,den1]=impinvar(num,den);Warning: The output is not correct/robust.Coeffs of B(s)/A(s) are real, but B(z)/A(z) has complex coeffs.Probable cause is rooting of high-order repeated poles in A(s).> In impinvar at 124>> [num2,den2]=bilinear(num,den,100);>> [h,w]=freqz(num1,den1);>> [h1,w1]=freqz(num2,den2);>> subplot(1,2,1);>> plot(w*fs/(2*pi),abs(h));>> subplot(1,2,2);>> plot(w1*fs/(2*pi),abs(h1));>> figure(1);>> subplot(1,2,1);>> zplane(num1,den1);>> subplot(1,2,2);>> zplane(num2,den2);%要求设计一chebyshev1低通数字滤波器,wp=30hz,ws=40hz,rp=0.5,rs=40,fs=100hz. >> wp=30;ws=40;rp=0.5;rs=40;fs=100;>> wp=30*2*pi;ws=40*2*pi;>> [n,wn]=cheb1ord(wp,ws,rp,rs,'s');>> [z,p,k]=cheb1ap(n,rp);>> [num,den]=zp2tf(z,p,k);[num1,den1]=impinvar(num,den);>> [num2,den2]=bilinear(num,den,100);[h,w]=freqz(num1,den1);[h1,w1]=freqz(num2,den2);subplot(1,2,1);plot(w*fs/(2*pi),abs(h));subplot(1,2,2);plot(w1*fs/(2*pi),abs(h1));figure(2);subplot(1,2,1);zplane(num1,den1);subplot(1,2,2);zplane(num2,den2);%要求设计一chebyshev2低通数字滤波器,wp=30hz,ws=40hz,rp=0.5,rs=40,fs=100hz. wp=30;ws=40;rp=0.5;rs=40;fs=100;wp=30*2*pi;ws=40*2*pi;[n,wn]=cheb2ord(wp,ws,rp,rs,'s');[z,p,k]=cheb2ap(n,rs);[num,den]=zp2tf(z,p,k);[num1,den1]=impinvar(num,den); [num2,den2]=bilinear(num,den,100); [h,w]=freqz(num1,den1);[h1,w1]=freqz(num2,den2);subplot(1,2,1);plot(w*fs/(2*pi),abs(h));subplot(1,2,2);plot(w1*fs/(2*pi),abs(h1));figure(2);subplot(1,2,1);zplane(num1,den1);subplot(1,2,2);zplane(num2,den2);Rp,rs增大,chebyshev2的幅频响应图会好wp=500;ws=600;rp=1;rs=60;fs=1000;wp=500*2*pi;ws=600*2*pi;不明白的是双线性变换法得倒的滤波器没有冲激不变法的好。
%要求设计一chebyshev2低通数字滤波器,wp=30hz,ws=40hz,rp=0.5,rs=40,fs=100hz. wp=30;ws=40;rp=0.5;rs=40;fs=100;wp=30*2*pi;ws=40*2*pi;[n,wn]=ellipord(wp,ws,rp,rs,'s');[z,p,k]=ellipap(n,rp,rs);[num,den]=zp2tf(z,p,k);[num1,den1]=impinvar(num,den);[num2,den2]=bilinear(num,den,100);[h,w]=freqz(num1,den1);[h1,w1]=freqz(num2,den2);subplot(1,2,1);plot(w*fs/(2*pi),abs(h));subplot(1,2,2);plot(w1*fs/(2*pi),abs(h1));figure(2);subplot(1,2,1);zplane(num1,den1);subplot(1,2,2);zplane(num2,den2);Butterworth带通wp1=10;wp2=20;ws1=5;ws2=25;fs=100;rp=0.5;rs=50; wp1=10*2*pi;wp2=20*2*pi;ws1=5*2*pi;ws2=25*2*pi; wp=[wp1,wp2];ws=[ws1,ws2];[n,wn]=buttord(wp,ws,rp,rs,'s');bw=wn(2)-wn(1);wo=sqrt(wn(1)*wn(2)); [z,p,k]=buttap(n);[a,b,c,d]=zp2ss(z,p,k);[a1,b1,c1,d1]=lp2bp(a,b,c,d,wo,bw); [num,den]=ss2tf(a1,b1,c1,d1);[h3,w3]=freqs(num,den);figure(4);plot(w3/(2*pi),abs(h3));title(‘模拟滤波器的幅频特性’);figure(5);plot(w3/(2*pi),angle(h3));title(‘模拟滤波器的相频特性’);figure(6);zplane(num,den);title(‘模拟滤波器的零极点图’);[num1,den1]=impinvar(num,den,fs);[h2,w2]=freqz(num1,den1);[num2,den2]=bilinear(num,den,fs);[h1,w1]=freqz(num2,den2);figure(1);subplot(1,2,1);plot(w2*fs/(2*pi),abs(h2));title(‘冲击不变法数字滤波器的幅频特性’); subplot(1,2,2);plot(w1*fs/(2*pi),abs(h1));title(‘双线性法数字滤波器的幅频特性’); figure(2);subplot(1,2,1);zplane(num1,den1);title(‘冲击不变法数字滤波器的零极点图’); subplot(1,2,2);zplane(num2,den2);title(‘双线性法数字滤波器的零极点图’); figure(3);subplot(1,2,1);plot(w2*fs/(2*pi),angle(h2));title(‘冲击不变法数字滤波器的相频特性’); subplot(1,2,2);plot(w1*fs/(2*pi),angle(h1));title(‘双线性法数字滤波器的相频特性’);Chebshev1通带wp1=10;wp2=20;ws1=5;ws2=25;fs=100;rp=0.5;rs=50; wp1=10*2*pi;wp2=20*2*pi;ws1=5*2*pi;ws2=25*2*pi; wp=[wp1,wp2];ws=[ws1,ws2];[n,wn]=cheb1ord(wp,ws,rp,rs,'s');bw=wn(2)-wn(1);wo=sqrt(wn(1)*wn(2));[z,p,k]=cheb1ap(n,rp);[a,b,c,d]=zp2ss(z,p,k);[a1,b1,c1,d1]=lp2bp(a,b,c,d,wo,bw);[num,den]=ss2tf(a1,b1,c1,d1);[num1,den1]=impinvar(num,den,fs);[h2,w2]=freqz(num1,den1);[num2,den2]=bilinear(num,den,fs);[h1,w1]=freqz(num2,den2);subplot(1,2,1);plot(w2*fs/(2*pi),abs(h2));subplot(1,2,2);plot(w1*fs/(2*pi),abs(h1));figure(1);subplot(1,2,1);zplane(num1,den1);subplot(1,2,2);zplane(num2,den2);Chebshev2通带wp1=10;wp2=20;ws1=5;ws2=25;fs=100;rp=0.5;rs=50; wp1=10*2*pi;wp2=20*2*pi;ws1=5*2*pi;ws2=25*2*pi; wp=[wp1,wp2];ws=[ws1,ws2];[n,wn]=cheb2ord(wp,ws,rp,rs,'s');bw=wn(2)-wn(1);wo=sqrt(wn(1)*wn(2));[z,p,k]=cheb2ap(n,rs);[a,b,c,d]=zp2ss(z,p,k);[a1,b1,c1,d1]=lp2bp(a,b,c,d,wo,bw);[num,den]=ss2tf(a1,b1,c1,d1);[num1,den1]=impinvar(num,den,fs);[h2,w2]=freqz(num1,den1);[num2,den2]=bilinear(num,den,fs);[h1,w1]=freqz(num2,den2);subplot(1,2,1);plot(w2*fs/(2*pi),abs(h2));subplot(1,2,2);plot(w1*fs/(2*pi),abs(h1));figure(1);subplot(1,2,1);zplane(num1,den1);subplot(1,2,2);zplane(num2,den2);椭圆通带wp1=10;wp2=20;ws1=5;ws2=25;fs=100;rp=0.5;rs=50; wp1=10*2*pi;wp2=20*2*pi;ws1=5*2*pi;ws2=25*2*pi; wp=[wp1,wp2];ws=[ws1,ws2];[n,wn]=ellipord(wp,ws,rp,rs,'s');bw=wn(2)-wn(1);wo=sqrt(wn(1)*wn(2));[z,p,k]=ellipap(n,rp,rs);[a,b,c,d]=zp2ss(z,p,k);[a1,b1,c1,d1]=lp2bp(a,b,c,d,wo,bw);[num,den]=ss2tf(a1,b1,c1,d1);[num1,den1]=impinvar(num,den,fs);[h2,w2]=freqz(num1,den1);[num2,den2]=bilinear(num,den,fs);[h1,w1]=freqz(num2,den2);subplot(1,2,1);plot(w2*fs/(2*pi),abs(h2));subplot(1,2,2);plot(w1*fs/(2*pi),abs(h1));figure(1);subplot(1,2,1);zplane(num1,den1);subplot(1,2,2);zplane(num2,den2);。