数字信号处理Matlab实现实例.docx

《数字信号处理Matlab实现实例.docx》由会员分享，可在线阅读，更多相关《数字信号处理Matlab实现实例.docx（20页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、实用标准文案数字信号处理 Matlab 实现实例第章离散时间信号与系统例 1-1 用 MATLAB 计算序列-20113和序列120-1的离散卷积。解 MATLAB 程序如下：a=-2 0 1 -1 3;b=1 2 0 -1;c=conv(a,b); M=length(c)-1; n=0:1:M;stem(n,c);xlabel(n); ylabel(幅度);图 1.1 给出了卷积结果的图形，求得的结果存放在数组 c 中为：-2-413151-3。例 1-2 用 MATLAB 计算差分方程当输入序列为时的输出结果。解 MATLAB 程序如下：N=41;a=0.8 -0.44 0.36 0.22

2、;b=1 0.7 -0.45 -0.6;x=1 zeros(1,N-1);精彩文档实用标准文案k=0:1:N-1;y=filter(a,b,x); stem(k,y) xlabel(n);ylabel(幅度)图 1.2 给出了该差分方程的前 41 个样点的输出，即该系统的单位脉冲响应。例 1-3 用 MATLAB 计算例 1-2 差分方程所对应的系统函数的DTFT。解 例 1-2 差分方程所对应的系统函数为：H (z) =0.8 - 0.44z-1 + 0.36z-2 + 0.02z-31+ 0.7z-1 - 0.45z-2 - 0.6z-3其DTFT 为H (e- jw ) =k=256;0

