LDX第5次课-信号处理算法Matlab仿真

第5讲MATLAB数字信号处理应用西安邮电大学通信与信息工程学院数字信号处理与Matlb仿真MATLAB桌面复习：Matlab编程2M文件介绍1)M文件编辑器M文件编辑器提供了一个进行文本编辑和M文件调试的图形用户界面，如图所示。

M文件编辑器3

新建一个M文件有以下三种途径：①菜单操作（FileNewScript）②命令操作（editM文件名）③命令按钮（快捷键）

打开已有的M文件有四种途径：①菜单操作（FileOpen）②命令操作（editM文件名）③命令按钮（快捷键）④双击M文件名M文件介绍4图所示为M文件编辑器新建/打开快捷键。

M文件编辑器新建/打开快捷键52)函数函数文件由function语句引导，5个部分。

函数题头：指函数的定义行，是函数语句的第一行，在该行中定义函数名、输入变量列表及输出变量列表等。

HI行：指函数帮助文本的第一行，为该函数文件的帮助主题，当使用lookfor命令时，可以查看到该行信息。

帮助信息：提供了函数的完整的帮助信息，包括HI之后

至第一个可执行的行或空行为止的所有注释语句。

函数体：指函数代码段，也是函数的主体部分。

注释部分：指对函数体中各语句的解释和说明文本，注

释语句是以符号%引导的。6脚本是若干命令或函数的集合，用于执行特定的功能。在MATLAB命令窗口直接输入此文件的主文件名，MATLAB可顺序执行此文件内的所有命令，不需要在其中输入参数，也不需要给出输出变量来接收处理结果。脚本的操作对象为MATLAB工作空间内的变量，脚本式M文件运行产生的所有变量都是全局变量，这些变量都驻留在MATLAB基本工作空间内，通过使用clear命令或者关闭MATLAB系统可以清除这些变量。

3)脚本7

脚本文件可以理解为简单的M文件，脚本文件中的变量都是全局变量。函数文件是在脚本文件的基础之上多添加了函数定义行，其代码组织结构和调用方式与对应的脚本文件截然不同。函数文件是以函数声明行“function”作为开始，如果没有特别声明，函数文件中的变量都是局部变量。函数运行完毕之后，其定义的变量将从工作区间中清除。而脚本文件只是将一系列相关的代码结合封装，没有输入参数和输出参数，即不自带参数，也不一定要返回结果。而多数函数文件一般都有输入和输出变量，并返回结果。8一、常用序列主要内容引言二、序列的基本运算三、滤波四、离散付立叶变换（DFT）及应用引言：模拟信号与其采样序列对模拟信号

以采样周期T进行等间隔采样,

得到采样序列x(n)它的每一个样本点与相应处的模拟信号值相等，而在其它时间位置处的值为零。10引言：模拟信号与其采样序列11通常，用式（1）的序列作为离散信号的研究模型，由于抽掉了模拟信号中的主要物理量——时间。在数字域处理信号时,通常都是按序号逐次进行运算和存储，这就是抽象为序列模型的好处。引言：模拟信号与其采样序列12一、常用的序列

序列的表示方法完全地表示序列x(n)要用x和n两个向量，x表示序列的值，n表示它的位置或时间次序。例如x(n)={2,1.2,-1.4,3,1,4,3.1,7}（下划双线表示n=0处的采样点）用MATLAB表示时，用两个向量：

n=[-3,-2,-1,0,1,2,3,4];x=[2,1.2,-1.4,3,1,4,3.1,7];可见n是顺序增加的整数序列，可写成n=ns:nf,序列的长度为Lx=nf-ns+1。本例n=-3:4;Lx=8。13有限长序列和无限长序列序列位置向量的起点和终点都是有限数称为有限长序列序列位置向量的起点向-∞

延伸为无限长序列，简称左序列。序列位置向量的终点向∞延伸也为无限长序列，简称右序列。工程中遇到的序列都属于有限序列，有些情况下可看成右序列，但不会有左序列。一、常用的序列14一、常用的序列1．单位取样序列2．单位阶跃序列15矩形序列复数指数序列正余弦序列

