數字信號處理

實驗指導書

實驗一時域離散信號的產生

一、實驗目的

1、了解常用時域離散信號及其特點；

2、掌握MATLAB程序的編程方法；

3、熟悉MATLAB函數的調用方法。

二、實驗原理

在時間軸上的離散點取值的信號，稱為離散時間信號。離散時間信號只在某些離散的

瞬時給出函數值，而在其他時刻無定義。它是時間上不連續按一定先后次序排列的一組數的

集合，稱為時間序列，用x(n)表示，n取整數代表時間的離散時刻。

在MATLAB中用向量來表示一個有限長度的序列。

常用離散信號：

1、單位抽樣序列

77=0

3(〃)='

nw0

2、單位階躍序列

/?>()

〃<0

3、實指數序列

x(〃)=aH

4、復指數序列

5、正(余)弦序列

x(n)=Umsin(d)on+0)

6、隨機序列

在利用計算機進行系統的研究時，經常需要產生隨機信號，MATLAB提供?個工具

函數rand來產生隨機信號。

7、周期序列

x(n)=x(n+N)

三、實驗用函數

1、stem

功能：繪制二維圖形。

調用格式：

stem(n,x)；n為橫軸，x為縱軸的線性圖形。

2、length

功能：計算某一變量的長度或采樣點數。

調用格式：

N=length(t)；”,算時間向量t的個數并賦給變量N。

3、axis

功能：限定圖形坐標的范圍。

調用格式：

axis([xl,x2,yl,y2])；橫坐標從xl—x2,縱坐標從yl—y2。

4、zeros

功能：產生一個全0序列。

調用格式：

x=zeros(l,n)；產生n個0的序列。

5、ones

功能：產生一個全1序列。

調用格式：

y=ones(l,n)；產生n個1的序列。

四、參考實例

例1.1用Matlab產生單位抽樣序列。

%先建立函數impseq(nl,n2,n0)

function[x,n]=impseq(n1,n2,n0)

n=[nl:n2];

x=[(n-nO)==OJ;

%編寫主程序調用該函數

[x,n]-impseq(-2,8,2);

2

stem(n,x)

程序運行結果如圖1-1所示:

圖1-1單位抽樣序列

例1.2實數指數序列(運算符’》”)

Matlab程序如下:

n=[0:10];

x=0.9.An;

stem(n,x)

程序運行結果如圖1-2所示

圖1-2實數指數序列

3

例1.3復數指數序列(x(n)=4-°」+加3)“(_]0</7<10))

Matlab程序如下：

n=[-IO:lO];alpha=-0.1+0.3*j;x=exp(alpha*n);

rcal_x=rcal(x);imagc_x=imag(x);

mag_x=abs(x);phase_x=angle(x);

subplot(2,2,1),stem(n,real_x)

subplot(2,2,2);stem(n,image_x)

subplot(2,2,3)；stem(n,mag_x)

subplot(2,2,4);stem(n,phase_x)

程序運行結果如圖1-3所示

圖1-3復數指數序列

例1.4正、余弦序列(x(n)=Umsin(4〃+0))

Matlab程序如下：

n=[0:10];

x=3*cos(0.1*pi*n+pi/3);

stem(n,x)

程序運行結果如圖1-4所示

4

igure1BB8

圖1-4正、余弦序列

例1.5隨機序列

rand(l,N)產生其元素在［0,I］之間均勻分布長度為N的隨機序列

