巴特沃斯带阻数字滤波器设计matlab程序及仿真图
设计巴特沃斯数字带通滤波器
设计巴特沃斯数字带通滤波器,要求通带范围为:0.25π rad ≤ω≤0.45π rad,通带最大衰减为3dB ,阻带范围为0≤ω≤0.15π rad 和0.55π rad ≤ω≤πrad ,阻带最小衰减为40dB 。
利用双线性变换设计,写出设计过程,并用MATLAB 绘出幅频和相频特性曲线。
设计思路及计算:(1)确定技术指标,求得数字边缘频率:10.25Pp rad ωπ=,20.45Pp rad ωπ=,3p a dB =10.15Ps rad ωπ=,20.55Ps rad ωπ=,40s a dB =(2)将数字带通滤波器的技术指标转换为模拟带通滤波器技术指标: 用双线性变换法,则2tan 2T ωΩ=,可得 112tan tan 0.1250.4142/2Pp Pp rad s T ωπΩ===222tan tan 0.2250.8541/2Pp Pp rad s T ωπΩ=== 112tan tan 0.0750.2401/2Ps Ps rad s T ωπΩ=== 222tan tan 0.275 1.1708/2Ps Ps rad s T ωπΩ=== (3)将带通滤波器的指标转换为模拟低通指标。
模拟低通归一化边界频率为:1Lp Ω=,()()2212221 1.9748Ps Pp Pp Ls Ps Pp Pp Ω-ΩΩΩ==ΩΩ-Ω(4)确定低通滤波器阶数N4020100.01s δ-==,()2211lg 1lg 10.01 6.76812lg 1.97482lg ss p N δ⎛⎫⎛⎫-- ⎪ ⎪⎝⎭⎝⎭≥==⎛⎫Ω ⎪ ⎪Ω⎝⎭取N =7。
(5)c c ΩΩ=Ω=1c Ω≈巴特沃兹模拟滤波器:(217)14711H (),()j K a k kk s p es p π++===-∏再由双线性变换即可得到所求。
代码实现:>> [N,Wn]=buttord([.25 .45],[.15 .55],3,40) N = 7 Wn = 0.2482 0.4525 >> [b,a]=butter(7,[.2482 .4525]) b = Columns 1 through 10 0.0001 0 -0.00070 0.0022 0 -0.0036 0 0.0036 0 Columns 11 through 15 -0.0022 0 0.00070 -0.0001 a = Columns 1 through 10 1.0000 -5.3094 16.2918 -34.7303 56.9401 -74.5112 80.0108 -71.1129 52.6364 -32.2233 Columns 11 through 15 16.1673 -6.4607 1.9827 -0.4217 0.0523>> [h,w]=freqz(b,a,100); >>subplot(211) >>h1=20*log10(abs(h)); >>plot(w/pi,h1); >>axis([0 1 -50 10]); >>subplot(212)>>plot(w/pi,angle(h))则滤波器传递函数为:24681012141234567890.00010.00070.00220.00360.00360.00220.00070.0001()1 5.309416.291834.730356.940174.511280.010871.112952.636432.223316.1673z z z z z z z H z z z z z z z z z z z ----------------⎛⎫-+- ⎪ ⎪+-+-⎝⎭=-+-+-+-+-+10111213146.4607 1.98270.42170.0523z z z z -----⎛⎫ ⎪ ⎪-+-+⎝⎭幅频和相频曲线:由上图可知,已满足设计要求0.10.20.30.40.50.60.70.80.91-50-40-30-20-1001000.10.20.30.40.50.60.70.80.91-4-2024。
用 MATLAB 设计巴特沃斯低通滤波器
用MATLAB 设计巴特沃斯低通滤波器1 巴特沃斯低通滤波器的特性一个理想低通滤波器的幅频特性如图3-80的阴影部分所示。
为了实现这个理想低通特性,需要在从0~ωC 的整个频带内增强增益,在ω>ωC 增益要降到0。
实际上,理想滤波器是不可能实现的。
图3-78是实际滤波器的幅频特性。
但是实际滤波器的特性愈接近理想特性愈好,巴特沃斯(Butterworth )滤波器就是解决这个问题的方法之一。
巴特沃斯滤波器以巴特沃斯函数来近似滤波器的系统函数,巴特沃斯的低通模平方函数为:221|()|1,2,,1(/)NC H j N j j ωωω==+ (3-138)式中以C ω是滤波器的电压-3dB 点或半功率点。
不同阶次的巴特沃斯滤波器特性如图3-79(a)所示。
4阶巴特沃斯滤波器的极点分布如图3-79(b)所示。
巴特沃斯滤波器幅频响应有以下特点:• 最大平坦性:在0=ω附近一段范围内是非常平直的,它以原点的最大平坦性来逼近理想低通滤波器。
• 通带、阻带下降的单调性。
这种滤波器具有良好的相频特性。
• 3dB 的不变性:随着N 的增加,频带边缘下降越陡峭,越接近理想特性。
但不管N 是多少,幅频特性都通过-3dB 点。
• 极点配置在半径为ωC 的圆上,并且均匀分布。
左半平面上的N 个极点是)(s H 的极点,右半平面上的N 个极点是)(s H -的极点。
2 巴特沃斯低通滤波器的实现为使巴特沃斯滤波器实用,我们必须能够实现它。
一个较好的方法是将巴特沃斯滤波器函数化成若干二阶节级联,其中每一节实现一对共轭复极点。
通过将极点以共轭复数的形式配对,对所有的每一个二阶节都具有实系数。
1图3-78 低通滤波器的幅频特性图3-80所示运算放大器电路为实现一对共轭极点提供了很好的方法。
电路的系统函数为202202121121122121)(1)11(1)(ωωω++=+++=s Qs C C R R s C R C R s C C R R s H (3-139)式中,ω0是S 平面原点与极点之间的距离,Q 被称为电路的“品质因数”,它提供了对响应峰值尖锐程度的一种度量。
利用MATLAB设计巴特沃斯低通数字滤波器
利用MATLAB设计巴特沃斯低通数字滤波器引言数字滤波器是数字信号处理中的重要组成部分,可以用于去除信号中的噪音和不需要的频率成分。
巴特沃斯滤波器是一种常见的数字滤波器,被广泛应用于信号处理领域。
本文将介绍如何利用MATLAB设计巴特沃斯低通数字滤波器,并给出详细的步骤和示例代码。
设计步骤利用MATLAB设计巴特沃斯低通数字滤波器主要包括以下步骤:1.设计滤波器的参数2.计算滤波器的传递函数3.绘制滤波器的幅频响应曲线4.通过频域图像观察滤波器的性能下面将分别介绍每个步骤的详细操作。
设计滤波器的参数巴特沃斯低通数字滤波器的参数包括截止频率和阶数。
截止频率决定了滤波器的通频带,阶数决定了滤波器的陡峭程度。
通过MATLAB的butter()函数可以方便地设计巴特沃斯低通数字滤波器。
该函数的参数为滤波器的阶数和截止频率。
示例代码如下:order = 4; % 阶数cutoff_freq = 0.4; % 截止频率[b, a] = butter(order, cutoff_freq);计算滤波器的传递函数通过设计参数计算得到滤波器的传递函数。
传递函数是一个复数,包括了滤波器的频率响应信息。
使用MATLAB的freqz()函数可以计算滤波器的传递函数。
该函数的参数为滤波器的系数b和a,以及频率取样点的数量。
示例代码如下:freq_points = 512; % 频率取样点数量[h, w] = freqz(b, a, freq_points);绘制滤波器的幅频响应曲线经过计算得到的传递函数能够提供滤波器的幅频响应信息。
通过绘制幅频响应曲线,可以直观地观察滤波器的频率特性。
使用MATLAB的plot()函数可以绘制滤波器的幅频响应曲线。
该函数的参数为频率点和传递函数的幅值。
示例代码如下:magnitude = abs(h); % 幅值plot(w/pi, magnitude);xlabel('归一化频率');ylabel('幅值');title('巴特沃斯低通数字滤波器幅频响应');通过频域图像观察滤波器的性能通过绘制滤波器的频域图像,可以直观地观察滤波器对不同频率的信号的响应情况。
(完整word版)巴特沃斯带阻数字滤波器设计matlab程序及仿真图 - 副本
fs=15000;T= 1/fs;rp=1;rs=40;wp1=0.11*pi;wp2=0.81*pi;ws1=0.31*pi;ws2=0.61*pi;%数字带阻滤波器技术指标wc1=(2/T)*tan(wp1/2);%频率预畸变wc2=(2/T)*tan(wp2/2);wr1=(2/T)*tan(ws1/2);wr2=(2/T)*tan(ws2/2);w0=sqrt(wc1*wc2);B=wc2-wc1;wp=1;%归一化通带截止频率ws=wp*(wr1*B) / (w0^2-wr1^2) ; %归一化阻带截止频率[N,wc]=buttord(wp,ws,rp,rs,'s')%求滤波器阶数和3dB截止频率[Z,P,K]=buttap(N)%设计模拟低通滤波器[Md,Nd]=zp2tf(Z,P,K)%将零极点形式转换为传输函数形式[M,N]=lp2bs(Md,Nd,w0,B)%对低通滤波器进行频率变换,转换为带阻滤波器[h,w]=freqs(M,N);%模拟带阻滤波器的幅频响应plot(w/(2*pi),abs(h));grid;xlabel('频率Hz');ylabel('幅度');title('模拟带阻滤波器');[b,a]=bilinear(M,N,15000)%对模拟滤波器双线性变换figure(1);freqz(b,a);[H,W]=freqz(b,a); %绘出频率响应;axis([0,1,-100,20]);figure(2);plot(W*fs/(2*pi),abs(H));grid on;xlabel('频率/Hz');ylabel('幅值');n=0:199;t=n/fs;x=sin(2*pi*400*t)+3*sin(2*pi*3000*t)+2*sin(2*pi*5000*t);figure(3);subplot(311);plot(t,x);axis([0,0.01,-5,5]);title('输入信号');grid on;y=filter(b,a,x);subplot(312);stem(y,'.');title('输出序列');grid on;ya=y*sinc(fs*(ones(length(n),1)*t-(n/fs)'*ones(1,length(t))));subplot(313);plot(t,ya);axis([0,0.01,-3,3]);title('输出波形');grid on;t=(0:100)/fs;figure(4)fs=1.5*10000;n=(0:100)/fs;f=sin(2*pi*400*t)+3*sin(2*pi*3000*t)+2*sin(2*pi*5000*t);y=fftfilt(b,x);[H1,f1]=freqz(f,[1]);[H2,f2]=freqz(y,[1]);f1=f1/pi*fs/2;f2=f2/pi*fs/2;subplot(2,1,1);plot(f1,abs(H1));title('输入信号的频谱');subplot(2,1,2);plot(f2,abs(H2));title('输出信号的频谱');基于Matlab 的带阻滤波器设计.10.20.30.40.50.60.70.80.91-800-600-400-2000N o r m a l i z e d Fre q u⨯π r a d /s a m p l e Ph a se(d e g r e e s )00.10.20.30.40.50.60.70.80.91-100-50N o r m a l i z e d Fr e q u⨯π r a d /s a m p l e M a g n i tu d e1000200030004000500060007000800000.20.40.60.811.21.4频率/Hz幅值00.0010.0020.0030.0040.0050.0060.0070.0080.0090.01-505输入信号020406080100120140160180200-22输出序列0.0010.0020.0030.0040.0050.0060.0070.0080.0090.01-202输出波形01000200030004000500060007000800050100150200输入信号的频谱010002000300040005000600070008000102030输出信号的频谱N =4wc =1.7947b =0.0186 -0.0410 0.1082 -0.1355 0.1810 -0.1355 0.1082 -0.0410 0.0186a =1.0000 -0.6707 -1.3750 0.5678 1.1964 -0.2996 -0.4631 0.0496 0.0762>。
巴特沃斯滤波器matlab实现
巴特沃斯滤波器matlab实现巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦,没有起伏,而在阻频带则逐渐下降为零。
在振幅的对数对角频率的波特图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大。
一阶巴特沃斯滤波器的衰减率为每倍频6分贝,每十倍频20分贝。
二阶巴特沃斯滤波器的衰减率为每倍频12分贝、三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。
巴特沃斯滤波器的振幅对角频率单调下降,并且也是唯一的无论阶数,振幅对角频率曲线都保持同样的形状的滤波器。
只不过滤波器阶数越高,在阻频带振幅衰减速度越快。
其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。
设计步骤如设计一个数字低通滤波器,其技术指标为:通带临界频率fp ,通带内衰减小于rp;阻带临界频率fs,阻带内衰减大于s;采样频率为FS1、将指标变为角频率wp=fp*2*pi;ws= fs*2*pi;2、将数字滤波器的频率指标{Wk}由wk=(2/T)tan(Wk/2)转换为模拟滤波器的频率指标{wk},由于是用双线性不变法设计,故先采取预畸变。
3、将高通指标转换为低通指标,进而设计高通的s域模型4、归一化处理由以上三式计算出N,查表可得模拟低通滤波器的阶数,从而由下式确定模拟高通滤波器的参数。
数字域指标变换成模拟域指标其程序为:fp = 400 fs= 300;Rp = 1; Rs = 20;wp =fp*2*pi;ws =fs*2*pi;FS=1000;T=1/FS;程序执行结果为:wp=2.5133e+003 ws=1.8850e+003 与实际计算结果相符。
数字滤波器的设计及其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仿真)
2012-5-9
13
图形:
0 Magnitude (dB) -20
-40
-60
0
100
200
300
400 500 600 Frequency (Hz)
700
800
900
1000
0 Phase (degrees) -20 -40 -60 -80 -100 0 100 200 300 400 500 600 Frequency (Hz) 700 800 900 1000
2012-5-9
24
得出结果: b= 0.2774 -0.0460 -0.0761 0.1008 -0.0571 0.1488 0.1262 -0.3949 -0.0783 0.0512 0.0270
a= 1.0000 -0.0185 0.2372 -0.0633 0.4144 0.1727 0.0393 0.0629 0.3588 0.0496 0.1016
2012-5-9
8
得出结果: N= 2 wn = 0.6630 b= 0.1578 -0.3155 a= 1.0000 0.6062 0.2373 0.1578
2012-5-9
9
图形:
0 Magnitude (dB) -50
-100
0
50
100
150
200 250 300 Frequency (Hz)
2012-5-9
18
设计一个带通Chebyshev II型数字滤波 器,满足:通带边界频率为100-200赫兹; 过渡带宽为30赫兹,通带波纹小于3dB, 阻带衰减大于30dB;采样频率为1000赫 兹。
2012-5-9
19
DSP实验4巴特沃斯滤波器的设计与实现(精)
DSP实验4巴特沃斯滤波器的设计与实现(精)实验四巴特沃斯数字滤波器的设计与实现1.数字滤波器的设计参数滤波器的4个重要的通带、阻带参数为:p f :通带截止频率(Hz ) s f :阻带起始频率(Hz )p R :通带内波动(dB ),即通带内所允许的最大衰减;s R :阻带内最小衰减设采样速率(即奈奎斯特速率)为N f ,将上述参数中的频率参数转化为归一化角频率参数:p ω:通带截止角频率(rad/s ),)2//(N p p f f =ω;s ω:阻带起始角频率(rad/s ),)2//(N s s f f =ω通过以上参数就可以进行离散滤波器的设计。
● 低通滤波器情况:采样频率为8000Hz ,要求通带截止频率为1500Hz ,阻带起始频率为2000Hz ,通带内波动3dB ,阻带内最小衰减为50dB ,则p ω=1500/4000,s ω=2000/4000,p R =3dB ,s R =50dB 。
● 高通滤波器情况:采样频率为8000Hz ,要求通带截止频率为1500Hz ,阻带起始频率为1000Hz ,通带内波动3dB ,阻带内最小衰减为65dB ,则p ω=1500/4000,s ω=1000/4000,p R =3dB ,s R =65dB 。
● 带通滤波器情况:采样频率为8000Hz ,要求通带截止频率为[800Hz ,1500Hz],阻带起始频率为[500Hz ,1800Hz],通带内波动3dB ,阻带内最小衰减为45dB ,则p ω=[800/4000,1500/4000],s ω=[500/4000,1800/4000],p R =3dB ,s R =45dB 。
● 带阻滤波器情况:采样频率为8000Hz ,要求通带截止频率为[800Hz ,1500Hz],阻带起始频率为[1000Hz ,1300Hz],通带内波动3dB ,阻带内最小衰减为55dB ,则p ω=[800/4000,1500/4000],sω=[1000/4000,1300/4000],p R =3dB ,s R =45dB 。
matlab巴特沃斯滤波器设计
数字IIR带阻滤波器的设计(基于巴特沃斯法)1、数字带阻IIR滤波器设计IIR数字滤波器在很多领域中有着广阔的应用。
与FIR数字滤波器相比,它可以用较低的阶数获得高选择性,所用存储单元少,经济而效率高,在相同门级规模和相同时钟速度下可以提供更好的带外衰减特性。
下面介绍一种设计实现IIR数字滤波器的方法。
设计步骤如下:步骤1:将数字滤波器H(z)的技术指标ωp和ωs,通过Ω=tan(ω/ 2)转变为模拟滤波器G(s)的技术指标Ωp和Ωs,作归一化处理后,得到ηp=1,ηs=Ωs/Ωp;步骤2:化解为模拟原型滤波器G(s)的技术指标;步骤3:设计模拟原型滤波器G(p);步骤4:将G(p)转换为模拟滤波器的转移函数G(s);步骤5:将G(s)转换成数字滤波器的转移函数H(z),s=(z–1)(z+1)。
所谓原型滤波器是指归一化的低通滤波器。
本节主要讨论通过IIR数字滤波器的原型转换设计法和IIR数字滤波器的直接设计方法来设计数字高通、带通及带阻滤波器其转换方法主要有3种:一是直接由模拟低通滤波器转换成数字高通、带通或带阻滤波器;二是先由模拟低通滤波器转换成模拟高通、带通或带阻滤波器,然后再把它转换成相应的数字滤波器;三是将模拟低通滤波器先转换成数字低通滤波器,再通过变量代换变换成高通、带通或带阻滤波器。
数字IIR 滤波器的设计结构如图:图1、IIR 滤波器的设计步骤2、带阻滤波器的设计规范:本设计中分别用h a (t)、H a (s)、Ha(j Ω)表示模拟滤波器的单位脉冲相应、系统函数、频率响应函数,三者的关系如下:Ha(s) = LT[ha(t)] =⎰∞∞-)(t ha e st -dtHa(j Ω)=FT[ha(t)] =⎰∞∞-)(t ha e t j Ω-dt可以用h a (t)、H a (s)、Ha(j Ω)中任一个描述模拟滤波器,也可以用线性常系数微分方程描述模拟滤波器。
但是设计模拟滤波器时,设计指标一般由幅频相应函数|Ha(jΩ)|给出,而模拟滤波器设计就是根据设计指标,求系统函数H a(s)。
MATLAB实现数字巴特沃斯高通IIR滤波器
目 录摘要摘要........................................................................ 2 Abstract ................................................................... 3 1 设计项目要求与说明设计项目要求与说明........................................................ 4 2 系统设计系统设计.................................................................. 4 2.1 设计思路设计思路............................................................ 4 2.2 设计方法对比设计方法对比........................................................ 5 2.3典型模拟滤波器比较典型模拟滤波器比较...................................................... 6 2.4 设计步骤设计步骤............................................................ 7 3 仿真程序的设计与调试仿真程序的设计与调试..................................................... 8 3.1 数字域指标变换成模拟域指标数字域指标变换成模拟域指标.......................................... 8 3.2 数字域频率进行预畸变数字域频率进行预畸变................................................ 8 3.3 模拟滤波器的设计模拟滤波器的设计.................................................... 8 3.4 模拟滤波器变成数字滤波器模拟滤波器变成数字滤波器........................................... 10 3.5 理论计算数字滤波器的仿真理论计算数字滤波器的仿真........................................... 13 4.4.程序调试中出现的问题程序调试中出现的问题程序调试中出现的问题..................................................... 14 5. 总结与体会总结与体会............................................................. 14 参考文献参考文献................................................................... 16 附录一附录一 总程序如下总程序如下总程序如下 ......................................................... 17 附录二附录二 设计数字滤波器函数总结设计数字滤波器函数总结设计数字滤波器函数总结 (19)摘要此报告重点介绍了用双线性不变法设计IIR 数字滤波器的基本流程,比较了各种设计方法的优缺点,总结了模拟滤波器的性能特征。
巴特沃斯滤波器matlab实现
巴特沃斯滤波器matlab实现巴特沃斯滤波器matlab实现巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦,没有起伏,而在阻频带则逐渐下降为零。
在振幅的对数对角频率的波特图上,从某一边界角频率开始,振幅随着角频率的增加而逐步减少,趋向负无穷大。
一阶巴特沃斯滤波器的衰减率为每倍频6分贝,每十倍频20分贝。
二阶巴特沃斯滤波器的衰减率为每倍频12分贝、三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。
巴特沃斯滤波器的振幅对角频率单调下降,并且也是唯一的无论阶数,振幅对角频率曲线都保持同样的形状的滤波器。
只不过滤波器阶数越高,在阻频带振幅衰减速度越快。
其他滤波器高阶的振幅对角频率图和低级数的振幅对角频率有不同的形状。
设计步骤如设计一个数字低通滤波器,其技术指标为:通带临界频率fp ,通带内衰减小于rp;阻带临界频率fs,阻带内衰减大于s;采样频率为FS1、将指标变为角频率wp=fp*2*pi;ws= fs*2*pi;2、将数字滤波器的频率指标{Wk}由wk=(2/T)tan(Wk/2)转换为模拟滤波器的频率指标{wk},由于是用双线性不变法设计,故先采取预畸变。
3、将高通指标转换为低通指标,进而设计高通的s域模型4、归一化处理由以上三式计算出N,查表可得模拟低通滤波器的阶数,从而由下式确定模拟高通滤波器的参数。
数字域指标变换成模拟域指标其程序为:fp = 400 fs= 300;Rp = 1; Rs = 20;wp =fp*2*pi;ws =fs*2*pi;FS=1000;T=1/FS;程序执行结果为:wp=2.5133e+003 ws=1.8850e+003 与实际计算结果相符。
MATLAB-巴特沃斯滤波器
MATLAB-巴特沃斯滤波器fp=0.2*pi;fs=0.3*pi;Rp=1;Rs=15;Fs=1000;Wp=2*pi*fp;Ws=2*pi*fs;Nn=256;n=(0:100-1);[N,Wn]=buttord(Wp,Ws,Rp,Rs,'s'); [b,a]=butter(N,Wn,'s');w=linspace(1,300,100)*2*pi;H=freqs(b,a,w);figure(1);plot(w/(2*pi),20*log10(abs(H)));title('巴特沃斯模拟滤低通波器幅频特性'); xlabel('频率/Hz');ylabel('幅度/db');[bz,az]=impinvar(b,a,Fs);figure(2);zplane(bz,az);title('冲击响应不变法零极点分布图') figure(3)Subplot(2,1,1);freqz(bz,az,Nn,Fs);grid on;title('巴特沃斯数字滤波器幅频特性'); xlabel('幅度/Hz');ylabel('频率/db'); Subplot(2,1,2);title('巴特沃斯数字滤波器相频特性'); xlabel('幅度/Hz');ylabel('频率/db');ws1=0.2*pi;wp1 = 0.35*pi;wp2 = 0.65*pi;wc = [wp1/pi,wp2/pi];B=wp1-ws1;N = ceil(8/0.15);n=0:N-1;%%rectangle窗window=rectwin(N);[h1,w]=freqz(window,1);subplot(121);stem(window);xlabel('n');title('rectangle窗函数');hn = fir1(N-1,wc,Boxcar(N));[h2,w]=freqz(hn,1,512);figure(1)subplot(122);plot(w/pi,20*log(abs(h2)/abs(h2(1)))); grid; xlabel('w/pi');ylabel('幅度(dB)');title('rectangle滤波器的幅度特性'); ;%%blackman窗window=hamming(N);[h1,w]=freqz(window,1);figure(2)subplot(121);stem(window);title('hamming窗函数');hn = fir1(N-1,wc,hamming(N));[h2,w]=freqz(hn,1,512);subplot(122);plot(w/pi,20*log(abs(h2)/abs(h2(1)))); grid; xlabel('w/pi');ylabel('幅度(dB)');title('hamming滤波器幅度特性'); %%hanning窗window=hanning(N);[h1,w]=freqz(window,1);figure(3)subplot(121);stem(window);xlabel('n');title('blackman窗函数');hn = fir1(N-1,wc,hanning(N));[h2,w]=freqz(hn,1,512);subplot(122);plot(w/pi,20*log(abs(h2)/abs(h2(1))));grid;xlabel('w/pi');ylabel('幅度(dB)');title('hanning滤波器的幅度特性'); %%blackman窗window=blackman(N);[h1,w]=freqz(window,1);figure(4)subplot(121);stem(window);xlabel('n');title('blackman窗函数');hn = fir1(N-1,wc,blackman(N));[h2,w]=freqz(hn,1,512);subplot(122);plot(w/pi,20*log(abs(h2)/abs(h2(1)))); grid; xlabel('w/pi');ylabel('幅度(dB)');title('blackman滤波器的幅度特性'); 第四题第五题。
基于MATLAB的巴特沃斯低通滤波器的设计
H (s)
1
(s p(1))( s p(2))(s p(n))
所以事实上 z 为空阵。上述零极点形式可以化为:
H (s)
sn
k bni s n1
b1s b0
其中 b0
n c
,令
c
1rad / s ,得到巴特沃斯滤波器归一化结果,如表
1
所示。
表 1 n 1 ~ 8阶的巴特沃斯滤波器系数
h(n)= ha(nT) 其中 T 是抽样周期。 如果令 Ha(s)是 ha(t)的拉普拉斯变换,H(z)为 h(n)的 z 变换,利用 抽样序列的 z 变换与模拟信号的拉普拉斯变换的关系,得:
X (z) zesT
1 T
Xa
k
(s
jks )
1 T
k
X
a
s
姓名: 班级: 学号: 时间:2011 年 6 月
设计题目
基于 MATLAB 的巴特沃斯低通滤波器的设计
设计要求
1. 通过实验加深对巴特沃斯低通滤波器基本原理的理解。 2.学习编写巴特沃斯低通滤波器的 MATLAB 仿真程序 3. 滤波器的性能指标如下:通带截止频率 fp=5kHz,通带最大衰减 p =2dB,阻带截止频率 fs=12kHz,阻带最小衰减 s =30dB
滤波器的性能指标如下:通带截止频率 fp=5kHz,通带最大衰减 p =2dB,
阻带截止频率 fs=12kHz,阻带最小衰减 s =30dB
3.1MATLAB 中所需函数
ATLAB 的信号处理工具箱提供了滤波器的函数 buttap、buttord、butter。 由[z,p,k] = buttap(n)函数可设计出 n 阶巴特沃斯低通滤波器原型,其传递函数 为
基于MATLAB的巴特沃斯滤波器.
数字信号处理课程设计2015年 6 月25 日目录一.设计目的: (3)二.设计要求: (3)三.设计内容: (4)3.1选择巴特涡斯低通数据滤波器及双线性变换法的原因 (4)3.2巴特沃思低通滤波器的基本原理 (4)3.3双线性变换法原理 (5)3.4数字滤波器设计流程图 (7)3.5数字滤波器的设计步骤 (7)四.用matlab实现巴特沃斯低通数字滤波器的仿真并分析 (9)4.1巴特沃斯低通数字滤波器技术指标的设置 (9)4.2用matlab实现巴特沃斯低通数字滤波器的仿真 (9)4.3波形图分析: (12)五.总结与体会 (13)六.附录参考文献 (14)2一.设计目的:该课程设计是测控技术与仪器专业的必修课,开设课程设计的目的使学生掌握数字信号处理的基本概念和基本理论,能够利用辅助工具进行FIR和IIR数字滤波器的设计,进行一维信号的频谱分析,并进行仿真验证。
加强实践教学环节,加强学生独立分析、解决问题的能力,培养学生动手能力和解决实际问题的能力,实现宽口径教育。
(1)理解低通滤波器的过滤方法。
(2)进一步熟悉低通滤波器的基本应用。
(3)用仿真工具matlab软件对设计的滤波器进行软件和硬件仿真。
(6)将对仿真结果进行比较,从而检验滤波器滤波性能的准确性。
二.设计要求:地震发生时,除了会产生地震波,还会由地层岩石在断裂、碰撞过程中所发生的震动产生次声波。
它的频率大约在每秒十赫兹到二十赫兹之间(可以用11Hz和15Hz的两个信号的和进行仿真,幅度可以分别设定为1、2)。
大气对次声波的吸收系数很小,因此它可以传播的很远,而且穿透性很强。
通过监测次声波信号可以监测地震的发生、强度等信息,因为自然界中广泛存在着各种次声波,这就对地震产生的次声波产生了干扰(可以用白噪声模拟,方差为5),需要采取一定的处理方法,才能检测到该信号,要求设计检测方案;并处理方法给出具体的软件(可以以51系列单片机、STM32F407、TMS320F28335或TMS320F6745为例)。
巴特沃斯带阻数字滤波器设计matlab程序及仿真图
fs=15000;T= 1/fs;rp=1;rs=40;wp1=0.11*pi;wp2=0.81*pi;ws1=0.31*pi;ws2=0.61*pi;%数字带阻滤波器技术指标wc1=(2/T)*tan(wp1/2);%频率预畸变wc2=(2/T)*tan(wp2/2);wr1=(2/T)*tan(ws1/2);wr2=(2/T)*tan(ws2/2);w0=sqrt(wc1*wc2);B=wc2-wc1;wp=1;%归一化通带截止频率ws=wp*(wr1*B) / (w0^2-wr1^2) ; %归一化阻带截止频率[N,wc]=buttord(wp,ws,rp,rs,'s')%求滤波器阶数和3dB截止频率[Z,P,K]=buttap(N)%设计模拟低通滤波器[Md,Nd]=zp2tf(Z,P,K)%将零极点形式转换为传输函数形式[M,N]=lp2bs(Md,Nd,w0,B)%对低通滤波器进行频率变换,转换为带阻滤波器[h,w]=freqs(M,N);%模拟带阻滤波器的幅频响应plot(w/(2*pi),abs(h));grid;xlabel('频率Hz');ylabel('幅度');title('模拟带阻滤波器');[b,a]=bilinear(M,N,15000)%对模拟滤波器双线性变换figure(1);freqz(b,a);[H,W]=freqz(b,a); %绘出频率响应;axis([0,1,-100,20]);figure(2);plot(W*fs/(2*pi),abs(H));grid on;xlabel('频率/Hz');ylabel('幅值');n=0:199;t=n/fs;x=sin(2*pi*400*t)+3*sin(2*pi*3000*t)+2*sin(2*pi*5000*t);figure(3);subplot(311);plot(t,x);axis([0,0.01,-5,5]);title('输入信号');grid on;y=filter(b,a,x);subplot(312);stem(y,'.');title('输出序列');grid on;ya=y*sinc(fs*(ones(length(n),1)*t-(n/fs)'*ones(1,length(t))));subplot(313);plot(t,ya);axis([0,0.01,-3,3]);title('输出波形');grid on;t=(0:100)/fs;figure(4)fs=1.5*10000;n=(0:100)/fs;f=sin(2*pi*400*t)+3*sin(2*pi*3000*t)+2*sin(2*pi*5000*t);y=fftfilt(b,x);[H1,f1]=freqz(f,[1]);[H2,f2]=freqz(y,[1]);f1=f1/pi*fs/2;f2=f2/pi*fs/2;subplot(2,1,1);plot(f1,abs(H1));title('输入信号的频谱');subplot(2,1,2);plot(f2,abs(H2));title('输出信号的频谱');00.10.20.30.40.50.60.70.80.91-800-600-400-200Normalized Frequency (⨯π rad/sample)P h a s e (d e g r e e s )00.10.20.30.40.50.60.70.80.91-100-50Normalized Frequency (⨯π rad/sample)M a g n i t u d e (d B )1000200030004000500060007000800000.20.40.60.811.21.4频率/Hz幅值00.0010.0020.0030.0040.0050.0060.0070.0080.0090.01-505输入信号020406080100120140160180200-22输出序列0.0010.0020.0030.0040.0050.0060.0070.0080.0090.01-202输出波形01000200030004000500060007000800050100150200输入信号的频谱010002000300040005000600070008000102030输出信号的频谱N =4wc =1.7947b =0.0186 -0.0410 0.1082 -0.1355 0.1810 -0.1355 0.1082 -0.0410 0.0186a =1.0000 -0.6707 -1.3750 0.5678 1.1964 -0.2996 -0.4631 0.0496 0.0762>。
巴特沃斯滤波器的设计与仿真
信号与系统课程设计论文摘要传统的数字滤波器的设计过程复杂,计算工作量大,滤波特性调整困难,影响了它的应用。
本文介绍了一种利用matlab提供的巴特沃斯滤波器设计函数“buttord”可求出所需的滤波阶数和3dB 截止频率的方法。
利用matlab设计滤波器设计函数,可以随时对比设计要求和滤波器特性调整参数,直观简便,极大的减轻了工作量,有利于滤波器设计的最优化。
关键词:巴特沃斯滤波器 Matlab 截止频率I信号与系统课程设计论文AbstractDesign for traditional digital filteris very complicated ,count also very complicated.it is very hard to readjust filtering character, influence apply,the main body of the book introduce use matlab provide butwosto design function“buttord” filteris 3dB end e matlab to design could compete reqire parameter burden,contribute to optimize.Keywords: filteris matlab end frequencyII信号与系统课程设计论文目录摘要 (Ⅰ)Abstrct第1章绪论 (1)1.1 课题背景 (1)第2章巴特沃斯滤波器的设计 (2)2.1 巴特沃斯滤波器阶数的选择 (2)2.2 巴特沃斯滤波器系数计算 (2)2.2.1 巴特沃斯低通滤波器系数计算 (2)2.2.2巴特沃斯高通滤波器系数计算 (3)2.2.3巴特沃斯带通滤波器系数计算 (3)2.2.4巴特沃斯带阻滤波器系数计算 (4)第3章巴特沃斯滤波器设计仿真 (5)3.1巴特沃斯滤波器设计仿真 (5)3.1.1 巴特沃斯低通滤波器实例仿真 (5)3.1.2巴特沃斯高通滤波器实例仿真 (6)3.1.3巴特沃斯带通滤波器实例仿真 (7)3.1.4巴特沃斯带阻滤波器实例仿真 (8)结论 (9)参考文献 (10)III第1章绪论1.1课题背景巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦,没有起伏,而在阻频带则逐渐下降为零。
基于Matlab的巴特沃斯IIR数字滤波器设计
பைடு நூலகம்
基于 Matlab 的巴特沃斯 IIR 数字滤波器设计
1.低通 Wp=0.2; Rp=1; Ws=0.3; Rs=15; [N,Wc]=buttord(Wp,Ws,Rp,Rs) [Bz,Az]=butter(N,Wc)
给出的滤波器的幅度和增益曲线。 参考曲线如下:
低通滤波器
ω/π
ω/π
幅度
增益 dB
给出的滤波器的幅度和增益曲线。 参考曲线如下:
带通滤波器
ω/π ω/π
增益 dB
幅度
4.带阻
Wp=[0.2,0.8]; Rp=1; Ws=[0.3,0.7]; Rs=15; [N,Wc]=buttord(Wp,Ws,Rp,Rs) [Bz,Az]=butter(N,Wc,'stop')
给出的滤波器的幅度和增益曲线。 参考曲线如下:
带阻滤波器
ω/π ω/π
某输入信号如下:
N=128; t=0:N-1; fs=1000; x1=sin(2*pi*50*t/fs); x2=sin(2*pi*150*t/fs); x3=sin(2*pi*250*t/fs); xn=x1+x2+x3;
对 xn 进行滤波处理: 1)设计低通滤波器,滤除 x2 和 x3,给出的滤波器的幅度、增益曲线和输出波形;(参 考过渡带:0.15π-0.25π) 2)设计高通滤波器,滤除 x1 和 x2,给出的滤波器的幅度、增益曲线和输出波形;(参 考过渡带:0.35π-0.45π) 3)设计带通滤波器,滤除 x1 和 x3,给出的滤波器的幅度、增益曲线和输出波形;(参 考过渡带:0.15π-0.25π和 0.35π-0.45π) 4)设计带阻滤波器,滤除 x2,给出的滤波器的幅度、增益曲线和输出波形;(参考 过渡带:0.15π-0.25π和 0.35π-0.45π) 5)调解 1)中的 as,使滤波器阶数 N=4; 6)调解 4)中的过渡带,改善输出波形;
基于matlab的数字滤波器设计及仿真
滤波器的指标
■ 滤波器的指标常在频域给出 ■ 数字滤波器的频率响应特性曲线:
H(ejw)= |H(ejw)|ej (w) |H(ejw)|为幅频特性函数,表示信号通过该 滤波器后各频率的衰减情况。
(w)为相频特性函数,表示信号通过滤波器 后在时间上的延时情况。
二、MATLAB环境下IIR数字滤波器的设计
四、两个有趣的实验
1、利用MATLAB进行声音滤波: 利用MATLAB读取wav格式的声音文件,对
信号进行离散傅里叶变换,得到其频谱图,给 信号加杂音,然后设计一个数字滤波器将杂音 滤除,得到原音,最后比较滤波器的滤波效果。
语音读取和加噪
■ clear;close all; ■ [Y,Fs,NBITS]=WAVREAD("WindowsXP.wav"); ■ n=length(Y); ■ f=0:Fs/n:Fs*(n-1)/n; ■ K=0:1/(n-1):1; ■ subplot(2,3,1);plot(K,Y);title("语音信号的时域波形"); ■ Y1=fft(Y,n); ■ subplot(2,3,3);plot(f,abs(Y1));title("语音信号的频谱"); ■ noise=0.02*sin(2*pi*4000*K); ■ subplot(2,3,2);plot(K(1:50),noise(1:50));title("噪声的时域波形"); ■ s=length(noise); ■ noise1=fft(noise,s); ■ subplot(2,3,5);plot(f,abs(noise1));title("噪声的频谱"); ■ A=Y"+noise; ■ subplot(2,3,4);plot(K,A);title("含噪语音的时域波形"); ■ A1=fft(A,s); ■ subplot(2,3,6);plot(f,abs(A1));title("含噪语音的频谱");
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
fs=15000;T= 1/fs;
rp=1;rs=40;
wp1=0.11*pi;wp2=0.81*pi;ws1=0.31*pi;ws2=0.61*pi;%数字带阻滤波器技术指标wc1=(2/T)*tan(wp1/2);%频率预畸变
wc2=(2/T)*tan(wp2/2);wr1=(2/T)*tan(ws1/2);wr2=(2/T)*tan(ws2/2);
w0=sqrt(wc1*wc2);B=wc2-wc1;
wp=1;%归一化通带截止频率
ws=wp*(wr1*B) / (w0^2-wr1^2) ; %归一化阻带截止频率
[N,wc]=buttord(wp,ws,rp,rs,'s')%求滤波器阶数和3dB截止频率
[Z,P,K]=buttap(N)%设计模拟低通滤波器
[Md,Nd]=zp2tf(Z,P,K)%将零极点形式转换为传输函数形式
[M,N]=lp2bs(Md,Nd,w0,B)%对低通滤波器进行频率变换,转换为带阻滤波器[h,w]=freqs(M,N);%模拟带阻滤波器的幅频响应
plot(w/(2*pi),abs(h));grid;
xlabel('频率Hz');ylabel('幅度');title('模拟带阻滤波器');
[b,a]=bilinear(M,N,15000)%对模拟滤波器双线性变换
figure(1);
freqz(b,a);[H,W]=freqz(b,a); %绘出频率响应;
axis([0,1,-100,20]);
figure(2);
plot(W*fs/(2*pi),abs(H));grid on;
xlabel('频率/Hz');
ylabel('幅值');
n=0:199;t=n/fs;
x=sin(2*pi*400*t)+3*sin(2*pi*3000*t)+2*sin(2*pi*5000*t);
figure(3);
subplot(311);plot(t,x);axis([0,0.01,-5,5]);
title('输入信号');grid on;
y=filter(b,a,x);
subplot(312);stem(y,'.');title('输出序列');grid on;
ya=y*sinc(fs*(ones(length(n),1)*t-(n/fs)'*ones(1,length(t))));
subplot(313);plot(t,ya);axis([0,0.01,-3,3]);title('输出波形');grid on;
t=(0:100)/fs;
figure(4)
fs=1.5*10000;
n=(0:100)/fs;
f=sin(2*pi*400*t)+3*sin(2*pi*3000*t)+2*sin(2*pi*5000*t);
y=fftfilt(b,x);
[H1,f1]=freqz(f,[1]);
[H2,f2]=freqz(y,[1]);
f1=f1/pi*fs/2;
f2=f2/pi*fs/2;
subplot(2,1,1);plot(f1,abs(H1));title('输入信号的频谱');
subplot(2,1,2);plot(f2,abs(H2));title('输出信号的频谱');
00.10.2
0.30.40.50.60.70.80.91
-800
-600-400-200
Normalized Frequency (⨯π rad/sample)
P h a s e (d e g r e e s )
00.10.2
0.30.40.50.60.70.80.91
-100
-50
Normalized Frequency (⨯π rad/sample)
M a g n i t u d e (d B )
100020003000
40005000600070008000
00.20.40.60.8
1
1.2
1.4
频率/Hz
幅值
00.0010.0020.0030.0040.0050.0060.0070.0080.009
0.01
-50
5输入信号
020406080100120140160180200
-2
2输出序列
0.0010.0020.0030.0040.0050.0060.0070.0080.009
0.01
-20
2输出波形
010002000300040005000600070008000
50100150
200输入信号的频谱
010002000300040005000600070008000
10
20
30输出信号的频谱
N =
4
wc =
1.7947
b =
0.0186 -0.0410 0.1082 -0.1355 0.1810 -0.1355 0.1082 -0.0410 0.0186
a =
1.0000 -0.6707 -1.3750 0.5678 1.1964 -0.2996 -0.4631 0.0496 0.0762
>。