x(n)=cos(ωn)+jsin(ωn)一、常用的序列16一、常用的序列17function[x,n]=impseq(np,ns,nf)%ns、nf分别为序列起始与终止位置,np为脉冲位置ifns>np|ns>nf|np>nf error（’输入位置参数不满足ns<=np<=nf’）else n=[ns:nf];x=[(n-np)==0];%核心语句end单位阶跃序列：把核心语句改为：x=[(n-np)>=0]例1：用Matlab生成脉冲序列18二、序列的基本运算序列的基本运算(a)序列相乘

(b)序列相加

(c) 倍乘(d)右移位(迟延)(e)左移位(提前)(f)复制191、序列相乘相加：将两个序列x1,x2的位置向量n1与n2对准，长度取两者之‘合’；形成一个新的包含n1与n2的公共位置向量n。其起点为两者中小的起点，终点为两者中大的终点，

n=min(ns1,ns2)::max(nf1,nf2);

将x1,x2补零，在n上延长，构成y1,y2;各个对应元素相乘或相加；y=y1.*y2;(或y=y1+y2;)二、序列的基本运算20%y(n)的长度n=min(min(n1),min(n2)):max(max(n1),max(n2));%y1,y2初始化y1=zeros(1,length(n));y2=y1;%把x1补零，扩展为具有y的长度的y1y1(find((n>=min(n1))&(n<=max(n1))))=x1;%把x2补零，扩展为具有y的长度的y2y2(find((n>=min(n2))&(n<=max(n2))))=x2;y=y1+y2;%序列相加.1、序列相乘相加：21例2：x1=[0,1,2,3,4,3,2,1,0];n1=-2:6; x2=[2,2,0,0,0,-2,-2];n2=2:8则 ny=-2:8; y1=[0,1,2,3,4,3,2,1,0,0,0]; y2=[0,0,0,0,2,2,0,0,0,-2,-2]; yadd=y1+y2=[0,1,2,3,6,5,2,1,0,-2,-2] ymult=y1.*y2=[0,0,0,0,6,9,0,0,0,0,0]1、序列相乘相加：22例2的计算图解23序列向量的值不变，只把位置向量右移k拍子程序[y,ny]=seqshift(x,nx,k)的核心语句：y=x;ny=nx+k左移时k取负值例3：把x=[3,7,-2,4];nx=-1:2;k=2代入则y=[3,7,-2,4];ny=1:42、序列移位y(n)=x(n-k)24将序列x以n=0点为中心，左右翻转。 例如：

x=[3,7,-2,4];nx=-1:2

则 y=[4,-2,7,3];ny=-2:1即x的值及其位置排序都要左右翻转子程序[y,ny]=seqfold(x,nx)的核心语句：

y=fliplr(x);

ny=-fliplr(nx)3、序列的折叠：y(n)=x(-n)25序列能量：其定义为

它相当于一个序列与自身共軛的相乘。

序列功率的计算与能量相仿。4、序列能量的计算26周期序列：对所有n,满足的称为周期序列。满足此式最小的N称为基本周期。给定x(n),生成周期序列

(xtide)的MATLAB方法：1.简单复制法：nxtide=[x,x,x,…]2.求余函数(mod)法：nxtide=0:K*N-1;xtide=x(mod(nxtide,N)+1)

5、周期序列及其产生27例4：设x=[2,3,4,5]，求将它延拓6个周期所得的序列。解：x=[2,3,4,5];N=length(x);K=6;ny=0:K*N－1；y=x(mod(ny,N)＋1)

5、周期序列及其产生286、序列的卷积运算序列x(n)和h(n)作卷积运算的定义此定义用图解表示时，先把x(k)和h(n-k)按k对齐。这时h(n-k)要解释为对k折叠并移位了n的h(k)。输出y(n)是把x(k)和h(n-k)覆盖区内k相同的样本乘积的求和而得。卷积输出的长度为两个输入长度之和减1：

Ly=Lx+Lh-129设x=[3,-3,7,0,-1,5,2];nx=[-4:2];h=[2,3,0,-5,2,1];nh=[-1:4];

求其卷积y(n)。解：用MATLAB函数conv可计算卷积