randn(l,N)產生均值為(J,方差為1,長度為N的高斯隨機序列

例L6周期序列

如何生成周期序列

1、將一個周期復制p次；

2、借助矩陣運算、matlab下標能力。先生成一個包含p列x(n)值的矩陣，然后用結構(:)

來把p列串接成?個長周期序列。因為這個結構只能用于列向，最后還需要做矩陣轉置獲得

所需序列。

Matlab程序如下:

x=[l,2,3];%一個x(n)

xn=x'*ones(1.3)%生成p列x(n)

xn=xn(:)'%將p列串接成長列序列并轉置

steni(xn)

程序運行的結果如圖1-5所示

5

圖1-5周期序歹J

五、實驗任務

1、調試部分例題程序，掌握Matlab基本操作方法。

2、編寫程序，完成下列函數波形：

1）利用zeros函數生成單位抽樣序列；

2）利用zeros函數和ones函數生成單位階躍序列;

六、實驗報告

1、簡述實驗目的、原理。

2、寫出上機調試通過的實驗任務的程序并描述其圖形曲線。

6

實驗二離散序列的基本運算

一、實驗目的

1、加強MATLAB運用。

2、了解離散時間序列在時域中的基本運算。

3、熟悉相關函數的使用方法，掌握離散序列運算程序的編寫方法。

二、實驗原理

離散序列的時域運算包括信號的相加、相乘，信號的時域變換包括信號的移位、反折、

倒相及尺度變換等。

在MATLAB中，序列的相加和相乘運算是兩個向量之間的運算，因此參加運算的兩個

序列必須具有相同的長度，否則不能直接進行運算，需要進行相應的處理后再進行運算。

三、實驗用函數

I、find

功能：尋找非零元素的索引號。

調用格式：

find((n>=min(n1))&(n<=max(n1)))：在符號關系運算條件的范圍內尋找非零元素的索引

號。

2、fliplr

功能：對矩陣行元素進行左右翻轉。

調用格式：

xl=fliplr(x)：將x的行元素進行左右翻轉，賦給變量xl。

四、實例

1、信號的時域變換

1)序列移位

將一個離散序列進行移位，形成新的序列：xl(n)=xin-m)o當m>0時，原序列向右移m

位，當m<0時，原序列向左移。

%建立移位函數(sigshifl(x,m,nO))

function[y,n]=sigshift(x,m,nO)

n=m+nO;

y二x；

7

2)序列反折

在這個運算中，x(n)以n=0為基準點，以縱軸為對稱軸反折得到一個新的序歹人

y(n)=|x(-n)|

在MATLAB中提供了fliplr函數實現序列反折。

%建立反折函數(sigfold(x,n))

function[y,n]-sigfbld(x,n)

y=fliplr(x);

n=-fliplr(n);

3)序列倒相

是求一個與原序列的向量值相反，對應的時間向量不變的新序列。

4)序列的尺度變換

通過對時間軸的放大或壓縮形成新的序列。

2、序列的算術運算

1)序列相加

序列相加是指兩個序列中相同序號的序列值逐項對應相加，形成新的序列。

參加運算的兩個序列的維數不同時

%建立通用函數

function[y,n]=sigadd(xLn1,x2,n2)

n=min(min(n1),min(n2)):max(max(n1),max(n2));

yl=zcros(l,lcngth(n));

y2=yl;

y1(find((n>=min(nl))&(n<=max(nl))==l))=x1;

y2(find((n>=min(n2))&(n<=max(n2))==l))=x2;

y=yl+y2;

1)序列相乘

序列相加是指兩個序列中相同序號的序列值逐項對應相乘，形成新的序列。

參加運算的兩個序列的維數不同時處理方法與序列相加相同。

8

五、實驗任務

1、理解序列運算的性質，r解函數語句的意義。

2、利用例題函數完成下列序列運算

1)已知X](n)=u(n+1)(-3<n<5);

X2(n)=u(n-3)(-4<n<7)

求：x(n)-xi(n)+X2(n)

2)己知xi(n尸3e-0.25n(-2<n<8)

X2(n)=u(n+1)(-3<n<6)

求：x(n)=xi(n)dX2(n)

六、實驗報告

1、簡述實驗目的和原理。

2、列寫上機調試通過的程序，并描繪其波形曲線。

9

實驗三離散卷積的原理及應用

一、實驗目的

1、通過實驗進一步理解卷積定理，了解卷積過程；

2、掌握應用線性卷積求解離散時間系統響應的基本方法。

二、實驗原理

對于線性移不變離散系統,任意的輸入信號x(n)可以用以〃)及其位移的線性組合來表

示，即

工(〃)=￡x(k)6(n-k)

當輸入為5(〃)時，系4的輸出y(n)=h(n),由系統的線性移不變性質可以得到系統對x(n)

的響應y(n)為

y(n)=Wx(k)h(n-k)

k=f

稱為離散系統的線性卷積，簡寫為

y(n)=x(n)*h(n)

也就是說，如果已知系統的沖激響應，將輸入信號與系統的沖激響應進行卷積運算，

即可求得系統的響應。

三、實驗用函數

1、卷積函數conv

功能：進行兩個序列的卷積運算。

調用格式：

y=conv(x,h)；用于求解兩有限長序列的卷積。

2、sum

功能：求各元素之和。

