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之后

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

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

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

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

3)腳本7

腳本文件可以理解為簡單的M文件，腳本文件中的變量都是全局變量。函數文件是在腳本文件的基礎之上多添加了函數定義行，其代碼組織結構和調用方式與對應的腳本文件截然不同。函數文件是以函數聲明行“function”作為開始，如果沒有特別聲明，函數文件中的變量都是局部變量。函數運行完畢之后，其定義的變量將從工作區(qū)間中清除。而腳本文件只是將一系列相關的代碼結合封裝，沒有輸入參數和輸出參數，即不自帶參數，也不一定要返回結果。而多數函數文件一般都有輸入和輸出變量，并返回結果。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)覆蓋區(qū)內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函數實現(xiàn)標準的一維信號卷積:

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

conv([111],[111])執(zhí)行后即得到y(tǒng)(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]為一原型序列（可看作雷達發(fā)射信號）；令y(n)為x(n)加入噪聲干擾并移位后的序列（可看作雷達回波信號）: y(n)=x(n-2)+w(n)

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

387、序列的相關

397.1時域離散信號的產生及時域處理40三、濾波數字濾波器的系統(tǒng)函數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函數來對離散信號進行濾波，表達信號通過系統(tǒng)后的結果。例如，一個單極點的低通濾波器系數如下：

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

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

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

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ω)在頻率區(qū)間(0≤ω<2π)的N個等間隔分布的點ωk=2πk/N(0≤k≤N-1)上也有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，但通常可直接調用FFT程序：正變換：X=fft(x,N) X=

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

ifft(x,N)

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

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

X1=fftshift(X) X1就是在[-π,π)內的頻譜向量。4、用工具箱函數計算DFT52傅里葉變換函數常用的FFT及反變換函數如下表所示。函數說明fft計算快速離散傅立葉變換fftshift調整fft函數的輸出順序，將零頻位置移到頻譜的中心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、用工具箱函數計算DFT55上述用頻域DFT相乘再求反變換方法3可以用下列方框圖來表示：這樣求出的是循環(huán)卷積?？梢宰C明，如果把N選得足夠大，N≥N1＋N2-1，進行N點循環(huá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點