y=conv(x,h) 不过它没有给出y的位置向量n.6、序列的卷积运算30MATLAB提供

conv函数实现标准的一维信号卷积:

例如，若系统h(n)为 >>h=[111]输入序列x(n)为 >>x=[111]则x(n)经过系统h(n)后的MATLAB实现为： >>conv(h,x)或

conv([111],[111])执行后即得到y(n)为 ans= 12321注意：使用conv函数时，h(n)和x(n)都必须是有限长的,，否则不能使用conv函数。6、序列的卷积运算31自编带有位置向量的卷积子程序convwthn.mfunction[y,ny]=convwthn(x,nx,h,nh)nys=nx(1)+nh(1);nyf=nx(end)+nh(end);%MATLAB允许用end表示最后一个下标

y=conv(x,h);ny=[nys:nyf];[y,ny]=convwthn(x,nx,h,nh)得到：y=635619-313018-27-192ny=-5-4-3-2-10123456序列卷积用MATLAB运算算例32卷积举例：已知序列：，求卷积（只需显示卷积结果的序列图）。

33n1=-2:2;xn=0.5.*n1;n2=0:3;hn=2.^n2;%点乘k0=n1(1)+n2(1);%计算序列yn非零样值的起点位置k1=length(n1)+length(n2)-2;;%计算序列yn非零样值的宽度nyn=k0:k0+k1;%确定序列yn非零样值的序号向量yn=conv(xn,hn)stem(nyn,yn,'b','*')xlabel('时间n'),ylabel('幅度')title(‘y(n)’)gridon;卷积举例：34卷积举例：357、序列的相关

序列互相关的定义：它提供了两个序列之间相关程度的度量。它的定义与卷积相仿，只是不必将一个序列折叠。相关输出长度也是两个输入长度之和减一，lrxy=lx+ly-1自相关的定义：36MATLAB提供的计算序列相关性的函数xcorr.m

rxy=xcorr(x,y)

它也不给出rxy的位置向量。

7、序列的相关

37例5、令x(n)=[3,5,-7,2,-1,-3,2]为一原型序列（可看作雷达发射信号）；令y(n)为x(n)加入噪声干扰并移位后的序列（可看作雷达回波信号）: y(n)=x(n-2)+w(n)

计算两者的相关性，可以判断被噪声淹没的回波的延迟，从而确定目标的距离。7、序列的相关

387、序列的相关

397.1时域离散信号的产生及时域处理40三、滤波数字滤波器的系统函数H(z)用如下式表示:在MATLAB中，用向量b，a来表示滤波器的系数b(i)和

a(i)。

411、滤波器分类当n

=

0，m≠0时，称为AR滤波器，即自回归（AutoRecurrence）滤波器，具无限冲激响应（IIR），也即其单位采样响应h(n)具无限长度；若m

=

0，a(1)≠0，称为MA滤波器，即滑动平均（MovingAverage）滤波器，其单位采样响应h(n)是有限长度，故称有限冲激响应（FIR）滤波器；如果n、m都大于零，称为ARMA滤波器，而其冲激响应也为IIR。422、filter函数

MATLAB提供了filter函数来对离散信号进行滤波，表达信号通过系统后的结果。例如，一个单极点的低通滤波器系数如下：

>>b=1; %分子系数向量b(i) >>a=[1-0.9]; %分母系数向量a(i)

如果用filter函数实现对信号x滤波，只要调用： >>y=filter(b,a,x);

就可给出输入x经过滤波以后的输出y。433、单位冲激响应数字滤波器的单位冲激响应定义为输入为单位样本序列时数字滤波器的响应，即：h(n)=T[δ(n)] 其中：444、单位冲激响应的MATLAB实现MATLAB近似实现单位采样信号的方法为：

imp=[1;zeros(p,1)]; %zeros(p,1)产生p个零元素组成的列向量，p是正整数。使用imp后，滤波器的冲激响应可近似得到为：