調用格式：

y=sum(x)；求序列xt＞各元素之和。

3、hold

功能：控制當前圖形窗口是否刷新的雙向切換開關。

調用格式：

holdon：使當前圖形窗口中的圖形保持旦不被刷新，準備接受繪制新的圖形。

holdoff：使當前圖形窗口中的圖形不具備不被刷新的性質「

10

4、impz

功能：求解數字系統的沖激響應。

調用格式：

[h,t]=impz(b,a)；求解數字系統的沖激響應h,取樣點數為缺省值。

lh,t]=impz(b,a,n)：求解數字系統的沖激響應h,取樣點數由n確定。

impz(b,a);在當前窗口用slem(t,h)函數繪制圖形。

5、dslep

功能：求解數字系統的階躍響應。

調用格式：

[h,t]=dstep(b,a)；求解數字系統的沖激響應h,取樣點數為缺省值。

[h,t]=dstep(b,a,n)；求解數字系統的沖激響應h,取樣點數由n確定。

dstep(b,a)；在當前窗口用stairs(t,h)函數繪制圖形。

四、參考實例

在利用Madab提供的卷積函數進行卷積運算時，主要是確定卷積結果的時間區間。conv

函數默認兩信號的時間序列從n=0開始，卷積結果對應的時間序列也從n=0開始。圻果信

號不是從0開始，則編程時必須用兩個數組確定一個信號，其中，一個數組是信號波形的幅

度樣值，另一個數組是其對應的時間向量。

%建立一個適用于信號從任意時間開始的通用函數

function[y,ny]=sconv(x,h,nx.nh,p)

y=conv(x,h):

nl=nx(l)+nh(l);%計算y的非零樣值的起點位置

n2=nx(length(x))+nh(length(h));%計算y的非零樣值的寬度

ny=nl:p:n2;%確定y的非零樣值的時間向量

五、實驗任務

已知一個IIR數字低通濾波器的系統函數公式為

0.1321+0.3963ZT+0.3963Z-+0.1321

21-0.34319z-,+0.604392-2-0.204072-3

輸入一個矩形信號序列x=square(n/5)(-2<n<10^),求該系統的響應。

II

六、實驗報告

1、簡述實驗目的、原理。

2、寫出上機調試通過的實驗任務的程序并描述其圖形曲線。

12

實驗四離散系統變換域分析一Z變換

一、實驗目的

1、通過實驗進一步加深對z變換的理解；

2、掌握進行Z變換的基本方法，學會利用部分分式法在Z變換中的應用。

3、掌握MATLAB提供的子函數。

二、實驗原理

在離散時間信號與系統中，變換域分析法是Z變換法和傅立葉變換法。Z變換在離散時

間系統中的作用就如同拉普拉斯變換在連續時間系統中的作用i樣，它把描述離散系統的差

分方程轉化為簡單的代數方程，使其求解大大簡化。

1、序列x(n)的Z變換定義如下：

X(z)=Z[x(n)]=之x(n)z~n

n="XJ

z變換存在的充要條件是上面的級數收斂。

Z變換的收斂域：

使序列x(n)的z變換X(z)收斂的所有z值的集合稱作X(z)的收斂域。X(z)收斂的充要條

件是絕對可和。

ZIx{n}z~n|=A/<oo

〃=-00

幾類序列的收斂域.：

有限長序列：

0<|z|<oo

右邊序列：

RI<|z|<oo

左邊序列：

0<|z|<R2

雙邊序列：

RI<|z|<R2

2、離散時間系統的z域分析

系統函數

13

y(z)=X(z)*H(z)

y(z)

H(z)=

x(z)

H(z)稱作線性移不變系統的系統函數，而目.在單位圓z上的系統函數就是系統的

頻率響應。

利用系統函數的極點分布確定系統因果性與穩定性

線性移不變系統穩定的充要條件是h(n)必須滿足絕對可和：Z|h(n)|<8。

z變換H(z)的收斂域由滿足冽h(n)z-n|<8的那些z值確定。如單位圓上收斂，此時則有

Z|h(n)|<8,即系統穩定：也就是說，收斂域包括單位圓的系統是穩定的。

因果系統的單位抽樣響應為因果序列，其收斂域為R+〈|Z|W8；而因果穩定系統的系

統函數收斂域為lW|z|W8,也就是說，其全部極點必Z頁在單位圓內。

