模拟滤波器与数字滤波器的Matlab转换

实验四 IIR 滤波器的实现(IIR 滤波器是无穷脉冲响应数字滤波器的简称)一、 讲课目的熟悉用双线性变换法设计IIR 数字滤波器和方式。

二、讲课内容（一）大体概念1、 数字滤波器是指输入、输出均为数字信号，通过必然运算关系，改变输入信号所含频率成份的相对照例或滤除某些频率成份的器件。

数字滤波器与模拟滤波器概念相同，只是信号的形式和实现滤波的方式不同。

一样数字滤波器从现实的网络结构或从单位脉冲响应分类,能够分成无穷脉冲脉冲响应(IIR)滤波器和有限脉冲响应(FIR)滤波器。

2、 IIR 滤波器设计的要紧方式是先设计低通模拟滤波器，进行频率变换，将其转换为相应的高通、带通模拟滤波器，再将模拟滤波器转换为相应的数字滤波器。

3、 N 阶模拟滤波器系统函数的一样形式为 11101110...()()...()M M M M N N N N b s b s b s b B s H s a s a s a s a A s ----++++==++++ 4、 N 阶数字模拟滤波器系统函数的一样形式为12(1)012112(1)121...()()1...()M M M M N N N N b b z b z b z b z B z H z a z a z a z a z A z ------------+++++==+++++（二）求滤波器频率响应的函数一、freqs 函数：求模拟滤波器的频率响应freqs(b,a,w)计算由向量w （rad/s ）指定的频率点上（或W 直接为采样点数）频率响应，其中b 和a 别离为模拟滤波器系统函数H(s)的分子与分母。

freqs 函数将自动绘出幅频和相频曲线。

实验4-1：系统传输函数为220.20.31()0.41s s H s s s ++=++的模拟滤波器，绘出其幅频和相频谱。

a=[1,,1];b=[,,1];w=logspace(-1,1); %产生从10-1到101之间地50个等间距点,即50个频率点运行得幅频特性和相频特性如下：二、freqz函数：求数字滤波器的频率响应freqz(b,a,w)计算由向量w指定的频率点上模拟滤波器H(z)的频率响应，其中b和a别离为数字滤波器系统函数H(z)的分子与分母。

ＸX ＸX 大学XXXX 学院实验名称 IIR 数字滤波器的设计实验目的:加深理解ＩIR 数字滤波器的时域特性和频域特性,掌握ＩＩR 数字滤波器的设计原理与设计方法,以及I ＩＲ数字滤波器的应用。

实验内容:IIR 数字滤波器一般为线性移不变的因果离散系统,N 阶IIR 数字滤波器的系统函数可以表达为-1z 的有理多项式，即 -1-1-2-M =0012-1-2-N -112=1z +z +z ++z (z)==1+z +z ++z 1+zM j j M N Ni i b b b b b H a a a a ∑∑ 式中：系数i a 至少有一个非零。

对于因果II Ｒ数据滤波器,应满足M N ≤。

ＩＩR 数字滤波器的设计主要通过成熟的模拟滤波器设计方法来实现。

首先在频域将数字滤波器设计指标转换为模拟滤波器设计指标,然后将任意的模拟滤波器为原型模拟低通滤波器指标，根据模拟滤波器的设计指标来设计出模拟低通滤波器(s)LP H ,然后又(s)LP H 经过相应的复频域转换得到H(s),最后又H(s ）经过脉冲响应不变法或双线性变换法得到所需要的III Ｒ数字滤波器H （z)。

由此可见，IIR 数字滤波器设计的重要环节是模拟滤波器的设计。

设计模拟低通滤波器的主要方法有Butterwor ｔ、Ch ｅby ｓhev 、和椭圆等滤波器设计方法。

实验步骤1.Butterw ｏｒt 数字滤波器设计(1) Bu ｔt ｅrwort 滤波器是通带阻带都单调衰减的滤波器。

