《MATLAB基础与应用教程(人民邮电出版社-蔡旭辉)第七章a.ppt》由会员分享，可在线阅读，更多相关《MATLAB基础与应用教程(人民邮电出版社-蔡旭辉)第七章a.ppt（65页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、第第 7章章 MATLAB在数字信号处在数字信号处理中的应用理中的应用17.1 时域离散信号的产生及时域处理时域离散信号的产生及时域处理时域离散信号用时域离散信号用x(n)表示，时间变量表示，时间变量n（表示采样位置）只（表示采样位置）只能取整数。因此，能取整数。因此，x(n)是一个离散序列，以后简称序列。是一个离散序列，以后简称序列。用一个向量用一个向量x不足以表示序列值不足以表示序列值x(n)。必须再用另一个等长。必须再用另一个等长的定位时间变量的定位时间变量n。x和和n同时使用才能完整地表示一个序同时使用才能完整地表示一个序列列，由于，由于n序列是按整数递增的，可简单地用其初值序列是按整
2、数递增的，可简单地用其初值ns决决定，因为它的终值定，因为它的终值nf取决于取决于ns 和和x的长度的长度length(x)，故可，故可写成：写成：n=ns:nf 或或n=ns:ns length(x)12单位脉冲序列单位脉冲序列(n-n0)的生成函数的生成函数impseqfunction x,n=impseq(n0,ns,nf)n=ns:nf;x=(n-n0)=0;单位阶跃序列单位阶跃序列u(n-n0)的生成函数的生成函数stepseqfunction x,n=stepseq(n0,ns,nf)n=ns:nf;x=(n-n0)=0;P173表表7.13例例7.1 序列的相加和相乘序列的相加和
3、相乘给出两个序列给出两个序列x1(n)和和x2(n)。x1=0,1,2,3,4,3,2,1,0;n1=-2:6;x2=2,2,0,0,0,-2,-2;n2=2:8;要求它们的和要求它们的和ya及乘积及乘积yp。解：编程的思路是把序列长度延拓到覆盖解：编程的思路是把序列长度延拓到覆盖n1和和n2的的范围，这样才能把两序列的时间变量对应起来，范围，这样才能把两序列的时间变量对应起来，然后进行对应元素的运算。然后进行对应元素的运算。4%第七章例第七章例7.1q701离散信号序列的相加和相乘离散信号序列的相加和相乘x1=0,1,2,3,4,3,2,1,0;ns1=-2;%给定给定x1及及ns1x2=2
4、,2,0,0,0,-2,-2;ns2=2;%给定给定x2及及ns2nf1=ns1+length(x1)-1;nf2=ns2+length(x2)-1;ny=min(ns1,ns2):max(nf1,nf2)%y(n)的时间变量的时间变量xa1=zeros(1,length(ny);xa2=xa1;%延拓序列初始化延拓序列初始化xa1(find(ny=ns1)&(ny=ns2)&(ny=0;%产生单位阶跃序列（产生单位阶跃序列（u(n-n1)）x1=(n-n1)=0-(n-n1-N)=0;%用阶跃序列差产生矩形序列用阶跃序列差产生矩形序列x2=(n=n1)&(n(N+n1);%用逻辑式产生矩形序
5、列用逻辑式产生矩形序列x3=exp(j*n*pi/8).*x2;%对复正弦序列加矩形窗（元素群乘）对复正弦序列加矩形窗（元素群乘）subplot(2,2,1);stem(n,x1,.);xlabel(n);ylabel(x1(n);%标注标注axis(n0,max(n),0,1);%定坐标范围定坐标范围subplot(2,2,3);stem(n,x2,.);xlabel(n);ylabel(x2(n);%标注标注axis(n0,max(n),0,1);%定坐标范围定坐标范围8subplot(2,2,2);stem(n,real(x3),.);xlabel(n);ylabel(x3(n)的实部的
6、实部);%标注标注line(n0,max(n),0,0);%画横轴画横轴axis(n0,max(n),-1,1);%定坐标范围定坐标范围subplot(2,2,4);stem(n,imag(x3),.);xlabel(n);ylabel(x3(n)的虚部的虚部);%标注标注line(n0,max(n),0,0);%画横轴画横轴axis(n0,max(n),-1,1);%定坐标范围定坐标范围结果结果输入序列起点输入序列起点:n0=-6输入序列长度输入序列长度:N=15输入位移输入位移:n1=3910例例7.3 序列的移位和周期延拓序列的移位和周期延拓已知已知，利用，利用MATLAB生成并图示生成
7、并图示 表示表示x(n)以以8为周为周期的延拓）和期的延拓）和 ，其中，其中N=24。解：周期延拓解：周期延拓方法方法1：利用矩阵乘法和冒号运算：利用矩阵乘法和冒号运算 x=1 2 3 4;y=x*ones(1,3);y1=(y(:)P176方法方法2：采用求余函数：采用求余函数mod y=x(mod(n,M)+1)可实现对可实现对x(n)以以M为周期的周期延拓。为周期的周期延拓。加加1是因为是因为MATLAB向量下标只能从向量下标只能从1开始。开始。11%第七章例第七章例7.3程序程序q703 序列序列x(n)的移位与周期延拓的移位与周期延拓%N:观察窗口长度观察窗口长度%M:序列序列x(n
8、)长度长度%m:移位样点数移位样点数clear;close allN=24;M=8;m=input(输入移位值：输入移位值：m=);if (m=N-M+1)%检验输入参数检验输入参数m是否合理是否合理 error(输入数据不在规定范围内！输入数据不在规定范围内！);end n=0:N-1;x1=(0.8).n;x2=(n=0)&(n=0;%产生输入信号产生输入信号x3(n)y3=filter(B,A,x3);%对对x3(n)的响应的响应x4=(n=0)&(n=0)&(n=nsx)&(n=nfx)=1)=x;%在对应的在对应的n处把处把xe赋值赋值xsubplot(2,1,1),stem(n,x
9、e,.),line(min(n(1),0),nf,0,0)%绘图绘图subplot(2,1,2),stem(n,y,.),line(min(n(1),0),nf,0,0)35由程序算出由程序算出nsy=-1留数：留数：r=-57.7581 204.7581极点：极点：p=0.7791 0.3209直接项：直接项：k=-150 -30代入代入得得如果不需要求出表达式，可用如果不需要求出表达式，可用filter函数求解。函数求解。36x(n)y(n)37例例7.9 离散时间傅里叶变换离散时间傅里叶变换取一个周期的正弦信号，作取一个周期的正弦信号，作8点采样，求它的连续频谱。然点采样，求它的连续频谱
10、。然后对该信号进行后对该信号进行N个周期延拓，再求它的连续频谱。把个周期延拓，再求它的连续频谱。把N无无限增大，比较分析其结果。限增大，比较分析其结果。解：离散时间傅里叶变换的定义：解：离散时间傅里叶变换的定义：MATLAB实现：设置一系列较密频率实现：设置一系列较密频率wi，求出一系列，求出一系列X(wi)X=x*exp(-j*n*w)。先把离散时间傅立叶变换写成一个子程序先把离散时间傅立叶变换写成一个子程序dtft.m，有了子程，有了子程序，本例就没有什么难度了。序，本例就没有什么难度了。38计算离散时间傅里叶变换的函数文件计算离散时间傅里叶变换的函数文件dtft.m：function X
11、=dtft(x,w)%计算离散时间傅立叶变换计算离散时间傅立叶变换%X 为在为在w频率点上的频率点上的DTFT数组数组%x 为输入离散序列，时间数组为增序整数，可用下标表示为输入离散序列，时间数组为增序整数，可用下标表示%w 为频率数组，由于为频率数组，由于exp(-j*n*w)以以2为周期，为周期，%所以其范围通常选为所以其范围通常选为-,或或0,2n=1:length(x);X=x*exp(-j*n*w);*注意：要恢复原来量纲，应把无量纲时间注意：要恢复原来量纲，应把无量纲时间n乘以采样周期乘以采样周期Ts，把无量纲频率把无量纲频率w（数字频率）乘以采样频率（数字频率）乘以采样频率Fs。
12、39%第七章例第七章例7.9程序程序q709%离散时间傅立叶变换离散时间傅立叶变换clear all;close all;disp(八点时间信号的离散时间傅立叶变换八点时间信号的离散时间傅立叶变换)x0=sin(2*pi*1:8/8)*5;%x0是是8点行向量点行向量dt=2*pi/8;%采样周期采样周期w=linspace(-2*pi,2*pi,1000)/dt;%w是是1000点行向量点行向量X0=dtft(x0,w)*dt;%求得频率响应求得频率响应X0subplot(3,1,1),plot(w,abs(X0),grid,shg%画图画图disp(重复重复N次的八点时间信号的离散时间傅立
13、叶变换次的八点时间信号的离散时间傅立叶变换)N=input(N=);%用键盘输入延拓周期数用键盘输入延拓周期数x1=reshape(x0*ones(1,N),1,N*length(x0);%延拓后的时域信号延拓后的时域信号x1X1=dtft(x1,w)*dt;%求求x1的频率响应的频率响应X1subplot(3,1,2),plot(w,abs(X1),grid,shg%画图画图disp(重复无穷次的八点信号的离散时间傅立叶变换重复无穷次的八点信号的离散时间傅立叶变换-傅立叶级数傅立叶级数)pause,X2=fft(x0)*dt;%离散傅立叶变换离散傅立叶变换w1=2*pi*0:length(x
14、0)-1/length(x0);%离散频点向量离散频点向量subplot(3,1,3),stem(-w1,w1,abs(X2),abs(X2),grid,shgaxis(min(w),max(w),0,max(abs(X2),grid,shg40 程序运行结果程序运行结果执行程序执行程序q709并按提示键入并按提示键入N=4，所得图形如图，所得图形如图7.10所示。所示。N取得愈大，其峰值愈大，宽度愈窄。当取得愈大，其峰值愈大，宽度愈窄。当N取得取得很大时，会出现内存不足的问题，这是用矩阵乘法做很大时，会出现内存不足的问题，这是用矩阵乘法做傅里叶变换的缺点。另外，因为那时峰值点处的宽度傅里叶变
15、换的缺点。另外，因为那时峰值点处的宽度很窄，也会出现所选频点对不上峰值点的问题。所以很窄，也会出现所选频点对不上峰值点的问题。所以对于对于N无限增大的情况，必须用无限增大的情况，必须用fft函数来求。这时用函数来求。这时用连续频谱也没有意义了。这里用同样的横坐标把几种连续频谱也没有意义了。这里用同样的横坐标把几种频谱进行对比，使读者更好地理解其关系。频谱进行对比，使读者更好地理解其关系。41图图7.1042例例7.10 时域采样频率与频谱混叠时域采样频率与频谱混叠分别以采样频率分别以采样频率fs=1000Hz，400Hz和和200Hz对对xa(t)进行等间隔采样，计算并图示三种采样频率下的采进
16、行等间隔采样，计算并图示三种采样频率下的采样信号及其幅频特性。样信号及其幅频特性。解：程序分别设定解：程序分别设定4种采样频率种采样频率fs=10kHz，1kHz，400Hz和和200Hz，对，对xa(t)进行采样，得到采样序列进行采样，得到采样序列xa(t)，xa1(n)，xa2(n)，xa3(n)，画出其幅度频谱。，画出其幅度频谱。采样时间区间均为采样时间区间均为0.1秒。为了便于比较，画出了幅秒。为了便于比较，画出了幅度归一化的幅频曲线，如图度归一化的幅频曲线，如图7.11所示。所示。43例例7.10 采样频率与频谱混叠采样频率与频谱混叠(续续)由于由于由以上关系式可见，由以上关系式可见
17、，采样信号的频谱函数是原模拟采样信号的频谱函数是原模拟信号频谱函数的周期延拓，延拓周期为信号频谱函数的周期延拓，延拓周期为2/T。如果如果以频率以频率f为自变量（为自变量（=2 f），则以采样频率），则以采样频率fs=1/T为延拓周期。为延拓周期。对频带限于对频带限于fc的模拟信号的模拟信号xa(t)，只，只有当有当fs2fc时，采样后时，采样后 才不会发生频谱混叠才不会发生频谱混叠失真。这就是著名的失真。这就是著名的采样定理采样定理 44%第七章例第七章例7.10程序程序q710%时域采样及其频谱时域采样及其频谱clear;close all;fs=10000;fs1=1000;fs2=40
18、0;fs3=200;%设置四种采样频率设置四种采样频率t=0:1/fs:0.1;%采集信号长度为采集信号长度为0.1秒秒A=444.128;a=50*sqrt(2)*pi;b=a;xa=A*exp(-a*t).*sin(b*t);k=0:511;f=fs*k/512;%由由wk=2k/512=2fT求得模拟频率求得模拟频率fXa=dtft(xa,2*pi*k/512);%近似模拟信号频谱近似模拟信号频谱T1=1/fs1;t1=0:T1:0.1;%采集信号长度为采集信号长度为0.1秒秒x1=A*exp(-a.*t1).*sin(b*t1);%1kHz采样序列采样序列x1(n)X1=dtft(x1
19、,2*pi*k/512);%x1(n)的的512点点dtftT2=1/fs2;t2=0:T2:0.1;%采集信号长度为采集信号长度为0.1秒秒x2=A*exp(-a.*t2).*sin(b.*t2);%400Hz采样序列采样序列x2(n)X2=dtft(x2,2*pi*k/512);%x2(n)的的512点点dtftT3=1/fs3;t3=0:T3:0.1;%采集信号长度为采集信号长度为0.1秒秒x3=A*exp(-a.*t3).*sin(b.*t3);%200Hz采样序列采样序列x3(n)X3=dtft(x3,2*pi*k/512);%x3(n)的的512点点dtft45figure(1);
20、subplot(2,2,1);plot(t,xa);axis(0,max(t),min(xa),max(xa);title(模拟信号模拟信号);xlabel(t(s);ylabel(Xa(t);line(0,max(t),0,0)subplot(2,2,2);plot(f,abs(Xa)/max(abs(Xa);title(模拟信号的幅度频谱模拟信号的幅度频谱);axis(0,500,0,1)xlabel(f(Hz);ylabel(|Xa(jf)|);subplot(2,2,3);stem(t1,x1,.);line(0,max(t1),0,0);axis(0,max(t1),min(x1),
21、max(x1)title(采样序列采样序列x1(n)（fs1=1kHz）);xlabel(n);ylabel(X1(n);f1=fs1*k/512;subplot(2,2,4);plot(f1,abs(X1)/max(abs(X1);title(x1(n)的幅度谱的幅度谱);xlabel(f(Hz);ylabel(|X1(jf)|);46figure(2);subplot(2,2,1);stem(t2,x2,.);line(0,max(t2),0,0);axis(0,max(t2),min(x2),max(x2);title(采样序列采样序列x2(n)（fs2=400Hz）);xlabel(n
22、);ylabel(X2(n);f=fs2*k/512;subplot(2,2,2);plot(f,abs(X2)/max(abs(X2);title(x2(n)的幅度谱的幅度谱);xlabel(f(Hz);ylabel(|X2(jf)|);subplot(2,2,3);stem(t3,x3,.);line(0,max(t3),0,0);axis(0,max(t3),min(x3),max(x3);title(采样序列采样序列x3(n)（fs3=200Hz）);xlabel(n);ylabel(X3(n);f=fs3*k/512;subplot(2,2,4);plot(f,abs(X3)/max
23、(abs(X3);title(x3(n)的幅度谱的幅度谱);xlabel(f(Hz);ylabel(|X3(jf)|)4748例例7.11 由离散序列恢复模拟信号由离散序列恢复模拟信号用时域内插公式用时域内插公式其中其中模拟用理想低通滤波器恢复的过程，观察恢复波形，模拟用理想低通滤波器恢复的过程，观察恢复波形，计算出最大恢复误差。计算出最大恢复误差。解：这个公式与卷积公式相像，可以用向量和矩阵乘解：这个公式与卷积公式相像，可以用向量和矩阵乘法来解决。法来解决。49例例7.11 由离散序列恢复模拟信号由离散序列恢复模拟信号xa=x*sinc(Fs*TNM)TNM=ones(length(n),1
24、)*t-n*T*ones(1,length(t);50%第七章例第七章例7.11程序程序q711%时域采样与重构时域采样与重构clear;close all;A=444.128;a=50*sqrt(2)*pi;b=a;for k=1:2 if k=1 Fs=400;elseif k=2 Fs=1000;end T=1/Fs;dt=T/3;%每个采样间隔每个采样间隔T上上g(t)取三个样点取三个样点 Tp=0.03;%重构时间区间为重构时间区间为0,0.03s t=0:dt:Tp;%生成序列生成序列t n=0:Tp/T;%生成序列生成序列n TNM=ones(length(n),1)*t-n*T
25、*ones(1,length(t);%生成生成TNM矩阵矩阵 x=A*exp(-a.*n*T).*sin(b*n*T);%生成模拟信号采样序列生成模拟信号采样序列x(n)xa=x*sinc(Fs*TNM);%内插公式内插公式 subplot(2,1,k);plot(t,xa);hold on axis(0,max(t),min(xa)-10,max(xa)+10);st1=sprintf(由由Fs=%d,Fs);st2=Hz的采样序列的采样序列x(n)重构的信号重构的信号;st=st1,st2;title(st),ylabel(xa(t);xo=A*exp(-a.*t).*sin(b*t);%
26、以以3Fs对原始模拟信号采样对原始模拟信号采样 stem(t,xo,.);line(0,max(t),0,0)emax2=max(abs(xa-xo)end51最大重构误差：最大重构误差：emax2=27.7015 emax2=9.943652例例7.12 梳状滤波器零极点和幅频特性梳状滤波器零极点和幅频特性梳状滤波器系统函数有如下两种类型。梳状滤波器系统函数有如下两种类型。FIR型：型：IIR型：型：freqz：数字滤波器频率特性计算和绘制函数：数字滤波器频率特性计算和绘制函数 zplane：H(z)的零的零-极点图绘制。极点图绘制。zplane(B,A)解：调用函数解：调用函数freqz和
27、和zplane 很容易写出程序很容易写出程序q712.m。53%第七章例第七章例7.12程序程序q712%梳状滤波器零极点和幅频特性梳状滤波器零极点和幅频特性clear;close allb=1,0,0,0,0,0,0,0,-1;a0=1;a1=1,0,0,0,0,0,0,0,-(0.8)8;a2=1,0,0,0,0,0,0,0,-(0.9)8;a3=1,0,0,0,0,0,0,0,-(0.98)8;H,w=freqz(b,a0);H1,w1=freqz(b,a1);H2,w2=freqz(b,a2);H3,w3=freqz(b,a3);54figure(1);subplot(2,2,1);z
28、plane(b,a0);title(FIR梳状滤波器零点图梳状滤波器零点图);subplot(2,2,2);zplane(b,a1);title(IIR梳状滤波器零、极点图梳状滤波器零、极点图,a=0.8);subplot(2,2,3);plot(w/pi,abs(H);title(FIR梳状滤波器幅频响应曲线梳状滤波器幅频响应曲线);ylabel(幅度幅度);xlabel(/);subplot(2,2,4);plot(w1/pi,abs(H1);title(IIR梳状滤波器幅频响应曲线梳状滤波器幅频响应曲线,a=0.8);ylabel(幅度幅度);xlabel(/);figure(2);su
29、bplot(2,2,1);zplane(b,a2);title(IIR梳状滤波器零、极点图梳状滤波器零、极点图,a=0.9);subplot(2,2,2);zplane(b,a3);title(IIR梳状滤波器零、极点图梳状滤波器零、极点图,a=0.98);subplot(2,2,3);plot(w2/pi,abs(H2);title(IIR梳状滤波器幅频响应曲线梳状滤波器幅频响应曲线,a=0.9);ylabel(幅度幅度);xlabel(/);subplot(2,2,4);plot(w3/pi,abs(H3);title(IIR梳状滤波器幅频响应曲线梳状滤波器幅频响应曲线,a=0.98);y
30、label(幅度幅度);xlabel(/)555657例例7.13 低通滤波及时域卷积定理低通滤波及时域卷积定理 输入信号输入信号 x(n)=cos(0.04 n)+cos(0.08 n)+cos(0.4 n)+0.3(n)，0n63 通过低通滤波器，计通过低通滤波器，计算滤波器对算滤波器对x(n)的响应输出的响应输出y(n)，并图示，并图示x(n)和和y(n)，观察滤波效果。，观察滤波效果。解：如前所述，只要求出解：如前所述，只要求出H(z)=B(z)/A(z)的分子和分的分子和分母多项式系数向量母多项式系数向量B和和A，则可调用滤波器直接，则可调用滤波器直接型型实现函数实现函数filter
31、对输入信号对输入信号x(n)进行滤波。进行滤波。y=filter(B,A,x)58%第七章例第七章例7.13程序程序q713%IIR滤波器实现及滤波器实现及FT的时域卷积定理验证的时域卷积定理验证clear;close all%产生输入信号产生输入信号x(n)n=0:255;N=4096;x=cos(0.04*pi*n)+cos(0.08*pi*n)+cos(0.4*pi*n);w=randn(size(x);%产生正态零均值噪声产生正态零均值噪声x=x+0.3*w;%求求H(z)分子分母多项式系数向量分子分母多项式系数向量B和和Ab=1,2,1;%(1+z-1)2 的展开系数的展开系数B=0
32、.0003738*conv(conv(b,b),b);%嵌套调用卷积函数嵌套调用卷积函数conva1=1,-1.2686,0.7051;a2=1,-1.0106,0.3583;a3=1,-0.9044,0.2155;A=conv(conv(a1,a2),a3);59%对对x(n)滤波滤波y=filter(B,A,x);%绘图绘图subplot(3,2,1);stem(x,.)axis(0,max(n)/4,min(x),max(x);line(0,max(n),0,0)title(输入信号输入信号x(n);xlabel(n);ylabel(x(n)subplot(3,2,3);stem(y,.
33、)axis(0,max(n)/4,min(y),max(y);line(0,max(n),0,0)title(输出信号输出信号 y(n);xlabel(n);ylabel(y(n)60X=fft(x,N);%计算计算x(n)的的N点离散点离散DFTY=fft(y,N);%计算计算y(n)的的N点离散点离散DFTk=0:N-1;f=2*k/N;subplot(3,2,2);plot(f,abs(X)title(输入信号输入信号x(n)的幅频曲线的幅频曲线);xlabel(/);ylabel(|FTx(n)|)axis(0,0.5,0,max(abs(X);subplot(3,2,4);plot(
34、f,abs(Y)title(输出信号输出信号y(n)的幅频曲线的幅频曲线);xlabel(/);ylabel(|FTy(n)|)axis(0,0.5,0,max(abs(Y);61H,f=freqz(B,A,N,whole);subplot(3,2,5);plot(f/pi,abs(H)title(滤波器幅频响应曲线滤波器幅频响应曲线);xlabel(/);ylabel(H幅度幅度)axis(0,0.5,0,max(abs(H);Ym=H.*X;subplot(3,2,6);plot(f/pi,abs(Ym)title(|FTx(n)FTh(n)|);xlabel(/);ylabel(Ym幅度
35、幅度)axis(0,0.5,0,max(abs(Ym);6263例例7.14 用符号运算工具箱解用符号运算工具箱解z变换问题变换问题解：无限长度时间序列的解：无限长度时间序列的z变换和逆变换和逆z变换都属于符变换都属于符号运算的范围。号运算的范围。MATLAB的的symbolic（符号运算）（符号运算）工具箱已提供了这种函数。如果读者已在计算机上工具箱已提供了这种函数。如果读者已在计算机上安装了这个工具箱，可以键入以下程序。安装了这个工具箱，可以键入以下程序。MATLAB程序程序q714.m 其特点是程序的开始要指定符号自变量其特点是程序的开始要指定符号自变量syms z n a N w0%规
36、定规定z,n,a为符号变量为符号变量64syms z n a N w0y1=an,Y1=simplify(ztrans(y1)y2=n,Y2=simplify(ztrans(y2)y3=n*an,Y3=simplify(ztrans(y3)y4=n*(n-1)/2,Y4=simplify(ztrans(y4)y5=exp(j*w0*n),Y5=simplify(ztrans(y5)y6=sin(w0*n),Y6=simplify(ztrans(y6)pauseX1=z/(z-1),x1=iztrans(X1)X2=-3*z-1/(2-5*z-1+2*z-2),x2=iztrans(X2)X3=z/(z-a),x3=iztrans(X3)X4=z/(z-1)2),x4=iztrans(X4)X5=z/(z-1)3),x5=iztrans(X5)X6=(1-z-N)/(1-z-1),x6=iztrans(X6)X7=z/(z-exp(j*w0),x7=iztrans(X7)65