LDX第5次課-信號處理算法Matlab仿真

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

M文件編輯器3

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

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

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

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

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

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

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

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

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

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

3)腳本7

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

以采樣周期T進行等間隔采樣,

得到采樣序列x(n)它的每一個樣本點與相應(yīng)處的模擬信號值相等，而在其它時間位置處的值為零。10引言：模擬信號與其采樣序列11通常，用式（1）的序列作為離散信號的研究模型，由于抽掉了模擬信號中的主要物理量——時間。在數(shù)字域處理信號時,通常都是按序號逐次進行運算和存儲，這就是抽象為序列模型的好處。引言：模擬信號與其采樣序列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是順序增加的整數(shù)序列，可寫成n=ns:nf,序列的長度為Lx=nf-ns+1。本例n=-3:4;Lx=8。13有限長序列和無限長序列序列位置向量的起點和終點都是有限數(shù)稱為有限長序列序列位置向量的起點向-∞

延伸為無限長序列，簡稱左序列。序列位置向量的終點向∞延伸也為無限長序列，簡稱右序列。工程中遇到的序列都屬于有限序列，有些情況下可看成右序列，但不會有左序列。一、常用的序列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ù)不滿足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)復(fù)制191、序列相乘相加：將兩個序列x1,x2的位置向量n1與n2對準(zhǔn)，長度取兩者之‘合’；形成一個新的包含n1與n2的公共位置向量n。其起點為兩者中小的起點，終點為兩者中大的終點，

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

將x1,x2補零，在n上延長，構(gòu)成y1,y2;各個對應(yīng)元素相乘或相加；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取負(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點為中心，左右翻轉(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=fliplr(x);

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

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

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

(xtide)的MATLAB方法：1.簡單復(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個周期所得的序列。解：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、序列的卷積運算序列x(n)和h(n)作卷積運算的定義此定義用圖解表示時，先把x(k)和h(n-k)按k對齊。這時h(n-k)要解釋為對k折疊并移位了n的h(k)。輸出y(n)是把x(k)和h(n-k)覆蓋區(qū)內(nèi)k相同的樣本乘積的求和而得。卷積輸出的長度為兩個輸入長度之和減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可計算卷積

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

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

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

conv([111],[111])執(zhí)行后即得到y(tǒng)(n)為 ans= 12321注意：使用conv函數(shù)時，h(n)和x(n)都必須是有限長的,，否則不能使用conv函數(shù)。6、序列的卷積運算31自編帶有位置向量的卷積子程序convwthn.mfunction[y,ny]=convwthn(x,nx,h,nh)nys=nx(1)+nh(1);nyf=nx(end)+nh(end);%MATLAB允許用end表示最后一個下標(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運算算例32卷積舉例：已知序列：，求卷積（只需顯示卷積結(jié)果的序列圖）。

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、序列的相關(guān)

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

rxy=xcorr(x,y)

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

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

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

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

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

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

a(i)。

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

=

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

=

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

MATLAB提供了filter函數(shù)來對離散信號進行濾波，表達信號通過系統(tǒng)后的結(jié)果。例如，一個單極點的低通濾波器系數(shù)如下：

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

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

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

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

h=filter(b,a,imp);stem(h);impz函數(shù)可以直接求出數(shù)字濾波器的單位沖激響應(yīng)，即：impz(b,a)該命令將同時繪出濾波器的單位沖激響應(yīng)（實驗課進行驗證）。45四、離散傅立葉變換(DFT)46定義：有限長序列x(n)(0≤n≤N-1)有N個樣本值。它的傅立葉變換X(jω)在頻率區(qū)間(0≤ω<2π)的N個等間隔分布的點ωk=2πk/N(0≤k≤N-1)上也有N個取樣值。這兩組同樣長的序列之間可以用簡單的關(guān)系聯(liá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，但通?？芍苯诱{(diào)用FFT程序：正變換：X=fft(x,N) X=

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

ifft(x,N)

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

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

X1=fftshift(X) X1就是在[-π,π)內(nèi)的頻譜向量。4、用工具箱函數(shù)計算DFT52傅里葉變換函數(shù)常用的FFT及反變換函數(shù)如下表所示。函數(shù)說明fft計算快速離散傅立葉變換fftshift調(diào)整fft函數(shù)的輸出順序，將零頻位置移到頻譜的中心abs求幅值angle求相角ifft計算離散傅立葉反變換5354FFT計算頻譜舉例例：選擇合適的變換區(qū)間長度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、用工具箱函數(shù)計算DFT55上述用頻域DFT相乘再求反變換方法3可以用下列方框圖來表示：這樣求出的是循環(huán)卷積?？梢宰C明，如果把N選得足夠大，N≥N1＋N2-1，進行N點循環(huán)卷積，就有可能避免頭尾混疊，結(jié)果就等于線性卷積。3.7用DFT計算線性卷積

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

計算快速卷積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：已知一個8點