三、實驗用函數

1、ztrans

功能：求無限長序列的z變換。

調用格式：

X=ztrans(x)；用于求解無限長序列的z變換。

2、iztrans

功能：求函數X(z)的z反變換。

調用格式：

x=iztrans(X)

3、syms

功能：定義符號對象。

調用格式：

sysmt,i,x：將變量t,i.x聲明為符號變量。

4、rcsiducz

功能：有理多項式的部分分式展開。

調用格式

[r,pxl=residuez(b.a)；求解以系數向量b,a表示的系統函數的部分分式展開。

5、zplane

14

功能：繪制零極點分布圖。

調用格式

zplane(z,p)：繪制由列向量z確定的零點、列向量p確定的極點構成的零極點分布圖。

zplanc(h,a)：繪制由行向量b,a構成的系統函數確定的零極點分布圖。

[hz,hp,ht]=zplane(z,p)：獲得三個句柄向量：hz為零點線句柄；hp為極點線句柄；ht為

坐標軸、單位圓及文本對象的句柄。

四、參考實例

1、求序列X(7?)=--——的Z變換

2

MATLAB程序

symsn

x=(n*(n-1))/2;

X=ztrans(x)

程序執行結果：

X=1/2*z*(z+1)/(z-1)A3-1/2*z/(z-1)A2

7

2、求函數X(z)=-------r的z反變換

(Z-D2

MATLAB程序

symsz

X=z/(z-l)A3;

x=iztrans(X)

程序執行結果：

x=-l/2*n+l/2*nA2

3、用部分分式法求X(z)=「一三------(|z|>l)z反變換

Z2-1.5Z+0.5

將表達式整理為：

X(z)=---------

；---------7

+0.5z-2

MATLAB程序

b=[l]；

a=[1,-1,5,0.5];

15

[r,p,cl=residuez(b.a)

程序執行結果：

r=

2

-1

P=

1.0000

0.5000

[]

多項式分解后為：

X(z)=^T------~~r

l-z-,I-0.5Z-1

X(z)的反變換:

x(〃)=2.(.)一(0.5)〃￡(〃)

Z?Z

4、已知一個離散系統的函數H(z)=--------------,輸入序列X(z)二——，求系統

z-1.5z+0.5Z-1

在變換域的響應Y(z)及時間域的響應y(n)o

MATLAB程序

symsz;

X=z./(z-l);

H=z.A2./(z.A2-1.5*z+0.5);

Y=X.*H

y=iztrans(Y)

程序執行結果

Y=zA3/(z-l)/(zA2-3/2*z+1/2)

y=(l/2)An+2*n

5、用部分分式法求系統函數的z反變換，并用圖形與impz求得的結果相比較。

已知系統函數：

16

w、0.1321-0.3963z-2+0.3962Z-4

H(z)=------------------------------------------—

1+0.34319Z-2+0.60439ZT

MATLAB程序：

b=[0.1321,0,-0.3963,0,0.3962];

a=[1,0,0.34319.0,0.60439];

lr,p,c]=rcsiducz(b,a);

N=40;n=0:N-l;

h=r(i)*p(I).An+rC2)*p(2).An+r(3)*p(3).An4-r(4)*p(4).An4-c(1).*[n==01;

subplot(1,2,1);stem(n,real(h))

subplot(l,2,2);impzib,a,n)

程序運行結果如圖4-1所示。

Eil?EditVisInsertloolsDesktopIindowHelp

D4Q昌4爭⑥要□目sO

ImpulseResponse

0.5r0.5

o

0.40.4

圖4-1

從圖像上可以看出，用部分分式法求系統函數的z反變換也是一種求解系統的沖激響應

的有效方法。

6、系統極點的位置對系統響應的影響

1）研究z右半平面的實數極點對系統響應的影響。

己知系統的零-極點增益模型分別為:

"z）二二

17

MATLAB程序:

zl=[O]*;pl=[0.85r；k=l;

[bl,al]=zp2tf(zl,pl,k);

subplot(2,l,l);zplane(bl,al);

subplot(2,1,2);impz(bl,a1,20);

ylabelC極點在單位圓內)

程序運行結果如圖4-2所示

o5

0

5

.o

-10123

RealPart

ImpulseResponse

$

o

n(samples)

圖4-2

2）研究z右半平面的復數極點對系統響應的影響

已知系統函數:

z(z-0.3)

H,(z)=

(Z-0.5-0.77)(Z-0.5+0.7J)

求這些系統的零極點分布圖以及系統的沖激響應，井判斷系統的穩定性。

MATLAB程序:

z1=10,0.3]'