3、.8 - 0.44e- jw + 0.36e- j 2w + 0.02e- j 3w1+ 0.7e- jw - 0.45e- j 2w - 0.6e- j 3w用MATLAB 计算的程序如下：num=0.8 -0.44 0.36 0.02;den=1 0.7 -0.45 -0.6;w=0:pi/k:pi; h=freqz(num,den,w); subplot(2,2,1);精彩文档实用标准文案plot(w/pi,real(h);grid title(实部)xlabel(omega/pi);ylabel(幅度) subplot(2,2,2); plot(w/pi,imag(h);gridtit

4、le( 虚 部 ) xlabel(omega/pi);ylabel(Amplitude) subplot(2,2,3);plot(w/pi,abs(h);grid title(幅度谱)xlabel(omega/pi);ylabel(幅值) subplot(2,2,4); plot(w/pi,angle(h);gridtitle( 相 位 谱 ) xlabel(omega/pi);ylabel(弧度)精彩文档实用标准文案第章离散傅里叶变换及其快速算法例 2-1 对连续的单一频率周期信号 按采样频率=20 和N =16，观察其DFT 结果的幅度谱。采样，截取长度N 分别选N解 此时离散序列，即k=

5、8。用 MATLAB 计算并作图，函数 fft 用于计算离散傅里叶变换DFT，程序如下： k=8;n1=0:1:19;xa1=sin(2*pi*n1/k); subplot(2,2,1) plot(n1,xa1)xlabel(t/T);ylabel(x(n); xk1=fft(xa1);xk1=abs(xk1); subplot(2,2,2)stem(n1,xk1) xlabel(k);ylabel(X(k); n2=0:1:15;xa2=sin(2*pi*n2/k); subplot(2,2,3) plot(n2,xa2)xlabel(t/T);ylabel(x(n); xk2=fft(xa

6、2);xk2=abs(xk2);stem(n2,xk2) xlabel(k);ylabel(X(k);精彩文档subplot(2,2,4)实用标准文案计算结果示于图 2.1，(a)和(b)分别是 N=20 时的截取信号和 DFT 结果， 由于截取了两个半周期，频谱出现泄漏；(c) 和(d) 分别是N=16 时的截取信号和DFT 结果， 由于截取了两个整周期，得到单一谱线的频谱。上述频谱的误差主要是由于时域中对信号的非整周期截断产生的频谱泄漏。例 2-2 用 FFT 计算两个序列的互相关函数。解 用 MATLAB 计算程序如下： x=1 3 -1 1 2 3 3 1;y=2 1 -1 1 2 0

7、 -1 3;k=length(x); xk=fft(x,2*k);yk=fft(y,2*k); rm=real(ifft(conj(xk).*yk); rm=rm(k+2:2*k) rm(1:k);m=(-k+1):(k-1);精彩文档实用标准文案stem(m,rm)xlabel(m); ylabel(幅度);其计算结果如图 2.2 所示。例 2-3 计算两个序列的的互相关函数，其中x(n)=2 3 5 2 1 1 0 0 12 3 5 3 0 1 2 0 1 2y(n)=x(n-4)+e(n), e(n)为一随机噪声，在MATLAB 中可以用随机函数rand 产生解 用 MATLAB 计算程

8、序如下：x=2 3 5 2 1 -1 0 0 12 3 5 3 0 -1 -2 0 1 2;y=0 0 0 0 2 3 5 2 1 -1 0 0 12 3 5 3 0 -1 -2 0 1 2;k=length(y); e=rand(1,k)-0.5; y=y+e; xk=fft(x,2*k); yk=fft(y,2*k);rm=real(ifft(conj(xk).*yk); rm=rm(k+2:2*k) rm(1:k);m=(-k+1):(k-1);stem(m,rm)xlabel(m); ylabel(幅度);精彩文档 实用标准文案计算结果如图 2.3(a)，我们看到最大值出现在 m=4

9、处，正好是 y(n)对于x(n)的延迟。2. 3(b)是 x(n)的自相关函数，他和y(n)的区别除时间位置外，形状也略不同，这是由于y(n)受到噪声的干扰。虹邓军立膘m立嚼虔1实度151;1邸11lt:铅劲0gg心 占 I |I心J.j_|扫笠-1叶.,.,”及立-10.1忑(a) (b)阳 2 3 延迟厗列的互相关函获(a)和自相关函 救(b)精彩文档实用标准文案第 3 章 无限长单位脉冲响应(IIR)滤波器的设计方法例 3-1 设采样周期T=250s（采样频率fs=4kHz），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB 边界频率为fc=1kHz。B,A=butte

10、r(3,2*pi*1000,s);num1,den1=impinvar(B,A,4000); h1,w=freqz(num1,den1); B,A=butter(3,2/0.00025,s);num2,den2=bilinear(B,A,4000); h2,w=freqz(num2,den2); f=w/pi*2000;plot(f,abs(h1),-.,f,abs(h2),-); grid;xlabel(频率/Hz ) ylabel(幅值/dB)程序中第一个butter 的边界频率 21000，为脉冲响应不变法原型低通滤波器的边界频率；第二个 butter 的边界频率 2/T=2/0.000

11、25，为双线性变换法原型低通滤波器的边界频率.图 3.1 给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。脉冲响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差， 并且不存在传输零点。同时，也看到双线性变换法，在z=-1 即=或 f=2000Hz 处有一个三阶传输零点，这个三阶零点正是模拟滤波器在=处的三阶传输零点通过映射形成的。例 3-2设计一数字高通滤波器，它的通带为 400500Hz，通带内容许有 0.5dB 的波精彩文档实用标准文案动，阻带内衰减在小于 317Hz 的频带内至少为 19dB，采样频率为 1,000Hz。wc=2*1000*ta

12、n(2*pi*400/(2*1000); wt=2*1000*tan(2*pi*317/(2*1000); N,wn=cheb1ord(wc,wt,0.5,19,s);B,A=cheby1(N,0.5,wn,high,s); num,den=bilinear(B,A,1000); h,w=freqz(num,den);f=w/pi*500; plot(f,20*log10(abs(h);axis(0,500,-80,10);grid; xlabel()ylabel(幅度/dB)图 3.2 给出了MATLAB 计算的结果，可以看到模拟滤波器在=处的三阶零点通过高通变换后出现在=0（z=1）处，这

13、正是高通滤波器所希望得到的。例 3-3 设计一巴特沃兹带通滤波器，其dB 边界频率分别为f =110kHz 和 f =90kHz，在阻带f321= 120kHz 处的最小衰减大于dB，采样频率f =400kHz。s精彩文档实用标准文案w1=2*400*tan(2*pi*90/(2*400); w2=2*400*tan(2*pi*110/(2*400); wr=2*400*tan(2*pi*120/(2*400); N,wn=buttord(w1 w2,0 wr,3,10,s); B,A=butter(N,wn,s); num,den=bilinear(B,A,400); h,w=freqz(n

14、um,den);f=w/pi*200; plot(f,20*log10(abs(h);axis(40,160,-30,10);grid;xlabel(频率/kHz) ylabel(幅度/dB)图 3.3 给出了MATLAB 计算的结果，可以看出数字滤波器将无穷远点的二阶零点映射为z= 1 的二阶零点，数字带通滤波器的极点数是模拟低通滤波器的极点数的两倍。例 3-4 一数字滤波器采样频率 fs= 1kHz，要求滤除 100Hz 的干扰，其dB 的边界频率为 95Hz 和 105Hz，原型归一化低通滤波器为w1=95/500; w2=105/500;B,A=butter(1,w1, w2,stop

15、); h,w=freqz(B,A);f=w/pi*500; plot(f,20*log10(abs(h);axis(50,150,-30,10);grid;xlabel(频率/Hz)精彩文档 实用标准文案ylabel(幅度/dB)图 3.4 为 MATLAB 的计算结果105 卜 t- + - + - 一一：一： ： ：一一一：二二二二oll- - lll； IIIIIIIII- - r - - T- - T 一 一 ， 一 一F- - - - - - - - - - IIII, 1 1 T 1 1 T 1 1T 1 1 + 1 1 + 1 1T 1 1 + 1 1 + 1 1T 1 1 +

16、1 1 + 1 1r+ +E8IIII- -,- ， ，一一一厂. - - - -I I I II I II I- - - - - IIII lL_ _ I I II I I I I Ir T 1 1 T 1 1,IIII,l 与，笠-;)击讼丙泊沁1心110 心1让1叩1扣频率I Hz图 3.4 巴特沃兹带阻慈波客精彩文档实用标准文案第 4 章 有限长单位脉冲响应(FIR)滤波器的设计方法例 2用凯塞窗设计一 FIR 低通滤波器，低通边界频率，阻带衰减不小于 50dB。，阻带边界频率解 首先由过渡带宽和阻带衰减来决定凯塞窗的N 和图 4.1 给出了以上设计的频率特性，(a) 为 N=30 直

17、接截取的频率特性(b)为凯塞窗设计的频率特性。凯塞窗设计对应的MATLAB 程序为：wn=kaiser(30,4.55); nn=0:1:29;alfa=(30-1)/2;hd=sin(0.4*pi*(nn-alfa)./(pi*(nn-alfa); h=hd.*wn;h1,w1=freqz(h,1); plot(w1/pi,20*log10(abs(h1); axis(0,1,-80,10);grid;xlabel(归一化频率/ ) ylabel(幅度/dB)精彩文档实用标准文案例 4-2 利用雷米兹交替算法，设计一个线性相位低通FIR 数字滤波器，其指标为：通带边界频率f =800Hz，阻

18、带边界f =1000Hz，通带波动阻带最小衰减At=40dB，cr采样频率f =4000Hz。s解在 MATLAB 中可以用remezord 和 remez 两个函数设计，其结果如图4.2，MATLAB 程序如下： fedge=800 1000;mval=1 0;dev=0.0559 0.01;fs=4000;N,fpts,mag,wt=remezord(fedge,mval,dev,fs); b=remez(N,fpts,mag,wt); h,w=freqz(b,1,256); plot(w*2000/pi,20*log10(abs(h);grid;xlabel(频率/Hz) ylabel(

19、幅度/dB)函数 remezord 中的数组fedge 为通带和阻带边界频率，数组mval 是两个边界处的幅值，而数组dev 是通带和阻带的波动，fs 是采样频率单位为Hz。精彩文档实用标准文案第 5 章 数字信号处理系统的实现例 5-1 求下列直接型系统函数的零、极点，并将它转换成二阶节形式解 用 MATLAB 计算程序如下： num=1 -0.1 -0.3 -0.3 -0.2;den=1 0.1 0.2 0.2 0.5;z,p,k=tf2zp(num,den); m=abs(p);disp(零点);disp(z);disp(极点);disp(p);disp(增益系数);disp(k); s

20、os=zp2sos(z,p,k);disp(二阶节);disp(real(sos); zplane(num,den)输入到“num”和“den”的分别为分子和分母多项式的系数。计算求得零、极点增益系数和二阶节的系数：零点0.9615-0.5730-0.1443+0.5850i-0.1443-0.5850i极点0.5276+0.6997i0.5276-0.6997i-0.5776+0.5635i-0.5776-0.5635i增益系数1二阶节1.0000-0.3885-0.55091.00001.15520.65111.00000.28850.36301.0000-1.05520.7679系统函数

21、的二阶节形式为：精彩文档实用标准文案极点图见图 5.1。例 5-2 分析五阶椭圆低通滤波器的量化效应，其截止频率为0.4，通带纹波为0.4dB，最小的阻带衰减为 50dB。对滤波器进行截尾处理时，使用函数a2dT.m.。解 用以下MATLAB 程序分析量化效应clf; b,a=ellip(5,0.4,50,0.4);h,w=freqz(b,a,512); g=20*log10(abs(h);bq=a2dT(b,5);aq=a2dT(a,5); hq,w=freqz(bq,aq,512);gq=20*log10(abs(hq); plot(w/pi,g,b,w/pi,gq,r:); grid;a

22、xis(0 1 -80 5);xlabel(omega/pi);精彩文档实用标准文案ylabel(Gain, dB);legend(量化前,量化后); figurez1,p1,k1 = tf2zp(b,a);z2,p2,k2 = tf2zp(bq,aq); zplaneplot(z1,z2,p1,p2,o,x,*,+);legend(量化前的零点,量化后的零点,量化前的极点,量化后的极点);图 5.1（a）表示系数是无限精度的理想滤波器的频率响应（以实线表示）以及当滤波器系数截尾到 5 位时的频率响应（以短线表示）。由图可知，系数量化对频带的边缘影响较大，经系数量化后，增加了通带的波纹幅度，减

23、小了过渡带宽，并且减小了最小的阻带衰减。图 5. 1（b）给出了系数量化以前和系数量化以后的椭圆低通滤波器的零极点位置。由图可知，系数的量化会使零极点的位置与它们的理想的标称位置相比发生显著的改变。在这个例子中，靠近虚轴的零点的位置变动最大，并且移向靠它最近的极点的位置。只要对程序稍作改变就可以分析舍入量化的影响。为了研究二进制数量化效应对数字滤波器的影响，首先需要将十进制表示的滤波器系数转换成二进制数并进行量化，二进制数的量化既可以通过截尾法也可以通过舍入法实现。我们提供了如下的两个 MATLAB 程序：a2dT.m 和 a2dR.m，这两段程序分别将向量 d 中的每一个数按二进制数进行截尾

24、或舍入量化，量化的精度是小数点以后保留 b 位，量化后返回的向量为beq。function beq = a2dT(d,b)精彩文档实用标准文案% beq = a2dT(d,b) 将十进制数利用截尾法得到b 位的二进制数，%然后将该二进制数再转换为十进制数m=1; d1=abs(d);while fix(d1)0 d1=abs(d)/(2m); m=m+1;end beq=fix(d1*2b);beq=sign(d).*beq.*2(m-b-1);function beq=a2dR(d,b)% beq=a2dR(d,b)将十进制数利用舍入法得到b 位的二进制数%然后将该二进制数再转换为十进制数m

25、=1; d1=abs(d);while fix(d1)0 d1=abs(d)/(2m); m=m+1;end beq=fix(d1*2b+.5);beq=sign(d).*beq.*2(m-b-1);精彩文档实用标准文案第 7 章 多采样率信号处理例 7-1 在时域上显示一个,信号频率为 0.042的正弦信号，然后以抽取因子 3 降采样率，并在时域上显示相应的结果，比较两者在时域上的特点。解 用 MATLAB 计算程序如下：M=3; %down-sampling factor=3; fo=0.042;%signal frequency=0.042;%generate the input sin

26、usoidal sequence n=0:N-1;m=0:N*M-1;x=sin(2*pi*fo*m);%generate the down-sampling squence y=x(1:M:length(x); subplot(2,1,1)stem(n,x(1:N);title(输入序列);xlabel(时间/n);ylabel(幅度); subplot(2,1,2) stem(n,y);title(输出序列,抽取因子为,num2str(M);xlabel(时间/n); ylabel(幅度);和降采样率后的情况精彩文档图 7.1,信号频率为 0.042的正弦信号实用标准文案例 7-2 用汉明

27、窗设计一长度为 32 的线性相位QMF 滤波器组。解 采用 MATLAB 设计，调用fir2 函数设计公共低通滤波器，参数缺省即为汉明窗，程序如下：b1=fir2(31,0,0.4,0.5,0.55,0.6,1,1,1,1,0.06,0,0);for k=1:32b2(k)=(-1)(k-1)*b1(k); endH1z,w=freqz(b1,1,256); h1=abs(H1z); g1=20*log10(h1); H2z,w=freqz(b2,1,256); h2=abs(H2z); g2=20*log10(h2); figure(1);plot(w/pi,g1,-,w/pi,g2,-);

28、 axis(0,1,-100,10);gridxlabel(omega/pi);ylabel(幅度,dB); sum=h1.*h1+h2.*h2;d=10*log10(sum);figure(2)plot(w/pi,d);grid; xlabel(omega/pi);ylabel(误差,dB); axis(0,1,-0.3,0.3);图 7.2(a)是一个 N=32 的汉明窗设计结果，图中实线表示=的低通频响，虚线表示它的镜像= 计方法的分析/综合滤波器组的整个频响+。图 7.2 (b)是基于这种设。从这个图可见，重建误精彩文档实用标准文案差小于0.05dB。由于汉明窗设 计的频率响应在通带中近乎是平坦的，因此最大重建误差发生在这个滤波器的通带边界和过渡带内精彩文档