LDX第5次課-信號(hào)處理算法Matlab仿真

第5講MATLAB數(shù)字信號(hào)處理應(yīng)用西安郵電大學(xué)通信與信息工程學(xué)院數(shù)字信號(hào)處理與Matlb仿真MATLAB桌面復(fù)習(xí)：Matlab編程2M文件介紹1)M文件編輯器M文件編輯器提供了一個(gè)進(jìn)行文本編輯和M文件調(diào)試的圖形用戶(hù)界面，如圖所示。

M文件編輯器3

新建一個(gè)M文件有以下三種途徑：①菜單操作（FileNewScript）②命令操作（editM文件名）③命令按鈕（快捷鍵）

打開(kāi)已有的M文件有四種途徑：①菜單操作（FileOpen）②命令操作（editM文件名）③命令按鈕（快捷鍵）④雙擊M文件名M文件介紹4圖所示為M文件編輯器新建/打開(kāi)快捷鍵。

M文件編輯器新建/打開(kāi)快捷鍵52)函數(shù)函數(shù)文件由function語(yǔ)句引導(dǎo)，5個(gè)部分。

函數(shù)題頭：指函數(shù)的定義行，是函數(shù)語(yǔ)句的第一行，在該行中定義函數(shù)名、輸入變量列表及輸出變量列表等。

HI行：指函數(shù)幫助文本的第一行，為該函數(shù)文件的幫助主題，當(dāng)使用lookfor命令時(shí)，可以查看到該行信息。

幫助信息：提供了函數(shù)的完整的幫助信息，包括HI之后

至第一個(gè)可執(zhí)行的行或空行為止的所有注釋語(yǔ)句。

函數(shù)體：指函數(shù)代碼段，也是函數(shù)的主體部分。

注釋部分：指對(duì)函數(shù)體中各語(yǔ)句的解釋和說(shuō)明文本，注

釋語(yǔ)句是以符號(hào)%引導(dǎo)的。6腳本是若干命令或函數(shù)的集合，用于執(zhí)行特定的功能。在MATLAB命令窗口直接輸入此文件的主文件名，MATLAB可順序執(zhí)行此文件內(nèi)的所有命令，不需要在其中輸入?yún)?shù)，也不需要給出輸出變量來(lái)接收處理結(jié)果。腳本的操作對(duì)象為MATLAB工作空間內(nèi)的變量，腳本式M文件運(yùn)行產(chǎn)生的所有變量都是全局變量，這些變量都駐留在MATLAB基本工作空間內(nèi)，通過(guò)使用clear命令或者關(guān)閉MATLAB系統(tǒng)可以清除這些變量。

3)腳本7

腳本文件可以理解為簡(jiǎn)單的M文件，腳本文件中的變量都是全局變量。函數(shù)文件是在腳本文件的基礎(chǔ)之上多添加了函數(shù)定義行，其代碼組織結(jié)構(gòu)和調(diào)用方式與對(duì)應(yīng)的腳本文件截然不同。函數(shù)文件是以函數(shù)聲明行“function”作為開(kāi)始，如果沒(méi)有特別聲明，函數(shù)文件中的變量都是局部變量。函數(shù)運(yùn)行完畢之后，其定義的變量將從工作區(qū)間中清除。而腳本文件只是將一系列相關(guān)的代碼結(jié)合封裝，沒(méi)有輸入?yún)?shù)和輸出參數(shù)，即不自帶參數(shù)，也不一定要返回結(jié)果。而多數(shù)函數(shù)文件一般都有輸入和輸出變量，并返回結(jié)果。8一、常用序列主要內(nèi)容引言二、序列的基本運(yùn)算三、濾波四、離散付立葉變換（DFT）及應(yīng)用引言：模擬信號(hào)與其采樣序列對(duì)模擬信號(hào)

以采樣周期T進(jìn)行等間隔采樣,

得到采樣序列x(n)它的每一個(gè)樣本點(diǎn)與相應(yīng)處的模擬信號(hào)值相等，而在其它時(shí)間位置處的值為零。10引言：模擬信號(hào)與其采樣序列11通常，用式（1）的序列作為離散信號(hào)的研究模型，由于抽掉了模擬信號(hào)中的主要物理量——時(shí)間。在數(shù)字域處理信號(hào)時(shí),通常都是按序號(hào)逐次進(jìn)行運(yùn)算和存儲(chǔ)，這就是抽象為序列模型的好處。引言：模擬信號(hào)與其采樣序列12一、常用的序列

序列的表示方法完全地表示序列x(n)要用x和n兩個(gè)向量，x表示序列的值，n表示它的位置或時(shí)間次序。例如x(n)={2,1.2,-1.4,3,1,4,3.1,7}（下劃雙線(xiàn)表示n=0處的采樣點(diǎn)）用MATLAB表示時(shí)，用兩個(gè)向量：