pl=[0.5+0.7jA5-0.7j]';

k=i;

[b1,a1]=zp2tf(z1,p1,k);

subplot(2,l,l);zpkuie(bl,al)；

subplot(2,l,2);impzibLal);

ylabelC極點在單位圓內，）

18

程序執行結果如圖4-3所示

=

e

d

n

e

w

一

o

RealPart

ImpulseResponse

102030405060

n(samples)

圖4-3

五、實驗任務

I、輸入并運行例題程序。

2、用部分分式法求解下列系統函數:

10+20Z-

X(z)=

l+8z-'4-l9z-24-12z-:

的z變換，寫出x(n)的表達式，并用圖形與impz求得的結果相比較。

3、已知離散時間系統函數分別為：

4-1.6Z-I-1.6Z-2+4Z-3

H,(z)=

1+0.4Z-,4-0.35Z-2-0.4Z-3

z-0.3

H(Z)=

2(z+i-j)(z+l+j)

求以上系統的零極點及零極點分布圖以及系統沖激響應，并判斷系統的因果穩定性。

六、實驗報告

1、簡述實驗日的、原理。

2、寫出上機調試通過的實驗任務的程序并描述其圖形曲線。

19

實驗五離散傅立葉級數

一、實驗目的

1、加深對離散周期序列傅里葉級數基本概念的理解。

2、掌握MATLAB求解周期序列傅里葉級數變換和逆變換的方法。

3、觀察離散周期序列的重復周期數對頻譜特性的影響。

二、實驗原理

離散時間序列x(n)滿足x(n)=x(n+rN),稱為離散周期序列，其中N為周期，x(n)為主值

序列。

周期序列可用離散傅里葉級數表示成

1N-\—kn

乳〃)=一Z又,=IDFS[X(k)]n=O/，..,N?l

N3

其中，元(3是周期序列離散傅里葉級數第K次諧波分鼠的系數，也稱為冏期序列的頻譜,

可表示為

N-\-i—kn

區(6=2只〃,N=z)FS[x(H)Jk=O,l，…,N-1

?=0

上面兩式是周期序列的一對傅里葉級數變換對。

.21

令叱\,二0'N,以上兩式可簡寫為:

N-\

寧(&)=。對[￡(〃)]=2雙〃)叱，

71=0

1N-\

乳〃)=/DFS[X(^)]=-X^(^AV

三、實驗用函數

1、mod

功能：模除求余。

調用格式：

mod(x,ni)：x整除m取正余數。

2、floor

功能：向-8舍入為整數。

20

調用格式:

floor(x)：將x向-8含入為整數。

四、實例

1、周期序列的傅里葉變換和逆變換

依據變換公式編寫通用函數

1)離散傅里葉級數正變換通用函數

functionxk=dfs(xn,N)

n=[O:l:N-l];%n的行向量

k=n;%k的行向量

WN=exp(-j*2*pi/N);%WN因子

nk=n'*k;%產生一個含nk值的N乘N維矩陣

WNnk=WN/nk;%DFS矩陣

xk=xn*WNnk;%DFS系數行向量

2)離散傅里葉級數逆變換通用函數

functionxn=idfs(xk,N)

n=[0:l:N-l];%n的行向量

k=n;%k的行向量

WN=exp(-j*2*pi/N);%WN因子

nk=n'*k;%產生一個含nk值的N乘N維矩陣

WNnk=WN/(-nk);%DFS矩陣

xn=xk*WNnk/N;%DFS系數行向量

例：已知一個周期性矩形序列的脈沖寬度占整個周期的1/4,一個周期的采樣點為16

點。用傅里葉級數求信號的幅度和相位頻譜；求傅里葉級數逆變換的圖形，與原信號性形進

行比較。

21

MATLAB程序

N=16;

xn=[ones(I,N/4),zeros(1,3*N/4)];

n=O:N-l;

xk=dfs(xn,N);

xnl=idfs(xk,N);

subplot(2,2,l);stem(n,xn);titlc('x(n)');

