matlab实现数字滤波器
用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 filter函数用法
matlab filter函数用法MATLAB是一种在科学和数学领域广泛使用的编程语言和交互式环境。
在MATLAB中，filter函数是一个非常有用的函数，用于数字信号处理和数据滤波。
本文将详细介绍MATLAB的filter函数用法，以提供更好的理解和应用指导。
1、filter函数介绍filter函数是MATLAB中的数字滤波器函数，可以根据给定的数字滤波器系数和输入信号对信号进行滤波处理。
filter函数根据第一个参数（滤波器系数）和第二个参数（输入信号），对输入信号进行滤波，从而生成滤波后的信号。
该函数的语法如下：y = filter(b,a,x)其中，参数b是滤波器系数的向量，参数a是滤波器的反馈系数（因果滤波器的a参数为1），参数x是输入信号向量。
返回y是滤波后的信号向量。
2、使用滤波器系数在使用filter函数时，滤波器系数是必不可少的。
滤波器系数定义了数字滤波器处理输入信号的方式。
在计算机上，滤波器是使用数字算法来模拟模拟电路滤波器。
数字滤波器中最常见的类型是有限长冲激响应（FIR）滤波器和无限长冲激响应（IIR）滤波器。
这两种滤波器都在MATLAB中有广泛的应用。
当使用滤波器系数时，必须理解数字滤波器的工作原理。
数字滤波器的本质是将数字信号的特定频率区域转换为另一个频率区域。
因此，当进行数字信号滤波时，需要在输入信号上定义一个频率范围，然后使用滤波器系数来修改信号的频率分布。
3、理解滤波器系数的类型滤波器系数的类型包括传递函数系数和差分方程系数。
传递函数系数是数字滤波器的频域表现形式，是将输入信号转化为输出信号的函数。
差分方程系数是数字滤波器的时域表现形式，它实际上是将输入信号转化为输出信号的线性差分方程。
在filter函数中，b和a参数都是传递函数系数，因此在使用filter函数时必须了解这两个参数。
传递函数系数通常包括设计滤波器时选择的标准频率一些完成滤波器设计的附加参数。
由于B和A代表不同的信号属性，因此在计算数字滤波器时，必须分开使用它们。
数字滤波器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;。
MATLAB中的数字滤波器设计与实现
MATLAB中的数字滤波器设计与实现数字滤波器在信号处理中具有重要的作用，可用于去除噪声、滤波信号以及提取特定频率的成分。
MATLAB作为一种强大的数学软件，提供了多种数字滤波器设计和实现的工具，为工程师和科学家们提供了便捷而高效的解决方案。
本文将介绍MATLAB中数字滤波器的设计原理和实现方法，帮助读者更好地理解数字滤波器在实际应用中的重要性。
1. 数字滤波器的基本原理数字滤波器是一种能够改变信号频谱特性的系统，常用于消除噪声、去除不需要的频率成分或者提取感兴趣的信号成分。
数字滤波器分为FIR（有限长冲击响应）和IIR（无限长冲击响应）两种类型。
FIR滤波器的冲击响应为有限长序列，实现简单且稳定；IIR滤波器的冲击响应为无限长序列，具备更好的频率响应特性。
在MATLAB中，我们可以通过不同的函数和工具箱来设计这些数字滤波器。
2. FIR数字滤波器的设计与实现FIR数字滤波器的设计主要通过窗函数和频域采样进行。
MATLAB提供了一系列用于FIR滤波器设计的函数，如fir1、fir2等。
其中，fir1函数可以使用窗函数方法设计低通、高通、带通和带阻滤波器；fir2函数则可以实现自定义的频率响应。
具体设计步骤为：选择合适的窗函数、确定滤波器阶数和截止频率、生成滤波器系数。
设计完成后，可以通过filter函数将滤波器应用到目标信号上。
3. IIR数字滤波器的设计与实现IIR数字滤波器的设计方法主要有脉冲响应不变法和双线性变换法。
MATLAB提供了butter、cheby1、cheby2、ellip等函数来方便地实现IIR滤波器设计。
这些函数可以通过选择滤波器类型、阶数、截止频率等参数来生成相应的滤波器系数。
与FIR滤波器不同的是，IIR滤波器具有反馈结构，在MATLAB中可以使用filter函数来实现。
4. 数字滤波器的性能评估与优化正确评估和优化数字滤波器的性能对于滤波器的应用至关重要。
MATLAB提供了多种函数和工具箱来评估滤波器的频率响应、相位响应、时域响应等，如freqz、grpdelay、impz等。
基于MATLAB的IIR数字滤波器设计与仿真
基于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
数字滤波器的设计matlab
数字滤波器的设计可以使用MATLAB中的fdatool或者直接使用代码实现。
具体步骤如下：
1. 确定滤波器的类型（低通、高通、带通、带阻）以及参数（截止频率、通带、阻带宽度等）。
2. 使用MATLAB中的fdatool进行滤波器的设计，可以手动调节参数并进行实时分析和优化，或者选择预设的滤波器模板进行修改。
3. 点击Export将滤波器参数导出至matlab代码中。
4. 使用生成的代码进行滤波器的实现，例如：
```
fs = 1000; %采样率
fc = 100; %截止频率
[b, a] = butter(6, fc/(fs/2)); %设计6阶巴特沃斯低通滤波器input_signal = sin(2*pi*50*(0:1/fs:1)); %输入信号
output_signal = filter(b, a, input_signal); %滤波结果
```
其中，b和a分别为滤波器的分子和分母系数，可以根据不同的设计方法得到，input_signal为需要进行滤波的信号，output_signal为滤波后的结果。
数字滤波器的设计及其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的编程，在编程方面一点点的有了进步。