n=[-3,-2,-1,0,1,2,3,4];x=[2,1.2,-1.4,3,1,4,3.1,7];可見(jiàn)n是順序增加的整數(shù)序列，可寫(xiě)成n=ns:nf,序列的長(zhǎng)度為L(zhǎng)x=nf-ns+1。本例n=-3:4;Lx=8。13有限長(zhǎng)序列和無(wú)限長(zhǎng)序列序列位置向量的起點(diǎn)和終點(diǎn)都是有限數(shù)稱(chēng)為有限長(zhǎng)序列序列位置向量的起點(diǎn)向-∞

延伸為無(wú)限長(zhǎng)序列，簡(jiǎn)稱(chēng)左序列。序列位置向量的終點(diǎn)向∞延伸也為無(wú)限長(zhǎng)序列，簡(jiǎn)稱(chēng)右序列。工程中遇到的序列都屬于有限序列，有些情況下可看成右序列，但不會(huì)有左序列。一、常用的序列14一、常用的序列1．單位取樣序列2．單位階躍序列15矩形序列復(fù)數(shù)指數(shù)序列正余弦序列

x(n)=cos(ωn)+jsin(ωn)一、常用的序列16一、常用的序列17function[x,n]=impseq(np,ns,nf)%ns、nf分別為序列起始與終止位置,np為脈沖位置ifns>np|ns>nf|np>nf error（’輸入位置參數(shù)不滿(mǎn)足ns<=np<=nf’）else n=[ns:nf];x=[(n-np)==0];%核心語(yǔ)句end單位階躍序列：把核心語(yǔ)句改為：x=[(n-np)>=0]例1：用Matlab生成脈沖序列18二、序列的基本運(yùn)算序列的基本運(yùn)算(a)序列相乘

(b)序列相加

(c) 倍乘(d)右移位(遲延)(e)左移位(提前)(f)復(fù)制191、序列相乘相加：將兩個(gè)序列x1,x2的位置向量n1與n2對(duì)準(zhǔn)，長(zhǎng)度取兩者之‘合’；形成一個(gè)新的包含n1與n2的公共位置向量n。其起點(diǎn)為兩者中小的起點(diǎn)，終點(diǎn)為兩者中大的終點(diǎn)，

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

將x1,x2補(bǔ)零，在n上延長(zhǎng)，構(gòu)成y1,y2;各個(gè)對(duì)應(yīng)元素相乘或相加；y=y1.*y2;(或y=y1+y2;)二、序列的基本運(yùn)算20%y(n)的長(zhǎng)度n=min(min(n1),min(n2)):max(max(n1),max(n2));%y1,y2初始化y1=zeros(1,length(n));y2=y1;%把x1補(bǔ)零，擴(kuò)展為具有y的長(zhǎng)度的y1y1(find((n>=min(n1))&(n<=max(n1))))=x1;%把x2補(bǔ)零，擴(kuò)展為具有y的長(zhǎng)度的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的計(jì)算圖解23序列向量的值不變，只把位置向量右移k拍子程序[y,ny]=seqshift(x,nx,k)的核心語(yǔ)句：y=x;ny=nx+k左移時(shí)k取負(fù)值例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點(diǎn)為中心，左右翻轉(zhuǎn)。 例如：

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

則 y=[4,-2,7,3];ny=-2:1即x的值及其位置排序都要左右翻轉(zhuǎn)子程序[y,ny]=seqfold(x,nx)的核心語(yǔ)句：

y=fliplr(x);

ny=-fliplr(nx)3、序列的折疊：y(n)=x(-n)25序列能量：其定義為

它相當(dāng)于一個(gè)序列與自身共軛的相乘。

序列功率的計(jì)算與能量相仿。4、序列能量的計(jì)算26周期序列：對(duì)所有n,滿(mǎn)足的稱(chēng)為周期序列。滿(mǎn)足此式最小的N稱(chēng)為基本周期。給定x(n),生成周期序列

(xtide)的MATLAB方法：1.簡(jiǎn)單復(fù)制法：nxtide=[x,x,x,…]2.求余函數(shù)(mod)法：nxtide=0:K*N-1;xtide=x(mod(nxtide,N)+1)

5、周期序列及其產(chǎn)生27例4：設(shè)x=[2,3,4,5]，求將它延拓6個(gè)周期所得的序列。解：x=[2,3,4,5];N=length(x);K=6;ny=0:K*N－1；y=x(mod(ny,N)＋1)

5、周期序列及其產(chǎn)生286、序列的卷積運(yùn)算序列x(n)和h(n)作卷積運(yùn)算的定義此定義用圖解表示時(shí)，先把x(k)和h(n-k)按k對(duì)齊。這時(shí)h(n-k)要解釋為對(duì)k折疊并移位了n的h(k)。輸出y(n)是把x(k)和h(n-k)覆蓋區(qū)內(nèi)k相同的樣本乘積的求和而得。卷積輸出的長(zhǎng)度為兩個(gè)輸入長(zhǎng)度之和減1：