h=filter(b,a,imp);stem(h);impz函数可以直接求出数字滤波器的单位冲激响应，即：impz(b,a)该命令将同时绘出滤波器的单位冲激响应（实验课进行验证）。45四、离散傅立叶变换(DFT)46定义：有限长序列x(n)(0≤n≤N-1)有N个样本值。它的傅立叶变换X(jω)在频率区间(0≤ω<2π)的N个等间隔分布的点ωk=2πk/N(0≤k≤N-1)上也有N个取样值。这两组同样长的序列之间可以用简单的关系联系起来：其中：1、DFT定义47它的反变换称为离散傅立叶反变换(IDFT)，可以表为：1、DFT定义48把DFT写成矩阵形式：记作：2、离散傅立叶变换的计算49傅立叶反变换(IDFT)的矩阵可表为：记作3、离散傅立叶反变换(IDFT)50若用矩阵式X=W*x及x=W-1*X求变换矩阵：W=dftmtx(N)，

反变换矩阵：W-1=conj(dftmtx(N))/N，但通常可直接调用FFT程序：正变换：X=fft(x,N) X=

fft(x)N默认为x的长度反变换：X=ifft(x)或X=

ifft(x,N)

4、用工具箱函数计算DFT51FFT程序求频谱的范围规范化为[0,2π)，对应的频谱下标为k=0,1,…N-1。设频率分辨率 dw=2π/N，则频率向量为

w=k*dw若要画出乃奎斯特频率范围[-π,π)内的频谱，需要把[π,2π)范围的频谱平移到[-π,0)内，对于这一运算，有专门的工具箱函数fftshift来实现。

X1=fftshift(X) X1就是在[-π,π)内的频谱向量。4、用工具箱函数计算DFT52傅里叶变换函数常用的FFT及反变换函数如下表所示。函数说明fft计算快速离散傅立叶变换fftshift调整fft函数的输出顺序，将零频位置移到频谱的中心abs求幅值angle求相角ifft计算离散傅立叶反变换5354FFT计算频谱举例例：选择合适的变换区间长度N，用DFT对下列信号进行谱分析，画出幅频、相频特性。解：(1)为周期序列，取一个周期用DFT谱分析%(1)周期序列,截取一个周期作谱分析N1=100;n1=0:N1-1;x1n=2*cos(0.2*pi*n1);%生成序列x1nX1k=fft(x1n,N1);%x1n的N1点DFTsubplot(1,2,1);stem(n1/N1*pi,abs(X1k),'.');title('X1k-Amplitude');subplot(1,2,2);stem(n1/N1*pi,angle(X1k),'.');title('X1k-Phrase');54例4.1求x=[2,-1,1,1]的离散傅立叶变换。并画到图上。

解：holdon Xd=fft([2,-1,1,1])%其DFT只有四点plot([0:3]*2*pi/4,abs(Xd),’o’)Xd1=fftshift(Xd);%画出移位后的FFT

plot([-2:1]*2*pi/4,abs(Xd1),’x’)4、用工具箱函数计算DFT55上述用频域DFT相乘再求反变换方法3可以用下列方框图来表示：这样求出的是循环卷积。可以证明，如果把N选得足够大，N≥N1＋N2-1，进行N点循环卷积，就有可能避免头尾混叠，结果就等于线性卷积。3.7用DFT计算线性卷积

56DFT应用举例—求卷积例:

计算快速卷积57DFT应用举例—求卷积%%fastconvolutionclear;closeall;clc;xn=sin(0.4*[1:1000]);hn=0.9.^[1:1000];M=length(xn);N=length(hn);nx=1:M;nh=1:N;L=pow2(nextpow2(M+N-1));%取L大于等于且最接近(M+N-1)的、2的正次幂tic,Xk=fft(xn,L);%L-pointFFT[x(n)]Hk=fft(hn,L);%L-pointFFT[x(n)]Yk=Xk.*Hk;%multiplyinfrequency:Y(k)yn=ifft(Yk,L);%L-pointIFFT;convolutionresultdisplay('(1)DFT');toc%time-out58DFT应用举例—求卷积subplot(221);plot(nx,xn,'.');

ylabel('x(n)');subplot(222);plot(nh,hn,'.');ylabel('h(n)');subplot(2,2,[34]);ny=1:L;stem(ny,real(yn),'.');ylabel('y(n)');%--------------------------------------------%tic,yn=conv(xn,hn);display('(2)function''conv''');toc59例1：已知一个8点