subplot(2,2,2);stem(n,abs(xn1));title('idfs(|X(k)|),);

subplot(2,2,3);stem(n,abs(xk));title('|X(k)r);

subplot(2,2,4)；stem(n,angle(xk));title('arg|X(k)r);

程序運行結果如圖4-1

圖4-1

2、周期重復次數對序列頻譜的影響

理論上講，周期序列不滿足絕對可積條件，要對周期序列進行分析，可以先取K個周

期進行處理，然后讓K無限增大，研究其極限情況。這樣可以觀察信號序列由非周期到周

22

期變換時，頻譜由連續譜逐漸向離散譜過渡的過程。

例：已知?個矩形序列的脈沖寬度占整個周期的1/2,?個周期的采樣點數為10,用傅

立葉級數變換求信號的重復周期數分別為1、4、7、10時的幅度頻譜。

MATLAB程序：

xn=[ones(1,5),zeros(1,5)];

Nx=lcngth(xn);

N\v=l000;dw=2*pi/Nw;

k=floor((-Nw/2+0.5):(Nw/2+0.5));

forr=0:3;

K=3*r+1;

nx=0:(K*Nx-l);

x=xn(mod(nx,Nx)+1);

Xk=x*(exp(-j*dw*nx**k))/K;

subplot(4,2,2*r+l);stcm(nx,x)

axis([0,K*Nx-l,0,11])；ylabclCxCn),);

subplot(4,2,2*r+2);piot(k*dw,abs(Xk))

axis([-4,4,0,1.1*max(abs(Xk))]);ylabel('X(k)');

end

程序運行結果如圖4-2

23

Figure1

FieEdftU?whsertToolsDesktopVMtxtowHe^>a

口團口昌洋aac⑥惠口日3國

圖4-2

從上圖可以看出，信號序列的周期數越多，則頻譜越是向幾個頻點集中，當信號周期

數趨于無窮人時，頻譜轉化為離散譜。

五、實驗任務

1、輸入并運行例題程序，熟悉基本指令的使用。

2、已知一個信號序列的主值為x(n)=[(M23,2,l,0]，顯示兩個周期的信號序列波形。

求解：用傅里葉級數求信號的幅度和相位頻譜；求傅里葉級數逆變換的圖形，與原信號

圖形進行比較。

六、實驗報告

1、簡述實驗日的、原理。

2、寫出上機調試通過的實驗任務的程序并描述其圖形曲線。

24

實驗六離散傅里葉變換

一、實驗目的

1、加深對離散傅里葉變換基本概念的理解。

2、了解有限長序列傅里葉變換與周期序列傅里葉級數的聯系。

3、熟悉相關函數的使用方法。

二、實驗原理

有限長序列的傅里葉變換和逆變換

對于非周期序列，在實際中常常使用有限長序列。有限長序列x(n)表示為

x(n)0<n<N-\

0n為其它值

x(n)是非周期序列，但可以理解為某一周期序列的主值序列。由離散傅立葉級數DFS

和IDFS引出有限長序列的離散傅立口I正、逆變換關系式。

13聲切

x(n)=IDFT[x(k)]=-Xx(k)eN

DFT與DFS的關系

比較兩者的變換對，可以看出兩者的區別僅僅是將周期序列換成了有限長序列。

有限長序列x(n)可以看作是周期序列雙〃)的?個周期；反之周期序列取〃)可以看作是

有限長序列x(n)以N為周期的周期延拓。

由于公式非常相似，在程序編寫上也基本一致。

三、實例

1、已知有限長序列x(n)為：

x(n)=[0,l,2,3,4,5,6,7,8,9],求x(n)的DFT和IDFT。要求

1)畫出序列傅里什變換對應的|X(k)|和arg[X(k)]圖形。

2)畫出原信號與傅里葉逆變換【DFT[X(k)]圖形進行比較。

25

MATLAB程序:

xn=[0,l,2,3,4,5,6,7,8,9];N=length(xn);

n=O:N-i;k=O:N-l;

xk=xn*exp(-j*2*pi/N).A(n,*k);

x=(xk*exp(j*2*pi/N).A(n,*k))/N;

subplot(2,2,1);stcm(n,xn);Ut]cCx(n)，);

subplot(2,2,2);stem(n,abs(x));title(*IDFT|X(k)r);

subplot(2,2,3);stem(k,abs(xk));title。|X(k)「);

subplot(2,2,4);stem(k,angle(Xk));title(4arg|X(k)r);

程序運行結果如圖5-1：