Ly=Lx+Lh-129設(shè)x=[3,-3,7,0,-1,5,2];nx=[-4:2];h=[2,3,0,-5,2,1];nh=[-1:4];

求其卷積y(n)。解：用MATLAB函數(shù)conv可計(jì)算卷積

y=conv(x,h) 不過(guò)它沒(méi)有給出y的位置向量n.6、序列的卷積運(yùn)算30MATLAB提供

conv函數(shù)實(shí)現(xiàn)標(biāo)準(zhǔn)的一維信號(hào)卷積:

例如，若系統(tǒng)h(n)為 >>h=[111]輸入序列x(n)為 >>x=[111]則x(n)經(jīng)過(guò)系統(tǒng)h(n)后的MATLAB實(shí)現(xiàn)為： >>conv(h,x)或

conv([111],[111])執(zhí)行后即得到y(tǒng)(n)為 ans= 12321注意：使用conv函數(shù)時(shí)，h(n)和x(n)都必須是有限長(zhǎng)的,，否則不能使用conv函數(shù)。6、序列的卷積運(yùn)算31自編帶有位置向量的卷積子程序convwthn.mfunction[y,ny]=convwthn(x,nx,h,nh)nys=nx(1)+nh(1);nyf=nx(end)+nh(end);%MATLAB允許用end表示最后一個(gè)下標(biāo)

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運(yùn)算算例32卷積舉例：已知序列：，求卷積（只需顯示卷積結(jié)果的序列圖）。

33n1=-2:2;xn=0.5.*n1;n2=0:3;hn=2.^n2;%點(diǎn)乘k0=n1(1)+n2(1);%計(jì)算序列yn非零樣值的起點(diǎn)位置k1=length(n1)+length(n2)-2;;%計(jì)算序列yn非零樣值的寬度nyn=k0:k0+k1;%確定序列yn非零樣值的序號(hào)向量yn=conv(xn,hn)stem(nyn,yn,'b','*')xlabel('時(shí)間n'),ylabel('幅度')title(‘y(n)’)gridon;卷積舉例：34卷積舉例：357、序列的相關(guān)

序列互相關(guān)的定義：它提供了兩個(gè)序列之間相關(guān)程度的度量。它的定義與卷積相仿，只是不必將一個(gè)序列折疊。相關(guān)輸出長(zhǎng)度也是兩個(gè)輸入長(zhǎng)度之和減一，lrxy=lx+ly-1自相關(guān)的定義：36MATLAB提供的計(jì)算序列相關(guān)性的函數(shù)xcorr.m

rxy=xcorr(x,y)

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

7、序列的相關(guān)

37例5、令x(n)=[3,5,-7,2,-1,-3,2]為一原型序列（可看作雷達(dá)發(fā)射信號(hào)）；令y(n)為x(n)加入噪聲干擾并移位后的序列（可看作雷達(dá)回波信號(hào)）: y(n)=x(n-2)+w(n)

計(jì)算兩者的相關(guān)性，可以判斷被噪聲淹沒(méi)的回波的延遲，從而確定目標(biāo)的距離。7、序列的相關(guān)

387、序列的相關(guān)

397.1時(shí)域離散信號(hào)的產(chǎn)生及時(shí)域處理40三、濾波數(shù)字濾波器的系統(tǒng)函數(shù)H(z)用如下式表示:在MATLAB中，用向量b，a來(lái)表示濾波器的系數(shù)b(i)和

a(i)。

411、濾波器分類(lèi)當(dāng)n

=

0，m≠0時(shí)，稱(chēng)為AR濾波器，即自回歸（AutoRecurrence）濾波器，具無(wú)限沖激響應(yīng)（IIR），也即其單位采樣響應(yīng)h(n)具無(wú)限長(zhǎng)度；若m

=

0，a(1)≠0，稱(chēng)為MA濾波器，即滑動(dòng)平均（MovingAverage）濾波器，其單位采樣響應(yīng)h(n)是有限長(zhǎng)度，故稱(chēng)有限沖激響應(yīng)（FIR）濾波器；如果n、m都大于零，稱(chēng)為ARMA濾波器，而其沖激響應(yīng)也為IIR。422、filter函數(shù)

MATLAB提供了filter函數(shù)來(lái)對(duì)離散信號(hào)進(jìn)行濾波，表達(dá)信號(hào)通過(guò)系統(tǒng)后的結(jié)果。例如，一個(gè)單極點(diǎn)的低通濾波器系數(shù)如下：

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

如果用filter函數(shù)實(shí)現(xiàn)對(duì)信號(hào)x濾波，只要調(diào)用： >>y=filter(b,a,x);

