《MATLAB在信号处理中的应用课件》由会员分享,可在线阅读,更多相关《MATLAB在信号处理中的应用课件(42页珍藏版)》请在金锄头文库上搜索。
1、MATLABMATLAB在信号处理中的应用1MATLAB在信号处理中的应用4.1 4.1 信号及其表示4.1.1连续时间信号的表示 连续时间信号:时间变化连续。如y=x(t) 离散时间信号( (序列) ):时间离散,如x(nT)=x(t)|t=nT.4.1.2工具箱中的信号产生函数函数名函数名功能功能函数名函数名功能功能sawtoothsawtooth产生锯齿波或三角波信号产生锯齿波或三角波信号pulstranpulstran产生冲激串产生冲激串squaresquare产生方波信号产生方波信号rectpulerectpule产生非周期的方波信号产生非周期的方波信号 sincsinc产生产生si
2、ncsinc函数波形函数波形tripulstripuls产生非周期的三角波信号产生非周期的三角波信号chirpchirp产生调频余弦信号产生调频余弦信号diricdiric产生产生DirichletDirichlet或周期或周期sincsinc函数函数gauspulsgauspuls产生高斯正弦脉冲信号产生高斯正弦脉冲信号gmonopulsgmonopuls产生高斯单脉冲信号产生高斯单脉冲信号vcovco电压控制振荡器电压控制振荡器2MATLAB在信号处理中的应用例:例:产生f=50Hz的锯齿波、三角波Fs=10000; %采样频率t=0:1/Fs:0.1; %采样间隔1/Fsf=50; %5
3、0Hzx1=sawtooth(2*pi*50*t,0);x2=sawtooth(2*pi*50*t,1);x3=sawtooth(2*pi*50*t,0.5);subplot(311); plot(t,x1);subplot(312); plot(t,x2);subplot(313); plot(t,x3)1 1、sawtoothsawtooth函数产生锯齿波或三角波用法:x=sawtooth(t)将产生周期为2的锯齿波。以02这个周期内为例,当t=0时,x=-1,当t=2时,x=1。x=sawtooth(t,width)width是0到1之间的标量。在0到2width区间内,x的值从-1线性
4、变化到1;在2width2区间内,x的值又从1线性变化到-1。sawtooth(t,1)和sawtooth(t)是等价的。3MATLAB在信号处理中的应用2、square函数函数产生矩形波例:例:产生50Hz占空比分别为20和50的矩形波。Fs=10000; %采样频率t=0:1/Fs:0.1; %采样间隔1/Fsf=50; %50Hzx1=square(2*pi*50*t,20);x2=square(2*pi*50*t,50);subplot(211); plot(t,x1);subplot(212); plot(t,x2);4MATLAB在信号处理中的应用3、sinc函数函数产生sinc波
5、形或sin(t)/(t)波形例:例:t=linspace(-10,+10,200);x=sinc(t);plot(t,x);5MATLAB在信号处理中的应用4.1.3离散时间信号的表示在MATLAB中,离散时间信号x(n)的表示:需用一个向量x表示序列幅值,用另一个等长的定位时间变量n,才能完整地表示一个序列。 例4-104-10绘制离散时间信号的棒状图。其中x(-1)=-1,x(0)=1,x(1)=2,x(2)=1,x(3)=0,x(4)=-1。MATLAB源程序为:n=-3:5;%n=-3:5;%定位时间变量x=0,0,-1,1,2,1,-1,0,0;x=0,0,-1,1,2,1,-1,0
6、,0;stem(n,x);grid;%stem(n,x);grid;%绘制棒状图line(-3,5,0,0);%line(-3,5,0,0);%画x x轴线xlabel(n);ylabel(xn)xlabel(n);ylabel(xn)运行结果如图4.11所示。图4.11离散时间信号图形line(起点横坐标,终点横坐标,起点纵坐标,终点纵坐标) 6MATLAB在信号处理中的应用4.1.44.1.4几种常用离散时间信号的表示1 1单位脉冲序列直接实现:x=zeros(1,N);x(1,nx=zeros(1,N);x(1,n0 0)=1;)=1;2 2单位阶跃序列直接实现:n=ns:nf;x=(n
7、-n0)=0;n=ns:nf;x=(n-n0)=0;x=zeros(1,5);x(1,3)=1;n=-2:2;stem(n,x)n=1:4;x=(n-2)=0;stem(n,x)7MATLAB在信号处理中的应用3 3实指数序列直接实现:n=ns:nf;x=a.n;n=ns:nf;x=a.n;4 4复指数序列直接实现:n=ns:nf;x=exp(sigema+jw)*n);n=ns:nf;x=exp(sigema+jw)*n);5 5正( (余) )弦序列直接实现:n=ns:nf;x=cos(w*n+sita);n=ns:nf;x=cos(w*n+sita);8MATLAB在信号处理中的应用4.
8、2 信号的基本运算4.2.14.2.1信号的相加与相乘信号的相加与相乘y(n)=x1(n)+x2(n)y(n)=x1(n)y(n)=x1(n)+x2(n)y(n)=x1(n)x2(n)x2(n)MATLABMATLAB实现:y=x1+x2;y=x1.*x2y=x1+x2;y=x1.*x24.2.24.2.2序列移位与周期延拓运算序列移位与周期延拓运算序列移位:y(n)=x(n-m)。MATLABMATLAB实现:y=x;ny=nx-my=x;ny=nx-m序列周期延拓:y(n)=x(n)M,MATLABMATLAB实现:ny=nxs:nxfny=nxs:nxf;y=x(mod(ny,M)+1)
9、y=x(mod(ny,M)+1)4.2.3 4.2.3 序列翻褶与序列累加运算序列翻褶与序列累加运算序列翻褶:y(n)=x(-n)。MATLAB可实现:y=fliplr(x)y=fliplr(x)序列累加的数学描述为:MATLABMATLAB实现:y=cumsum(x)y=cumsum(x)9MATLAB在信号处理中的应用4.2.4 两序列的卷积运算两序列卷积运算:MATLABMATLAB实现:y=conv(x1,x2)y=conv(x1,x2)。序列x1(n)x1(n)和x2(n)x2(n)必须长度有限。4.2.5 两序列的相关运算两序列相关运算:。MATLABMATLAB实现:y=xcor
10、r(x1,x2)y=xcorr(x1,x2)。10MATLAB在信号处理中的应用已知离散信号x(n)和h(n),求y(n)=x(n)*h(n),并用图形表示。Nh=20;Nx=10;m=5;Nh=20;Nx=10;m=5;% %设定NxNx,NhNh和位移值mmn=0:Nh-1;h1=(0.9).n;n=0:Nh-1;h1=(0.9).n;% %产生h1(n)h1(n)h2=h1;h2=h1;nx=0:Nx-1;x1=ones(1,Nx);nx=0:Nx-1;x1=ones(1,Nx);% %产生x1(n)x1(n)x2=zeros(1,Nx+m);x2=zeros(1,Nx+m);fork=
11、m+1:m+Nxfork=m+1:m+Nx% %产生x2(n)=x1(n-m)x2(n)=x1(n-m)x2(k)=x1(k-m);x2(k)=x1(k-m);endend% %产生x2(n)x2(n)11MATLAB在信号处理中的应用y1=conv(x1,h1);%计算y1(n)=x1(n)*h1(n)y2=conv(x2,h2);%计算y2(n)=x2(n)*h2(n)subplot(3,2,1)stem(nx,x1,.)axis(03001.2),title(x1(n)%绘图(以下省略)12MATLAB在信号处理中的应用4.4 线性时不变系统4.4.1 4.4.1 系统的描述系统的描述系
12、统传递函数在Matlab中,传递函数用分子、分母两个多项式的系数表示,系数为降幂排列降幂排列。13MATLAB在信号处理中的应用线性时不变(LTI)系统的常用表示方法1、传递函数表示法、传递函数表示法在Matlab中,传递函数用分子、分母两个多项式的系数表示,系数为降幂排列降幂排列。分子(Numerator) :B=b(1) b(2) b(m+1)分母(Denominator):A=a(1) a(2) b(n+1)例:num=1 0.2 1;den=1 0.5 1;14MATLAB在信号处理中的应用2、零极点模型表示法、零极点模型表示法在Matlab中,增益系数、零点向量、极点向量用三个列向量
13、三个列向量表示。增益系数(Gain) :k零点向量(Zero):z=z1 z2 zn极点向量(Pole):p=p1 p2 pnsys=zpk(z,p,k) %获得零-极点模型表达式k=3;z=234;p=567;sys=zpk(z,p,k,0.1)Zero/pole/gain:3(z-2)(z-3)(z-4)-(z-5)(z-6)(z-7)Samplingtime:0.115MATLAB在信号处理中的应用线性系统模型的变换函数函数名功能说明函数名功能说明ss2tfss2tf状态空间模型转换为传状态空间模型转换为传递函数模型递函数模型zp2tfzp2tf零极点增益模型转换为传递函零极点增益模型转
14、换为传递函数模型数模型ss2zpss2zp状态空间模型转换为零状态空间模型转换为零极点增益模型极点增益模型zp2sszp2ss零极点增益模型转换为状态空零极点增益模型转换为状态空间模型间模型ss2sosss2sos状态空间模型转换为二状态空间模型转换为二次分式模型次分式模型zp2soszp2sos零极点增益模型转换为二次分零极点增益模型转换为二次分式模型式模型tf2sstf2ss传递函数模型转换为状传递函数模型转换为状态空间模型态空间模型sos2tfsos2tf二次分式模型转换为传递函数模二次分式模型转换为传递函数模型型tf2zptf2zp传递函数模型转换为零传递函数模型转换为零极点增益模型极
15、点增益模型sos2zpsos2zp二次分式模型转换为零极点增二次分式模型转换为零极点增益模型益模型tf2sostf2sos传递函数模型转换为二传递函数模型转换为二次分式模型次分式模型sos2sssos2ss二次分式模型转换为状态空间模二次分式模型转换为状态空间模型型4.4.2系统模型的转换函数16MATLAB在信号处理中的应用 例4-18 4-18 求离散时间系统的零、极点向量和增益系数。在命令窗口输入:在命令窗口输入:num=2,3;den=1,0.4,1;num=2,3;den=1,0.4,1;num,den=eqtflength(num,den);num,den=eqtflength(n
16、um,den);%使长度相等z,p,k=tf2zp(num,den)z,p,k=tf2zp(num,den)屏幕显示为z=0-1.5000p=-0.2000+0.9798i-0.2000-0.9798ik=217MATLAB在信号处理中的应用4.5 线性时不变系统的响应4.5.14.5.1线性时不变系统的时域响应线性时不变系统的时域响应1 1连续LTILTI系统的响应2 2离散LTILTI系统的响应用MATLAB中的卷积函数conv()conv()来实现。用MATLAB中的卷积函数conv()conv()来实现。18MATLAB在信号处理中的应用4.5.2 LTI系统的单位冲激响应1.1.求连
17、续LTILTI系统的单位冲激响应函数impulse()impulse()格式:Y,T=impulse(sys)Y,T=impulse(sys)或impulse(sys)impulse(sys)功能:返回系统的响应Y和时间向量T,自动选择仿真的时间范围。其中sys可为系统传递函数、零极增益模型或状态空间模型。2.2.求离散系统的单位冲激响应函数dimpulse()dimpulse()格式:y,x=dimpulse(num,den)y,x=dimpulse(num,den)功能:返回项式传递函数的单位冲激响应y向量和时间状态历史记录x向量。19MATLAB在信号处理中的应用1、impulse函数函
18、数求连续连续系统的单位冲击响应。b=23;b=23;a=10.41;a=10.41;sys=tf(b,a)sys=tf(b,a) Transferfunction:Transferfunction:2s+32s+3-s2+0.4s+1s2+0.4s+1 yt=impulse(sys)yt=impulse(sys) ; ;plot(t,y)plot(t,y)20MATLAB在信号处理中的应用*2、impz函数函数求离散离散系统(数字滤波器)的单位冲击响应。(注:(注:Matlab 7.0不再支持不再支持dimpulse函数)函数)21MATLAB在信号处理中的应用4.5.3 时域响应的其它函数1
19、.1.求连续LTILTI系统的零输入响应函数initial()initial()格式:y,t,x=initial(a,b,c,d,x0)y,t,x=initial(a,b,c,d,x0)功能:计算出连续时间LTI系统由于初始状态x0所引起的零输入响应y。其中x为状态记录,t为仿真所用的采样时间向量。2.2.求离散系统的零输入响应函数dinitial()dinitial()格式:y,x,n=dinitial(a,b,c,d,x0)y,x,n=dinitial(a,b,c,d,x0)功能:计算离散时间LTI系统由初始状态x0所引起的零输入响应y和状态响应响应x,取样点数由函数自动选取。n为仿真所用
20、的点数。3.3.求连续系统的单位阶跃响应函数step()step()格式:Y,T=step(sys)Y,T=step(sys)功能:返回系统的单位阶跃响应Y和仿真所用的时间向量T,自动选择仿真的时间范围。其中sys可为系统传递函数(TF)、零极增益模型(ZPK)或状态空间模型(SS)。4.4.求离散系统的单位阶跃响应函数dstep()dstep()格式:y,x=dstep(num,den)y,x=dstep(num,den)功能:返回多项式传递函数G(z)=num(z)/den(z)表示的系统单位阶跃响应。22MATLAB在信号处理中的应用4.6线性时不变系统的频率响应1 1求模拟滤波器HaH
21、a( (s s) )的频率响应函数freqs()freqs()格式:H Hfreqs(B,A,W)freqs(B,A,W)功能:计算由向量W(rad/s)指定的频率点上模拟滤器系统函数Ha(s)的频率响应Ha(j),结果存于H向量中。 例4-314-31已知某模拟滤波器的系统函数求该模拟滤波器的频率响应。MATLAB源程序如下。B=1;A=12.61313.41422.61311;B=1;A=12.61313.41422.61311;W=0:0.1:2*pi*5;W=0:0.1:2*pi*5;freqs(B,A,W)freqs(B,A,W)图4.304.30模拟滤波器的频率响应23MATLAB
22、在信号处理中的应用 例4-324-32已知某滤波器的系统函数为求该滤波器的频率响应。MATLAB源程序为:B=100000001;B=100000001;A=1;A=1;freqz(B,A)freqz(B,A)该程序运行所绘出的幅频与相频性曲线如图4.31所示。图4.31滤波器幅度和相位曲线2 2求数字滤波器H H( (z z) )的频率响应函数freqz()freqz()格式:H=freqz(B,A,W)H=freqz(B,A,W)功能:计算由向量W(rad)指定的数字频率点上(通常指0,范围)在H(z)的频率响应H(ejw)。24MATLAB在信号处理中的应用3 3滤波函数filterfi
23、lter格式:y=filter(B,A,x)y=filter(B,A,x)功能:对向量x中的数据进行滤波处理,即差分方程求解,产生输出序列向量y。B和A分别为数字滤波器系统函数H(z)的分子和分母多项式系数向量。 例4-334-33设系统差分方程为MATLAB源程序为:B=1;A=1,-0.8;B=1;A=1,-0.8;N=0:31;x=0.8.n;N=0:31;x=0.8.n;y=filter(B,A,x);y=filter(B,A,x);subplot(2,1,1);stem(x)subplot(2,1,1);stem(x)subplot(2,1,2);stem(y)subplot(2,1
24、,2);stem(y)该程序运行所得结果如图4.32所示。,求该系统对信号的响应。图4.32系统对信号的响应25MATLAB在信号处理中的应用4.7傅里叶(Fourier)变换4.7.14.7.1连续时间、连续频率傅里叶变换连续时间、连续频率傅里叶变换4.7.2 4.7.2 连续时间、离散频率傅里叶级数连续时间、离散频率傅里叶级数正变换:逆变换:正变换:逆变换:26MATLAB在信号处理中的应用4.7.3 4.7.3 时间离散、连续频率序列傅里叶变换时间离散、连续频率序列傅里叶变换4.7.4 4.7.4 离散时间、离散频率离散傅里叶级数离散时间、离散频率离散傅里叶级数4.7.5离散时间、离散频
25、率离散傅里叶变换(DFT)正变换:逆变换:正变换:逆变换:正变换:逆变换:27MATLAB在信号处理中的应用1 1一维快速正傅里叶变换函数fftfft格式:X=fft(x,N)X=fft(x,N)功能:采用FFT算法计算序列向量x的N点DFT变换,当N缺省时,fft函数自动按x的长度计算DFT。当N为2整数次幂时,fft按基-2算法计算,否则用混合算法。2 2一维快速逆傅里叶变换函数ifftifft格式:x=ifft(X,N)x=ifft(X,N)功能:采用FFT算法计算序列向量X的N点IDFT变换。 例4-364-36用快速傅里叶变换FFT计算下面两个序列的卷积。,并测试直接卷积和快速卷积的
26、时间。图4.354.35快速卷积框图28MATLAB在信号处理中的应用MATLABMATLAB程序( (部分) ):%线性卷积xn=sin(0.4*1:15);xn=sin(0.4*1:15);%对序列x(n)赋值,M=15hn=0.9.(1:20);hn=0.9.(1:20);%对序列h(n)赋值,N=20yn=conv(xn,hn);yn=conv(xn,hn);%直接调用函数conv计算卷积%园周卷积L=pow2(nextpow2(M+N-1);L=pow2(nextpow2(M+N-1);Xk=fft(xn,L);Xk=fft(xn,L);Hk=fft(hn,L);Hk=fft(hn,
27、L);Yk=Xk.*Hk;Yk=Xk.*Hk;yn=ifft(Yk,L);yn=ifft(Yk,L);图4.36x(n),h(n)及其线性卷积波形29MATLAB在信号处理中的应用4.8 IIR数字滤波器的设计方法1.1.数字滤波器的频率响应函数幅度响应:相位响应:图4.37理想低通、高通、带通、带阻数字滤波器幅度特性30MATLAB在信号处理中的应用2.2.滤波器的技术指标幅度响应指标、相位响应指标图4.38数字低通滤波器的幅度特性通带要求: 阻带要求: 通带最大衰减:阻带最小衰减:31MATLAB在信号处理中的应用4.8.14.8.1冲激响应不变法冲激响应不变法2.MATLAB2.MATL
28、AB信号处理工箱中的专用函数impinvar()impinvar():格式:BZ,AZ=impinvar(B,A,Fs)BZ,AZ=impinvar(B,A,Fs)功能:把具有B,A模拟滤波器传递函数模型转换成采样频率为Fs(Hz)的数字滤波器的传递函数模型BZ,AZ。采样频率Fs的默认值为Fs=1。1.1.冲激响应不变法设计IIRIIR数字滤波器的基本原理: 例4-374-37MATLAB源程序如下:num=1;num=1;%模拟滤波器系统函数的分子den=1,sqrt(5),2,sqrt(2),1;den=1,sqrt(5),2,sqrt(2),1;%模拟滤波器系统函数的分母num1,de
29、n1=impinvar(num,den)num1,den1=impinvar(num,den)%求数字低通滤波器的系统函数程序的执行结果如下:num1=-0.00000.09420.21580.0311num1=-0.00000.09420.21580.0311den1=1.0000-2.00321.9982-0.76120.1069den1=1.0000-2.00321.9982-0.76120.106932MATLAB在信号处理中的应用MATLABMATLAB信号处理工具箱中的专用双线性变换函数bilinear()bilinear()格式:numd,dendnumd,dendbilinea
30、r(num,den,Fs)bilinear(num,den,Fs)功能:把模拟滤波器的传递函数模型转换成数字滤波器的传递函数模型。4.8.2双线性变换法双线性变换利用频率变换关系: 例4-384-38MATLAB源程序如下:num=1;num=1;%模拟滤波器系统函数的分子den=1,sqrt(3),sqrt(2),1;den=1,sqrt(3),sqrt(2),1;%模拟滤波器系统函数的分母num1,den1=bilinear(num,den,1)num1,den1=bilinear(num,den,1)%求数字滤波器的传递函数运算的结果如下:num1=0.05330.15990.15990
31、.0533num1=0.05330.15990.15990.0533den1=1.0000-1.33820.9193-0.1546den1=1.0000-1.33820.9193-0.154633MATLAB在信号处理中的应用4.8.3 IIR4.8.3 IIR数字滤波器的频率变换设计法数字滤波器的频率变换设计法1.IIR1.IIR数字滤波器的频率变换设计法的基本原理 根据滤波器设计要求,设计模拟原型低通滤波器,然后进行频率变换,将其转换为相应的模拟滤波器(高通、带通等),最后利用冲激响应不变法或双线性变换法,将模拟滤波器数字化成相应的数字滤波器。图4.39IIR数字滤波器MATLAB设计步骤
32、流程图34MATLAB在信号处理中的应用1 1MATLABMATLAB的典型设计利用在MATLAB设计IIR数字滤波器可分以下几步来实现(1)按一定规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标;(2)根据转换后的技术指标使用滤波器阶数函数,确定滤波器的最小阶数N和截止频率Wc;(3)利用最小阶数N产生模拟低通滤波原型;(4)利用截止频率Wc把模拟低通滤波器原型转换成模拟低通、高通、带通或带阻滤波器;(5)利用冲激响应不变法或双线性不变法把模拟滤波器转换成数字滤波器。 例4-394-39设计一个数字信号处理系统,它的采样率为Fs100Hz,希望在该系统中设计一个Butterworth
33、型高通数字滤波器,使其通带中允许的最小衰减为0.5dB,阻带内的最小衰减为40dB,通带上限临界频率为30Hz,阻带下限临界频率为40Hz。35MATLAB在信号处理中的应用MATLAB源程序设计如下:%把数字滤波器的频率特征转换成模拟滤波器的频率特征wp=30*2*pi;ws=40*2*pi;rp=0.5;rs=40;Fs=100;wp=30*2*pi;ws=40*2*pi;rp=0.5;rs=40;Fs=100;N,Wc=buttord(wp,ws,rp,rs,s);N,Wc=buttord(wp,ws,rp,rs,s);%选择滤波器的最小阶数Z,P,K=buttap(N);Z,P,K=b
34、uttap(N);%创建Butterworth低通滤波器原型A,B,C,D=zp2ss(Z,P,K);A,B,C,D=zp2ss(Z,P,K);%零极点增益模型转换为状态空间模型AT,BT,CT,DT=lp2hp(A,B,C,D,Wc);AT,BT,CT,DT=lp2hp(A,B,C,D,Wc);%实现低通向高通的转变num1,den1=ss2tf(AT,BT,CT,DT);num1,den1=ss2tf(AT,BT,CT,DT);%状态空间模型转换为传递函数模型%运用双线性变换法把模拟滤波器转换成数字滤波器num2,den2=bilinear(num1,den1,100);num2,den2
35、=bilinear(num1,den1,100);H,W=freqz(num2,den2);H,W=freqz(num2,den2);%求频率响应plot(W*Fs/(2*pi),abs(H);grid;plot(W*Fs/(2*pi),abs(H);grid;%绘出频率响应曲线xlabel(xlabel(频率/Hz);/Hz);ylabel(ylabel(幅值) )程序运行结果如图4.40所示。36MATLAB在信号处理中的应用2 2MATLABMATLAB的直接设计图4.39IIR4.39IIR数字滤波器MATLABMATLAB设计步骤流程图 例4-414-41试设计一个带阻IIR数字滤波
36、器,其具体的要求是:通带的截止频率:wp1650Hz、wp2850Hz;阻带的截止频率:ws1700Hz、ws2800Hz;通带内的最大衰减为rp0.1dB;阻带内的最小衰减为rs50dB;采样频率为Fs2000Hz。MATLABMATLAB源程序设计如下:wp1=650;wp2=850;ws1=700;ws2=800;rp=0.1;rs=50;Fs=2000;wp1=650;wp2=850;ws1=700;ws2=800;rp=0.1;rs=50;Fs=2000;wp=wp1,wp2/(Fs/2);ws=ws1,ws2/(Fs/2);wp=wp1,wp2/(Fs/2);ws=ws1,ws2/
37、(Fs/2);%利用Nyquist频率频率归一化N,wc=ellipord(wp,ws,rp,rs,z);N,wc=ellipord(wp,ws,rp,rs,z);%求滤波器阶数num,den=ellip(N,rp,rs,wc,stop);num,den=ellip(N,rp,rs,wc,stop);%求滤波器传递函数H,W=freqz(num,den);H,W=freqz(num,den);%绘出频率响应曲线plot(W*Fs/(2*pi),abs(H);grid;plot(W*Fs/(2*pi),abs(H);grid;xlabel(xlabel(频率/Hz);ylabel(/Hz);yl
38、abel(幅值) )该程序运行后的幅频响应曲线如图4.42所示。37MATLAB在信号处理中的应用4.9 FIR数字滤波器设计格式:w=boxcar(M)w=boxcar(M)功能:返回M点矩形窗序列。MATLABMATLAB信号处理工具箱中的窗函数法设计FIRFIR数字滤波器的专用命令fir1()fir1()。格式:B Bfir1(N,wc)fir1(N,wc)功能:设计一个具有线性相位的N阶(N点)的低通FIR数字滤波器,返回的向量B为滤波器的系数(单位冲激响应序列),其长度为N+1。4.9.14.9.1窗函数设计法窗函数设计法窗函数设计的基本原理:h h( (n n)=)=w w( (n
39、 n) )h hd d( (n n) )w(n)为窗函数,hd(n)理想数字滤波器的单位冲激响应。在MATLAB信号处理工具箱中为用户提供了BoxcarBoxcar(矩形)、BartletBartlet(巴特利特)、HanningHanning(汉宁)等窗函数,这些窗函数的调用格式相同。FIR数字滤波器的单位冲激响应h(n)满足偶(奇)对称h(n)=h(N-n-1)或h(n)=-h(N-n-1)FIR数字滤波器具有线性相位:或38MATLAB在信号处理中的应用 例4-434-43用矩形窗设计线性相位FIR低通滤波器。该滤波器的通带截止频率wc=pi/4,单位脉冲响h(n)的长度M=21。并绘出
40、h(n)及其幅度响应特性曲线。MATLAB源程序为:M=21;wc=pi/4;M=21;wc=pi/4; %理想低通滤波器参数n=0:M-1;r=(M-1)/2;n=0:M-1;r=(M-1)/2;nr=n-r+eps*(n-r)=0);nr=n-r+eps*(n-r)=0);hdn=sin(wc*nr)/pi./nr;hdn=sin(wc*nr)/pi./nr; %计算理想低通单位脉冲响应hd(n)ifrem(M,2)=0,hdn(r+1)=wc/pi;endifrem(M,2)=0,hdn(r+1)=wc/pi;end%M为奇数时,处理n=r点的0/0型wn1=boxcar(M);wn1=
41、boxcar(M);%矩形窗hn1=hdn.*wn1;hn1=hdn.*wn1;%加窗subplot(2,1,1);stem(n,hn1,.);line(0,20,0,0);subplot(2,1,1);stem(n,hn1,.);line(0,20,0,0);xlabel(n),ylabel(h(n),title(xlabel(n),ylabel(h(n),title(矩形窗设计的h(n);h(n);hw1=fft(hn1,512);w1=2*0:511/512;hw1=fft(hn1,512);w1=2*0:511/512;%求频谱subplot(2,1,2),plot(w1,20*log
42、10(abs(hw1)subplot(2,1,2),plot(w1,20*log10(abs(hw1)xlabel(w/pi),ylabel(xlabel(w/pi),ylabel(幅度(dB);title(dB);title(幅度特性(dB);(dB);程序运行结果如图4.44所示。39MATLAB在信号处理中的应用4.9.24.9.2频率抽样法频率抽样法1.1.频率抽样法的基本原理对所期望的滤波器的频率响应Hd(ejw)进行等间隔采样获得H(k),利用h(n)=IDFTH(k)求得FIR的单位冲激响应。2.MATLAB2.MATLAB信号处理工具箱中的频率抽样法专用函数命令fir2()fi
43、r2()格式:B Bfir2(N,F,A)fir2(N,F,A)功能:设计一个N阶的FIR数字滤波器,其频率响应由向量F和A指定,滤波器的系数(单位冲激响应)返回在向量B中,长度为N+1。向量F和A分别指定滤波器的采样点的频率及其幅值,F中的频率必须在0.0到1.0之间,1.0对应于采样频率的一半。它们必须按递增的顺序从0.0开始到1.0为结束。40MATLAB在信号处理中的应用 例4-474-47试用频率抽样法设计一个FIR低通滤波器,该滤波器的截止频率为0.5pi,频率抽样点数为33。MATLAB源程序为:N=32;N=32;F=0:1/32:1;F=0:1/32:1;%设置抽样点的频率,
44、抽样频率必须含0和1。A=ones(1,16),zeros(1,N-15);A=ones(1,16),zeros(1,N-15);%设置抽样点相应的幅值B=fir2(N,F,A);B=fir2(N,F,A);freqz(B);freqz(B);%绘制滤波器的幅相频曲线figure(2);stem(B,.);figure(2);stem(B,.);%绘制单位冲激响应的实部line(0,35,0,0);xlabel(n);ylabel(h(n);line(0,35,0,0);xlabel(n);ylabel(h(n);图4.494.49滤波器的频率响应和单位冲激响应序列41MATLAB在信号处理中的应用4.9.3 MATLAB的其它相关函数1最小二乘逼近法设计线性相位FIR滤波器函数fircls()fircls()2有限制条件的最小二乘逼近法设计低通和高通FIR数字滤波器函数fircls1()fircls1()3最小二乘逼近法设计线性相位FIR数字滤波器函数firls()firls()4升余弦FIR滤波器设计函数firrcos()firrcos()5Parks-McClellan优化等纹波FIR滤波器设计函数remez()remez()42MATLAB在信号处理中的应用