1.1市場與市場營銷1.2我國汽車市場的發展與現狀復習思考題實驗13離散傅里葉變換的性質一、實驗目的

(1)加深對離散傅里葉變換(DFT)基本性質的理解。

(2)了解有限長序列傅里葉變換(DFT)性質的研究方法。

(3)掌握用MATLAB語言進行離散傅里葉變換性質分析時程序編寫的方法。二、實驗原理

1.線性性質

如果兩個有限長序列分別為x1(n)和x2(n)，長度分別為N1和N2，且

y(n)＝ax1(n)＋bx2(n)(a、b均為常數)

則該y(n)的N點DFT為

Y(k)＝DFT［y(n)］＝aX1(k)＋bX2(k)0≤k≤N－1

其中：N＝max［N1，N2］，X1(k)和X2(k)分別為x1(n)和x2(n)的N點DFT。

例13-1

已知x1(n)＝［0，1，2，4］，x2(n)＝［1，0，1，0，1］，求：

(1)y(n)＝2x1(n)＋3x2(n)，再由y(n)的N點DFT獲得Y(k)；

(2)由x1(n)、x2(n)求X1(k)、X2(k)，再求Y(k)＝2X1(k)＋3X2(k)。

用圖形分別表示以上結果，將兩種方法求得的Y(k)進行比較，由此驗證有限長序列傅里葉變換(DFT)的線性性質。

解

MATLAB程序如下：

xn1＝［0，1，2，4］；%建立xn1序列

xn2＝［1，0，1，0，1］；%建立xn2序列

N1＝length(xn1)；N2＝length(xn2)；

N＝max(N1，N2)；%確定N

ifN1>N2xn2＝［xn2，zeros(1，N1－N2)］；%對長度 短的序列補0

elseifN2>N1xn1＝［xn1，zeros(1，N2－N1)］；

end

yn＝2*xn1＋3*xn2；%計算yn

n＝0：N－1；k＝0：N－1；

Yk1＝yn*(exp(－j*2*pi/N)).^(n￠*k)；%求yn的N點DFT

Xk1＝xn1*(exp(－j*2*pi/N)).^(n￠*k)；%求xn1的N點DFT

Xk2＝xn2*(exp(－j*2*pi/N)).^(n￠*k)；%求xn2的N點DFT

Yk2＝2*Xk1＋3*Xk2；%由Xk1、Xk2求Yk

以上程序作圖部分省略。用兩種方法求得的Y(k)結果一致，如下所示：

Yk＝

23.0000 －7.5902＋1.5388i 3.5902－0.3633i

3.5902＋0.3633i－7.5902－1.5388i

運行結果如圖13-1所示。

圖13-1例13-1有限長序列的傅里葉變換的線性性質

2.循環移位性質

如果有限長序列為x(n)，長度為N，將x(n)左移m位，則

y(n)＝x((n＋m)N)RN(n)

x(n)左移m位的過程可由以下步驟獲得：

(1)將x(n)以N為周期進行周期延拓，得到 ＝x((n)N)；

(2)將 左移m位，得到 ；

(3)取 的主值序列，得到x(n)循環移位序列y(n)。有限長序列的移位也稱為循環移位，原因是將x(n)左移m位時，移出的m位又依次從右端進入主值區。下面舉例說明。

例13-2

已知有限長序列x(n)＝［1，2，3，4，5，6］，求x(n)左移2位成為新的向量y(n)，并畫出循環移位的中間過程。

解MATLAB程序如下：

xn＝［1，2，3，4，5，6］；%建立xn序列

Nx＝length(xn)；nx＝0：Nx－1；

nx1＝－Nx：2*Nx－1；%設立周期延拓的范圍

x1＝xn(mod(nx1，Nx)＋1)；%建立周期延拓序列

ny1＝nx1－2；y1＝x1；%將x1左移2位，得到y1

RN＝(nx1>＝0)&(nx1<Nx)；%在x1的位置向量nx1上設置主值窗

RN1＝(ny1>＝0)&(ny1<Nx)；%在y1的位置向量ny1上設置主值窗

subplot(4，1，1)，stem(nx1，RN.*x1)；%畫出x1的主值部分

subplot(4，1，2)，stem(nx1，x1)；%畫出x1

subplot(4，1，3)，stem(ny1，y1)；%畫出y1