就可給出輸入x經(jīng)過(guò)濾波以后的輸出y。433、單位沖激響應(yīng)數(shù)字濾波器的單位沖激響應(yīng)定義為輸入為單位樣本序列時(shí)數(shù)字濾波器的響應(yīng)，即：h(n)=T[δ(n)] 其中：444、單位沖激響應(yīng)的MATLAB實(shí)現(xiàn)MATLAB近似實(shí)現(xiàn)單位采樣信號(hào)的方法為：

imp=[1;zeros(p,1)]; %zeros(p,1)產(chǎn)生p個(gè)零元素組成的列向量，p是正整數(shù)。使用imp后，濾波器的沖激響應(yīng)可近似得到為：

h=filter(b,a,imp);stem(h);impz函數(shù)可以直接求出數(shù)字濾波器的單位沖激響應(yīng)，即：impz(b,a)該命令將同時(shí)繪出濾波器的單位沖激響應(yīng)（實(shí)驗(yàn)課進(jìn)行驗(yàn)證）。45四、離散傅立葉變換(DFT)46定義：有限長(zhǎng)序列x(n)(0≤n≤N-1)有N個(gè)樣本值。它的傅立葉變換X(jω)在頻率區(qū)間(0≤ω<2π)的N個(gè)等間隔分布的點(diǎn)ωk=2πk/N(0≤k≤N-1)上也有N個(gè)取樣值。這兩組同樣長(zhǎng)的序列之間可以用簡(jiǎn)單的關(guān)系聯(lián)系起來(lái)：其中：1、DFT定義47它的反變換稱(chēng)為離散傅立葉反變換(IDFT)，可以表為：1、DFT定義48把DFT寫(xiě)成矩陣形式：記作：2、離散傅立葉變換的計(jì)算49傅立葉反變換(IDFT)的矩陣可表為：記作3、離散傅立葉反變換(IDFT)50若用矩陣式X=W*x及x=W-1*X求變換矩陣：W=dftmtx(N)，

反變換矩陣：W-1=conj(dftmtx(N))/N，但通?？芍苯诱{(diào)用FFT程序：正變換：X=fft(x,N) X=

fft(x)N默認(rèn)為x的長(zhǎng)度反變換：X=ifft(x)或X=

ifft(x,N)

4、用工具箱函數(shù)計(jì)算DFT51FFT程序求頻譜的范圍規(guī)范化為[0,2π)，對(duì)應(yīng)的頻譜下標(biāo)為k=0,1,…N-1。設(shè)頻率分辨率 dw=2π/N，則頻率向量為

w=k*dw若要畫(huà)出乃奎斯特頻率范圍[-π,π)內(nèi)的頻譜，需要把[π,2π)范圍的頻譜平移到[-π,0)內(nèi)，對(duì)于這一運(yùn)算，有專(zhuān)門(mén)的工具箱函數(shù)fftshift來(lái)實(shí)現(xiàn)。

X1=fftshift(X) X1就是在[-π,π)內(nèi)的頻譜向量。4、用工具箱函數(shù)計(jì)算DFT52傅里葉變換函數(shù)常用的FFT及反變換函數(shù)如下表所示。函數(shù)說(shuō)明fft計(jì)算快速離散傅立葉變換fftshift調(diào)整fft函數(shù)的輸出順序，將零頻位置移到頻譜的中心abs求幅值angle求相角ifft計(jì)算離散傅立葉反變換5354FFT計(jì)算頻譜舉例例：選擇合適的變換區(qū)間長(zhǎng)度N，用DFT對(duì)下列信號(hào)進(jìn)行譜分析，畫(huà)出幅頻、相頻特性。解：(1)為周期序列，取一個(gè)周期用DFT譜分析%(1)周期序列,截取一個(gè)周期作譜分析N1=100;n1=0:N1-1;x1n=2*cos(0.2*pi*n1);%生成序列x1nX1k=fft(x1n,N1);%x1n的N1點(diǎn)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]的離散傅立葉變換。并畫(huà)到圖上。

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

plot([-2:1]*2*pi/4,abs(Xd1),’x’)4、用工具箱函數(shù)計(jì)算DFT55上述用頻域DFT相乘再求反變換方法3可以用下列方框圖來(lái)表示：這樣求出的是循環(huán)卷積?？梢宰C明，如果把N選得足夠大，N≥N1＋N2-1，進(jìn)行N點(diǎn)循環(huán)卷積，就有可能避免頭尾混疊，結(jié)果就等于線(xiàn)性卷積。3.7用DFT計(jì)算線(xiàn)性卷積

56DFT應(yīng)用舉例—求卷積例:

計(jì)算快速卷積57DFT應(yīng)用舉例—求卷積%%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應(yīng)用舉例—求卷積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：已知一個(gè)8點(diǎn)