圖5-1

由上圖可看出，與周期序列不同的是，有限長序列本身是僅有N點的離散序列，相當

于周期序列的土值部分。因此，其頻譜也對應序列的土值部分，是含N點的離散序列,

26

2、有限長序列DFT與周期序列DFS的聯系

已知周期序列的主值x(n尸。1,2,3,4,5],求x(n)周期重復次數為4次時的DFS。要求

1)畫出原主值序列和信號周期序列；

2)畫出序列傅里葉變換對的圖形。

MATLAB程序：

xn=[O,l23,4,5];N=lcngth(xn);

n=0:4*N-l;k=0:4*N-l;

xnl=xn(mod(n,N)+1);

xk=xnl*exp(-j*2*pi/Ni.A(n,*k);

subplot(2,2,l);stem(xn);title('原主值信號x(n)');

subplot(2,2,2);stemin,xnl);title('周期序列信號')；

subplot(2,2,3);stem(k,abs(xk));title('|X(k)r);

subpiot(2,2,4);stem(k,angle(xk));title('arg|X(k)r);

程序運行結果如圖5-2：

圖5-2

27

與上一個例題比較，有限長序列x(n)可以看成是周期序列的一個周期，反之，周期序列

可以看成是有限長序列以N為周期的周期延拓。領域上的情況也是相同的。從這個意義上

說，周期序列只是有限個序列值有意義。

四、實驗任務

1.輸入并運行例題程序，熟悉基木指令的使用“

2、驗證離散傅里葉變換的線性性質。

有兩個有限長序列分別為xl(n)和x2(n),長度分別為N1和N2,且y(n)=ax1(n)+bx2(n),

(a,b均為常數)，則該y(n)的N點DFT為

Y(k)=DFT[y(n)]=aXl(k)+bX2(k)(0<=k<=N-l)

其中：N=max(NLN2),Xl(k)和X2(k)分別為xl(n)和x2(n)的N點DFT。

已知序列:xl(n>[0,l,2,4J

x2(n)=ll,0,l,0JJ

五、實驗報告

1、簡述實驗目的、原理。

2、寫出上機調試通過的實驗任務的程序并描述其圖形曲線。

28

實驗七快速傅里葉變換

一、實驗目的

1、加深對快速傅里葉變換基本理論的理解。

2、了解用MATLAB語言進行快速傅里葉變換的方法。

3、掌握常用函數的調用方法。

二、實驗原理

DFT是在時域和頻域均為離散序列的變換方法，它適用于有限長序列。但如果按照變

換公式進行運算的話，當序列長度很大時，將占用很大的內存空間，運算時間也會很長，無

法實時處理問題。

快速傅里葉變換是用于提高DFT運算的高速運算方法的統稱，FFT是其中的一種，FFT

不是一種新的變換形式，它僅僅只是一種快速算法。FFT主要有時域抽取算法和頻域抽取算

法，基本思想是將一個長度為N的序列分解成多個短序列再進行運算，如基2算法、基4

算法等等，從而可以大大縮短運算時間。

三、實驗用函數

1、fft

功能：一維基2快速傅里葉變換。

調用格式：

y=fft(x)：利用FFT算法計算矢量x的離散傅里葉變換。當x為2的事次方時，采用高

速基2FFT算法，否則為稍慢的混合算法。

y=fft(x,N)：采用n點FFT。當x的長度小于N時，FFT函數在x的尾部補零，以構成

N點數據，當x的長度大于N時，FFT函數在x的尾部截斷，以構成N點數據。

2、ifft

功能：一維基2快速傅里葉逆變換。

調用格式：

與FFT調用方法相同，只需改換函數名。

3、fftshifl

功能：對FFT的輸出進行重新排列，將零頻分量移到頻譜的中心。

調用格式：

y=fft?;hift(x)I對FFT的輸出進行重新排列,將零頻分量移到頻譜的中心°

29

四、實例

1、用MATLAB工具箱函數FFT進行頻譜分析時福要注意：

1)函數fft的返回值y的數據結構的對稱性

若已知序列x=[4,3,2,6,7,8,9,0],求X(k)=DFT[x(n)]

利用函數fft計算，其MATLAB程序如下：

N-8,

n=O:N-l;

xn=[4,3,2,6,7,8,9,0];

xk二fft(xn)'

程序運行結果如下：

xk=

39.0000

-10.7782-6.2929i

0+5.0()00i

4.7782+7.707li

5.0000

4.7782-7.707li

0-5.0000i

-10.7782+6.2929i

由程序運行結果可見，xk的第一行元素對應頻率值為0,第五行元素對應頻率為萊奎

斯特(Nyquist)頻率，即標準頻率值為I。因此第一行至第五行對應的標準頻率為而

第五行至第八行對應的是負頻率，其K(x)值是以Nyquist頻率為軸對稱。

一般情況，對于N點的x(n)序列的FFT是N點的復數序列，其點n=N/2+l對應Nyquist

頻率，作譜分析時僅取序列X(k)的前一半即可，其后一半序列和前一半是對稱的。

2)頻率計算

若N點序列x(n)(n=0,l,…,N-1)是在采樣頻率fs(Hz)下獲得的。它的FFT也是N點序列，

即X(k)(k=0,l…則第K點對應實際頻率值為：

f=k*fs/N

3)作FFT分析時，幅值大小與FFT選擇的點數有關，但不影響分析結果。

30

2、已知信號由15Hz幅值0.5的正弦信號和40Hz幅值2的正弦信號組成，數據采樣頻

率為100Hz,試繪制N=128點DFT的幅頻圖。

MATLAB程序如下：

fs=100;

N=128;

n-0;N-I;

l=n/fs;

x=0.5*sin(2*pi*l5*l)+2*sin(2*pi*40*l);

y=fft(x,N);

f=(0:length(y)-1)'*fs/length(y);

mag=abs(y);

stem(f,mag);

title(*N=128點')

程序運行的結果如圖6/

圖6-1

如圖所示，由于信號采樣頻率為100Hz,故其萊奎斯特頻率為50Hz,圖中整個頻譜圖

是以萊奎斯特頻率為軸對稱的。因此利用FFTXj信號作譜分析時，只要考察0~Nyquist頻率

范圍的幅頻特性就可以了,

31

3、利用FFT進行功率譜的噪聲分析

已知帶有測量噪聲信號x(t)=sin(2/z\")+sin(2乃f2t)+2"(，)其中f1=50Hz,

f2=120Hz,。⑺為均值為零、方差為1的隨機信號，采樣頻率為1000Hz,數據點數N=512。

試繪制信號的功率譜圖。

MATLAB程序如下

t-00001:0.6,

x=sin(2*pi*50*l)+sin(2*pi*120*1);

y=x+2*randn(1,length^));

Y=fft(y,5⑵；

P=Y.*conj(Y)/512;%求功率

f=1000*(0:255)/512;

subplot(2,l,l);

plot(y);

subplot(2,1,2);

plol(f,P(1:256));

程序運行結果如圖6-2

圖6-2

32

4、序列長度和FFT的長度對信號頻譜的影響。

已知信號x(7)=0.5sin(2不工，)+2sin(2;rf/)

其中fl=15Hz,f2=40Hz,采樣頻率為100Hz.

在下列情況卜繪制其幅頻譜。

Ndata=32,Nfft=32;

Ndata-32,Nfft-l28.

MATLAB程序如下

fs=100;

Ndata=32;Nfft=32;

n=0:Ndata-l;

t=n/fs;

x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);

y=fft(x,Nfft);

mag=abs(y);

f=(0:length(y)-1)**fs/length(y);

subplot(2,l,l)

plot(f(1:Nfft/2),mag(1:Nfft/2))

titleCNdata=32,Nfft=32')

程序運行結果如圖6-3所示

33

5、快速卷積的FFT算法

在MATLAB實現卷積的函數為CONV,對于N值較小的向量，這是十分有效的。對于

N值較大的向晟卷積可用FFT加快計算速度。

由DFT性質可知，若DFT[xl(n)]=Xl(k),DFT[x2(n)]=X2(k)則

王(〃)*W(〃)=IDFT'Xi(k)?X2(k)]=IDFT[DFT[xx(n)VDFT[x2(n)]\

若DFT和IDFT均采用FFT和IFFT算法，可提高卷積速度。

五、實驗任務

1、輸入并運行例題程序，熟悉基本指令的使用。

2、比較定義式計算傅里葉變換和用快速算法計算傅里葉變換所用時間。

3、比較卷積函數與快速卷積運算所用時間。

(提示：clock函數讀取瞬時時鐘

elime⑴,12)函數計