subplot(4，1，4)，stem(ny1，RN1.*y1)；%畫出y1的主值部分

運行結果如圖13-2所示。

圖13-2例13-2有限長序列的循環移位

3.循環折疊性質

如果要把有限長N點序列x(n)直接進行折疊，則x的下標(－n)將不在0≤n≤N－1區域內。但根據有限長序列傅里葉變換隱含的周期性，可以對變量(－n)進行N求余運算。即在MATLAB中，序列x(n)的折疊可以由y＝x(mod(－nx，N)＋1)得到。

有限長N點序列x(n)的循環折疊序列y(n)定義為

可以想像成，序列x(n)以反時針方向等間隔放置在一個圓周上，則x(－n)是將x(n)沿著圓周順時針方向等間隔放置。

循環折疊性質同樣適用于頻域。經循環折疊后，序列的DFT由下式給出：

就是說，在時域循環折疊后的函數，其對應的DFT在頻域也作循環折疊，并取X(k)的共軛。

例13-3

求x(n)＝［1，2，3，4，5，6，7］，循環長度分別取N＝7，N＝10。

(1)畫出x(n)的圖形；

(2)畫出x(－n)的圖形。

解

MATLAB程序如下：

x1＝［1，2，3，4，5，6，7］；%建立x(n)，N ＝7序列

N1＝length(x1)；n1＝0：N1－1；

y1＝x1(mod(－n1，N1)＋1)；%建立x(－n)，N＝7序列

N2＝10；

x2＝［x1，zeros(1，N2－N1)］；%建立x(n)，N＝10 序列

n2＝0：N2－1；

y2＝x2(mod(－n2，N2)＋1)；%建立x(－n)，N＝10序列

subplot(2，2，1)，stem(n1，x1，￠k￠)；%畫x(n)，N＝7

title(￠x(n)，N＝7￠)；

subplot(2，2，3)，stem(n1，y1，￠k￠)；%畫x(－n)，N ＝7

title(￠x(－n)，N＝7￠)；

subplot(2，2，2)，stem(n2，x2，￠k￠)；%畫x(n)，N＝10

title(￠x(n)，N＝10￠)；

subplot(2，2，4)，stem(n2，y2，￠k￠)；%畫x(－n)，N ＝10

title(￠x(－n)，N＝10￠)；

運行結果如圖13-3所示。

圖13-3例13-3離散序列的循環折疊

例13-4

如例13-3求x(n)＝［1，2，3，4，5，6，7］，循環長度取N＝7。求證：在時域循環折疊后的函數x(－n)，其對應的DFT在頻域也作循環折疊，并取X(k)的共軛。

解

MATLAB程序如下：

x1＝［1，2，3，4，5，6，7］；%建立x(n)，N＝ 7序列

N＝length(x1)；

n＝0：N－1；k＝0：N－1；

y1＝x1(mod(－n，N)＋1)；%建立x(－n)，N＝7序列

Xk＝x1*exp(－j*2*pi/N).^(n￠*k)%求x(n)的DFT

Yk＝y1*exp(－j*2*pi/N).^(n￠*k)%求x(－n)的DFT

運行結果：

Xk＝

28.0000－3.5000＋7.2678i－3.5000＋2.7912i－3.5000＋0.7989i

－3.5000－0.7989i－3.5000－2.7912i－3.5000－7.2678i

Yk＝

28.0000－3.5000－7.2678i－3.5000－2.7912i－3.5000－0.7989i

－3.5000＋0.7989i－3.5000＋2.7912i－3.5000＋7.2678i

4.時域和頻域循環卷積特性

離散傅里葉變換的循環卷積特性也稱為圓周卷積，分為時域卷積和頻域卷積兩類。

1)時域循環卷積

假定x(n)、h(n)都是N點序列，則時域循環卷積的結果y(n)也是N點序列：

若x(n)、h(n)和y(n)的DFT分別為X(k)、H(k)和Y(k)，則

Y(k)＝X(k)H(k)

2)頻域循環卷積

利用時域和頻域的對稱性，可以得到頻域卷積特性。若

y(n)＝x(n)h(n)

則

下面重點討論時域循環卷積。時域循環卷積的方法有多種：

方法1：直接使用時域循環卷積。

由于有限長序列可以看成是周期序列的主值，因此，時域圓周卷積的結果可以由對應的周期序列卷積和取主值部分獲得，請參考例11-4。