调用b ｕｔｔord 函数可以确定巴特沃斯滤波器的阶数，其格式为:[N,Omegac ］=bu ｔt ｏｒｄ(Ｏｍegap,Ｏｍe ｇas,Rp,As ，’s ’)。

其中,输入参数Rp,Ａs 分别为通带最大衰减和阻带最小衰减，以d Ｂ为单位；Om ｅg ａp,Ｏmegas 分别为通带截止频率和阻带截止频率,‘s ’说明所设计的是模拟滤波器。

输出参数为滤波器的阶数，Omeｇac为3ｄB截止频率。

基于MATLAB信号处理工具箱的数字滤波器设计与仿真•简介：传统的数字滤波器的设计过程复杂，计算工作量大，滤波特性调整困难，影响了它的应用。

本文介绍了一种利用MATLAB信号处理工具箱（Signal Processing Toolbox）快速有效的设计由软件组成的常规数字滤波器的设计方法。

给出了使用MATLAB语言进行程序设计和利用信号处理工具箱的FDATool工具进行界面设计的详细步骤。

利用MATLAB设计滤波器，可以随时对比设计要求和滤波器特性调整参数，直观简便，极大的减轻了工作量，有利于滤波器设计的最优化。

本文还介绍了如何利用MATLAB环境下的仿真软件Simulink对所设计的滤波器进行模拟仿真。

•关键字：数字滤波器 MATLAB FIR IIR•一、引言：在电力系统微机保护和二次控制中，很多信号的处理与分析都是基于对正弦基波和某些整次谐波的分析，而系统电压电流信号（尤其是故障瞬变过程）中混有各种复杂成分，所以滤波器一直是电力系统二次装置的关键部件【1】。

目前微机保护和二次信号处理软件主要采用数字滤波器。

传统的数字滤波器设计使用繁琐的公式计算，改变参数后需要重新计算，在设计滤波器尤其是高阶滤波器时工作量很大。

利用MATLAB信号处理工具箱（Signal Proce ssing Toolbox）可以快速有效的实现数字滤波器的设计与仿真。

1 数字滤波器及传统设计方法数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。

数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。

IIR数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配。

所以IIR滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。

FIR数字滤波器的单位脉冲响应是有限长序列。

基于MATLAB的IIR数字滤波器设计与仿真一、概述在现代数字信号处理领域中，数字滤波器扮演着至关重要的角色。

其通过对输入信号的特定频率成分进行增强或抑制，实现对信号的有效处理。

无限脉冲响应（IIR）数字滤波器因其设计灵活、实现简单且性能优良等特点，得到了广泛的应用。

本文旨在基于MATLAB平台，对IIR数字滤波器的设计与仿真进行深入研究，以期为相关领域的研究与应用提供有益的参考。

IIR数字滤波器具有无限长的单位脉冲响应，这使得其在处理信号时能够展现出优秀的性能。

与有限脉冲响应（FIR）滤波器相比，IIR滤波器在实现相同性能时所需的阶数更低，从而减少了计算复杂度和存储空间。

在需要对信号进行高效处理的场合，IIR滤波器具有显著的优势。

MATLAB作为一款功能强大的数学软件，提供了丰富的函数和工具箱，使得数字滤波器的设计与仿真变得简单而高效。

通过MATLAB，我们可以方便地实现IIR滤波器的设计、分析和优化，从而满足不同应用场景的需求。

本文将首先介绍IIR数字滤波器的基本原理和特性，然后详细阐述基于MATLAB的IIR数字滤波器的设计方法和步骤。

接着，我们将通过仿真实验验证所设计滤波器的性能，并对其结果进行分析和讨论。

本文将总结IIR数字滤波器设计与仿真的关键技术和注意事项，为相关领域的研究人员和工程师提供有益的参考和启示。

1. IIR数字滤波器概述IIR（Infinite Impulse Response）数字滤波器是数字信号处理中常用的一类滤波器，它基于差分方程实现信号的滤波处理。

与FIR （Finite Impulse Response）滤波器不同，IIR滤波器具有无限长的单位脉冲响应，这意味着其输出不仅与当前和过去的输入信号有关，还与过去的输出信号有关。

这种特性使得IIR滤波器在实现相同的滤波效果时，通常具有更低的计算复杂度，从而提高了处理效率。

IIR滤波器的设计灵活多样，可以根据不同的需求实现低通、高通、带通和带阻等多种滤波功能。

用MATLAB信号处理工具箱进行FIR滤波器设计的三种方法摘要介绍了利用MATLAB信号处理工具箱进行FIR滤波器设计的三种方法：程序设计法、FDATool设计法和SPTool设计法，给出了详细的设计步骤，并将设计的滤波器应用到一个混和正弦波信号，以验证滤波器的性能。

关键词 MATLAB，数字滤波器，有限冲激响应，窗函数，仿真1 前言数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。

根据其单位冲激响应函数的时域特性可分为两类：无限冲激响应（IIR）滤波器和有限冲激响应（FIR）滤波器。

与IIR滤波器相比，FIR的实现是非递归的，总是稳定的；更重要的是，FIR滤波器在满足幅频响应要求的同时，可以获得严格的线性相位特性。

因此，它在高保真的信号处理，如数字音频、图像处理、数据传输、生物医学等领域得到广泛应用。

2 FIR滤波器的窗函数设计法 FIR滤波器的设计方法有许多种，如窗函数设计法、频率采样设计法和最优化设计法等。

窗函数设计法的基本原理是用一定宽度窗函数截取无限脉冲响应序列获得有限长的脉冲响应序列，主要设计步骤为：(1) 通过傅里叶逆变换获得理想滤波器的单位脉冲响应hd(n)。

(2) 由性能指标确定窗函数W(n)和窗口长度N。

(3) 求得实际滤波器的单位脉冲响应h(n)， h(n)即为所设计FIR滤波器系数向量b(n)。

(4) 检验滤波器性能。

本文将针对一个含有5Hz、15Hz和30Hz的混和正弦波信号，设计一个FIR带通滤波器，给出利用MATLAB实现的三种方法：程序设计法、 FDATool设计法和SPTool设计法。

参数要求：采样频率fs=100Hz，通带下限截止频率fc1=10 Hz，通带上限截止频率fc2=20 Hz，过渡带宽6 Hz，通阻带波动0.01，采用凯塞窗设计。

2 程序设计法MATLAB信号处理工具箱提供了各种窗函数、滤波器设计函数和滤波器实现函数。

模拟有源滤波器设计的MATLAB实现课程设计摘要本文阐述了滤波器的基本概念,介绍了模拟有源滤波器的设计原理和逼近理论。

其中包括巴特沃思逼近、切比雪夫I型逼近、切比雪夫Ⅱ型逼近、椭圆函数逼近和贝塞尔逼近。

并研究了模拟有源滤波器的设计流程及性能测试。

综合了传统的硬件设计方法与软件编程技术,由MATLAB仿真出了各种滤波器逼近技术的幅频特性曲线并进行了实例分析。

对巴特沃思滤波器实例的研究仿真,由程序快速的得到了最小阶数和截止频率,取代了传统繁复的计算;方便的实现了由模拟低通滤波器向高通、低通和带阻滤波器的转换;对四运放复杂电路进行了设计仿真,通过求取其不同点的输出传递函数,模拟了二阶低通、高通、带通和带阻滤波器的幅频特性曲线并得到了较好的仿真结果。

关键词:模拟有源滤波器;逼近理论;幅频特性;MATLAB程序设计ABSTRACTThis paper describes the basic concepts of filters introduced analog active filter design principles and approximation theory. Including Butterworth approximation, Chebyshev type I approximation, Chebyshev Ⅱ type approximation, elliptic function approximation and Bezier approximation. The first author studied the analog active filter design process and performance testing.Combines the traditional hardware design methods and software programming techniques, the MATLAB simulation of a variety of filter approximation technique of amplitude-frequency characteristic curve and an illustrative example. Butterworth filter instance on simulation studies, by the program to quickly get the minimum order and the cutoff frequency, replacing the traditional complex calculations; convenientlyachieved by the analog low-pass filter to the high-pass, low-pass and band-stop filter the conversion; complex on quad op amp circuit design simulation, through its different points strike output transfer function to simulate the second-order low-pass, high pass, band pass and band-stop filter frequency characteristic curve and get better simulation results.Keywords:Analog and Active Filter;Theory of Approximation, Amplitude-frequency characteristics,MATLAB program design Compiling?and?organizing?dataAfter?you?have?established?the?purpose?of?the?report,?you?need?to compileandorganizetheinformationneededtosupportit.Thegathe ring?of?information?may?have?given?you?a?lot?of?materials,?but?you?ne ed?to?be?a?firm?editor?and?retain?only?the?essential?data?and?throw?o ut?the?rest.Consider?your?readers,?and?think?about?how?much?background?inform ation?they?will?need.Writing?the?reportA?report?consists?of?five?parts:?title,?introduction,?findings,?c onclusion?and?recommendations.1.?TitleThis?should?run?no?longer?than?one?line.2.?IntroductionThis?tells?the?reader?what?the?purpose?and?objective?of?the?repor t?are.?It?might?also?give?the?reader?some?background?information?on?t he?subject.The?purpose/objective/aim?of?this?report?is?toThis?report?aims?to/is?intended?to3.?FindingsThis?is?the?main?part?of?the?report.?It?tells?the?reader?what?you havefoundinyourinformationgathering.It?is?important?to?differentiate?between?fact?and?opinion.?Which? of?the?following?phrases?report?facts?and?which?report?opinions?We?found?that?could?be?both?fact?and?opinionIt?clearly?shows?that?factIt?was?found?thatfactWe?discovered?thatfactWe?observed?thatopinionThere?is?evidence?thatfact4.?ConclusionThis?part?tells?the?reader?about?the?results?of?the?report?based? on?the?findings.It's?concluded/decided/agreed/felt?thatIt?can?be?seen?thatIn?conclusionNo?conclusions?were?reached?regarding.关于一事未得出任何结论We?can?conclude?that5.?RecommendationsFinally,?recommendations?are?made?on?what?future?actions?need?to? be?taken.Based?on?our?findings,?we?would?recommend?thatIt?is?recommend-ed/proposed/suggested?that.It?seems?to?me?thattitle,?introduction,?findings,?conclusion?and?recommendations.1 Introduction滤波是信号处理的一种最基本而重要的技术,利用滤波可从复杂的信号中提取所需要的信号,抑制不需要的部分。

设计低通数字滤波器，要求在通带内频率低于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编程绘制各种窗函数的形状。

实验四IIR数字滤波器设计及软件实现一、实验指导1．实验目的（1）熟悉用双线性变换法设计IIR数字滤波器的原理与方法；（2）学会调用MATLAB信号处理工具箱中滤波器设计函数（或滤波器设计分析工具fdatool）设计各种IIR数字滤波器，学会根据滤波需求确定滤波器指标参数。

（3）掌握IIR数字滤波器的MATLAB实现方法。

（3）通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念。

2．实验原理设计IIR数字滤波器一般采用间接法（脉冲响应不变法和双线性变换法），应用最广泛的是双线性变换法。

基本设计过程是：①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标；②设计过渡模拟滤波器；③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。

MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。

第六章介绍的滤波器设计函数butter、cheby1 、cheby2 和ellip可以分别被调用来直接设计巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。

本实验要求读者调用如上函数直接设计IIR数字滤波器。

本实验的数字滤波器的MATLAB实现是指调用MATLAB信号处理工具箱函数filter对给定的输入信号x(n)进行滤波，得到滤波后的输出信号y(n）。

3. 实验内容及步骤（1）调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，该函数还会自动绘图显示st的时域波形和幅频特性曲线，如图10.4.1所示。

由图可见，三路信号时域混叠无法在时域分离。

但频域是分离的，所以可以通过滤波的方法在频域分离，这就是本实验的目的。

图10.4.1 三路调幅信号st的时域波形和幅频特性曲线（2）要求将st中三路调幅信号分离，通过观察st的幅频特性曲线，分别确定可以分离st 中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率。

西南石油大学实验报告一实验目的：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的编程，在编程方面一点点的有了进步。

开始↓读入数字滤波器技术指标↓将指标转换成归一化模拟低通滤波器的指标↓设计归一化的模拟低通滤波器阶数N 和截止频率↓模拟域频率变换,将H(P变换成模拟带通滤波器H(s ↓用双线性变换法将H(s转换成数字带通滤波器H(z ↓输入信号后显示相关结果求相应的幅频响应与相频响应↓50100150-202tx 1(tx1的波形50100150-202tx 2(tx2的波形50100150-202t x (t输入信号x 的波形10203040-0.01-0.00500.0050.01ty滤波器输出y 的波形clc;clear all ;结束%数字滤波器的技术指标Rp = 1; % 通带最大衰减Rs = 40;% 阻带最小衰减OmegaS1_1=350; % 通带截止频率OmegaS1_2=550;% 通带截止频率OmegaP1_1=400; % 阻带截止频率OmegaP1_2=500;% 阻带截止频率Fp=2000; % 抽样频率Wp1=2*pi*OmegaP1_1/Fp; % 模数频率变换Wp2=2*pi*OmegaP1_2/Fp;Ws1=2*pi*OmegaS1_1/Fp;Ws2=2*pi*OmegaS1_2/Fp;OmegaP1=2*Fp*tan(Wp1/2; % 非线性变换OmegaP2=2*Fp*tan(Wp2/2; % 非线性变换OmegaS1=2*Fp*tan(Ws1/2; % 非线性变换OmegaS2=2*Fp*tan(Ws2/2; % 非线性变换OmegaP0=sqrt(OmegaP1*OmegaP2;% 等效中心频率Bw=OmegaP2-OmegaP1; % 带通滤波器的通带宽度Eta_P0=OmegaP0/Bw; % 归一化处理Eta_P1=OmegaP1/Bw; % 归一化处理Eta_P2=OmegaP2/Bw; % 归一化处理Eta_S1=OmegaS1/Bw; % 归一化处理Eta_S2=OmegaS2/Bw; % 归一化处理Lemta_P_EquivalentLowPass=Eta_P2/(Eta_P2^2-Eta_P0^2; % 转换成低通参数Lemta_S1_EquivalentLowPass=-Eta_S1/(Eta_S1^2-Eta_P0^2; % 转换成低通参数Lemta_S2_EquivalentLowPass=Eta_S2/(Eta_S2^2-Eta_P0^2; % 转换成低通参数Lemta_S_EquivalentLowPass=min(Lemta_S1_EquivalentLowPass,Lemta_S2_EquivalentLowPass; % 取最小值% E求滤波器阶数[N, Wn]=cheb2ord(Lemta_P_EquivalentLowPass, Lemta_S_EquivalentLowPass, Rp, Rs,'s';% 滤波器设计[num1,den1]=cheby2(N,Rs,Wn,'s';[num2,den2]=lp2bp(num1,den1,OmegaP0,Bw;[num,den]=bilinear(num2,den2,Fp;[Z,P,K]=cheb1ap(N,Rp;w=linspace(1,1000,100*2*pi;[M1,N1]=zp2tf(Z,P,K; %将零极点形式转换为传输函数形式[M,N]=lp2bp(M1,N1,OmegaP0,Bw; %对低通滤波器进行频率变换转换为带通滤波器% 计算增益响应w = 0:pi/255:pi;h = freqz(num,den,w;g = 20*log10(abs(h;%绘制切比雪夫带通滤波器幅频特性figure;plot(w/pi,g;gridaxis([0 1 -60 5];xlabel('\频率/\pi'; ylabel('增益/dB'; title('切比雪夫II型带通滤波器幅频响应'; %Plot the poles and zeros[z,p,k]=tf2zp(num,den;figure;zplane(z,p; %绘制传输函数零极点title('?传输函数的零极点'f1=450;f2=600;t=0:0.0001:1x1=sin(2*pi*f1*t;x2=sin(2*pi*f2*t;x=x1+x2;figure;subplot(2,2,1%绘制x1的波形plot(x1;grid on;axis([0,50*pi,-3,3];xlabel('t';ylabel('x1(t';title('x1的波形';subplot(2,2,2%绘制x2的波形plot(x2;grid on;axis([0,50*pi,-3,3];xlabel('t';ylabel('x2(t';title('x2的波形';subplot(2,2,3%绘制输入x的波形plot(x;grid on;axis([0,50*pi,-3,3];xlabel('t';ylabel('x(t';title('输入信号x的波形'%X=fft(x;y=filter(num,den,x;%数字滤波器输出subplot(2,2,4; plot(real(y;grid on;axis([0,15*pi,-0.01,0.01];xlabel('t';ylabel('y';title('滤波器输出y的波形';附录:PPpppp5. 用双线性变换法设计IIR数字带通滤波器例21-3采用双线性变换法设计一个切比雪夫Ⅰ型数字带通滤波器,要求:通带wp1=0.3p,wp2=0.7p,Rp=1 dB;阻带ws1=0.2p,ws2=0.8p,As=20 dB解程序如下:wp1=0.4*pi;wp2=0.5*pi;ws1=0.35*pi;ws2=0.55*pi;Rp=1;As=40;T=0.0005;Fs=1/T;Omgp1=(2/T*tan(wp1/2;Omgp2=(2/T*tan(wp2/2;Omgp=[Omgp1,Omgp2];Omgs1=(2/T*tan(ws1/2;Omgs2=(2/T*tan(ws2/2;Omgs=[Omgs1,Omgs2];bw=Omgp2-Omgp1;w0=sqrt(Omgp1*Omgp2;bw=Omgs2-Omgs1;w0=sqrt(Omgs1*Omgs2; %[ZK(]模拟滤波器阻带带宽和中心频率[n,Omgn]=cheb2ord(Omgp,Omgs,Rp,As,'s' %计算阶数n和截止频率[z0,p0,k0]=cheb2ap(n,As; %设计归一化的模拟原型滤波器[n,Omgn]=cheb1ord(Omgp,Omgs,Rp,As,'s'[z0,p0,k0]=cheb1ap(n,Rp;ba1=k0*real(poly(z0;aa1=real(poly(p0;[ba,aa]=lp2bp(ba1,aa1,w0,bw;[bd,ad]=bilinear(ba,aa,Fs[H,w]=freqz(bd,ad;dbH=20*log10((abs(H+eps/max(abs(H;subplot(2,2,1,plot(w/2/pi*Fs,abs(H,'k';ylabel('|H|';title('幅度响应';axis([0,Fs/2,0,1.1];set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2];set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1];grid subplot(2,2,2,plot(w/2/pi*Fs,angle(H/pi*180,'k';ylabel('\phi';title('相位响应';axis([0,Fs/2,-180,180];set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2];set(gca,'YTickMode','manual','YTick',[-180,0,180];grid subplot(2,2,3,plot(w/2/pi*Fs,dbH;title('幅度响应( dB';axis([0,Fs/2,-40,5];ylabel('dB';xlabel('频率(\pi';set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2];set(gca,'YTickMode','manual','YTick',[-50,-20,-1,0];gridsubplot(2,2,4,zplane(bd,ad;axis([-1.1,1.1,-1.1,1.1];title('零极图';程序运行结果如下:n = 3Omgn =1.0e+003 * 1.0191 3.9252bd =0.0736 0.0000 -0.2208 0.0000 0.2208 -0.0000 -0.0736ad =1.0000 0.0000 0.9761 0.0000 0.8568 0.0000 0.2919 采用双线性变换法设计一个切比雪夫Ⅱ型数字带通滤波器,其它条件不变,则需要修改下面几句程序:bw=Omgs2-Omgs1;w0=sqrt(Omgs1*Omgs2; %[ZK(]模拟滤波器阻带带宽和中心频率[n,Omgn]=cheb2ord(Omgp,Omgs,Rp,As,'s' %计算阶数n和截止频率[z0,p0,k0]=cheb2ap(n,As; %设计归一化的模拟原型滤波器采用阻带截止频率来计算W0和BW,是因为切比雪夫Ⅱ型模拟低通原型是以阻带衰减As为主要设计指标的。