方法2：用頻域DFT相乘再求逆變換。

即先分別求x1(n)、x2(n)的DFTX1(k)、X2(k)，再求Y(k)的IDFT獲得y(n)。基本思路如圖13-4所示。

圖13-4用DFT實現循環卷積的框圖方法3：用FFT和IFFT進行循環卷積。

基本思路同方法2，但直接使用了MATLAB提供的fft和ifft子函數來實現。見后面的快速傅里葉變換實驗。例13-5將例11-4已知的兩個時域周期序列分別取主值，得到x1＝［1，1，1，0，0，0］，x2＝［0，1，2，3，0，0］，求時域循環卷積y(n)并用圖形表示。

解本例采用方法2。程序如下(作圖程序部分省略)：

xn1＝［0，1，2，3，0，0］；%建立x1(n)序列

xn2＝［1，1，1，0，0，0］； %建立x2(n)序列

N＝length(xn1)；

n＝0：N－1；k＝0：N－1；

Xk1＝xn1*(exp(－j*2*pi/N)).^(n￠*k)；%由x1(n)的DFT求X1(k)

Xk2＝xn2*(exp(－j*2*pi/N)).^(n￠*k)；%由x2(n)的DFT 求X2(k)

Yk＝Xk1.*Xk2；%Y(k)＝X1(k)X2(k)

yn＝Yk*(exp(j*2*pi/N)).^(n￠*k)/N；%由Y(k)的IDFT求 y(n)

yn＝abs(yn)%取模值，消除DFT帶來的微小復數影響

得到：

yn＝

0.00001.00003.00006.00005.00003.0000

運行結果如圖13-5所示。由y(n)圖形可見，與例11-4主值區域的卷積結果相同。

圖13-5例13-5離散序列時域循環卷積的結果

5.循環對稱性

由于序列x(n)及其離散傅里葉變換X(k)的定義在主值為0～N－1的區間，因此DFT的循環對稱性對時間序列是指關于n＝0和n＝N/2的對稱性，對頻譜序列是關于數字頻率為0和p的對稱性。

本實驗重點分析實序列的循環對稱性。

實序列x(n)可以分解為循環偶序列xe(n)和循環奇序列xo(n)：

x(n)＝xe(n)＋xo(n) 0≤n≤N－1其中：

設DFT［x(n)］＝X(k)＝Re［X(k)］＋j*Im［X(k)］，則有

即實序列中的偶序列xe(n)對應于x(n)的離散傅里葉變換X(k)的實部，而實序列中的奇序列xo(n)對應于x(n)的離散傅里葉變換X(k)的虛部。

例13-6

已知一個定義在主值區間的實序列x＝［ones(1，4)，zeros(1，4)］，試將其分解成為偶對稱序列和奇對稱序列，并求它們的DFT，驗證離散傅里葉變換的循環對稱性。

解程序如下(作圖程序省略)：

x＝［ones(1，5)，zeros(1，5)］%建立x(n)序列

N＝length(x)；

n＝0：N－1；k＝0：N－1；

xr＝x(mod(－n，N)＋1)；%求x(－n)

xe＝0.5*(x＋xr)%求x(n)的偶序列

xo＝0.5*(x－xr)%求x(n)的奇序列

X＝x*(exp(－j*2*pi/N)).^(n￠*k)；%由x(n)的DFT求X(k)

Xe＝xe*(exp(－j*2*pi/N)).^(n￠*k)；%由xe(n)的DFT求 Xe(k)

Xo＝xo*(exp(－j*2*pi/N)).^(n￠*k)；%由xo(n)的DFT求 Xo(k)

error1＝(max(abs(real(X)－Xe)))%計算X(k)的實部與 Xe(k)的差值

error2＝(max(abs(j*imag(X)－Xo)))%計算X(k)的虛部與 Xo(k)的差值

運行結果顯示：

x＝

1

1

1

1

0

0

0

0

xe＝

1.0000

0.50000.50000.5000

00.50000.50000.5000

xo＝

00.5000

0.50000.5000

0－0.5000－0.5000－0.5000

error1＝7.2173e－015

error2＝7.4517e－015

程序執行結果如圖13-6所示。

圖13-6例13-6驗證離散實序列的循環對稱性由以上輸出數據和圖形可知：

(1)xe(n)具有循環對稱性。對稱中心