陳章寶,鄧運生

(蚌埠學(xué)院電子與電氣工程學(xué)院,安徽 蚌埠 233030)

1 問題的提出

基于離散傅里葉變換(DFT)和快速傅里葉變換(FFT)的信號頻譜分析、信號處理算法是數(shù)字信號處理課程的重要教學(xué)內(nèi)容。傅里葉變換是數(shù)學(xué)史上十大經(jīng)典的算法之一,是信號處理類課程如“信號與系統(tǒng)”“數(shù)字信號處理”“數(shù)字圖像處理”“語音信號處理”理論分析和應(yīng)用開發(fā)不可或缺的重要工具[1]。DFT在信號及其頻譜分析、系統(tǒng)及其頻率特性分析、系統(tǒng)響應(yīng)求解等教學(xué)中被廣泛應(yīng)用[2],是進行工程應(yīng)用開發(fā)的基礎(chǔ)。王靜等[3]從信號時頻域的對應(yīng)關(guān)系入手,分析了傅里葉變換各個性質(zhì)之間的關(guān)系。姜恩華等[4]借助時域采樣定理和頻域采樣定理,分析了各種信號傅里葉變換之間的關(guān)系和推理過程。趙潔等[5]以二維圖像信號為分析對象,給出了圖像頻譜與內(nèi)容之間對應(yīng)關(guān)系的教學(xué)案例設(shè)計。

在數(shù)字信號處理課程中對傅里葉變換的介紹非常詳盡透徹,課程內(nèi)容包括四種傅里葉變換(CTFT、CTFS、DTFT、DTFS)、離散傅里葉變換(DFT)、快速傅里葉變換(FFT),算法的代碼實現(xiàn)到工程應(yīng)用[6],構(gòu)成了一個完整的體系。本文從連續(xù)周期信號的傅里葉級數(shù)引入四種信號(連續(xù)周期、連續(xù)非周期、離散周期、離散非周期)的傅里葉變換,再到工程中普遍應(yīng)用的DFT及其快速算法FFT,給出其理論推導(dǎo)和代碼實現(xiàn)的教學(xué)案例,最后給出傅里葉變換在語音信號頻譜分析、數(shù)字圖像信號頻譜分析的教學(xué)案例設(shè)計。

2 由傅里葉級數(shù)引入四種傅里葉變換

2.1 連續(xù)周期信號的傅里葉級數(shù)

首先,從連續(xù)周期信號這一簡單概念開始,“連續(xù)周期信號(方波信號、三角波信號等)可以分解為不同頻率(直流分量、基波分量、各次諧波分量)正弦信號的加權(quán)疊加”是學(xué)生初次接觸的傅里葉變換的概念。以周期方波信號(周期T=1 s,脈沖寬度為0.5,幅度A=2)為例,如圖1(a)所示,其三角函數(shù)形式的傅里葉級數(shù)展開式見式(1),其模擬角頻率ω0=2×π/t,各正弦項的系數(shù)即為其傅里葉級數(shù)。

(a)方波信號

(1)

借助MATLAB進行可視化,方波信號及其合成圖如圖1所示。

通過仿真案例分析,學(xué)生可以直觀地看到方波信號可以分解為直流分量、基波分量和各次諧波分量,以及合成原信號過程的可視化效果。隨著諧波次數(shù)的增加,合成波形越來越接近于原信號波形。此案例也是介紹Gibbs現(xiàn)象的最佳案例,如圖1(b)和(c)所示,在信號不連續(xù)點處,有限項諧波合成出現(xiàn)了約9%的起伏峰值,說明在連續(xù)周期信號的不連續(xù)點處,只能滿足能量意義下的收斂,不滿足一致收斂。

2.2 四種傅里葉變換

傅里葉變換給出了信號的時域和頻域之間的變換關(guān)系,嚴(yán)格意義上說,自然界的信號皆為連續(xù)非周期信號,某些具有節(jié)奏感或頻率相對固定的信號也可以看作是周期信號,現(xiàn)代信號處理系統(tǒng)都是由計算機實現(xiàn)的,數(shù)字信號處理系統(tǒng)需要對自然信號進行采樣,得到離散信號并進行處理。根據(jù)信號在時域的四種表現(xiàn)形式,對應(yīng)有四種不同的傅里葉變換,分別為連續(xù)時間傅里葉變換(Continuous Time Fourier Transform,CTFT)、連續(xù)時間傅里葉級數(shù)(Continuous Time Fourier Series,CTFS)、離散時間傅里葉變換(Discrete Time Fourier Transform,DTFT)、離散時間傅里葉級數(shù)(Discrete Time Fourier Series,DTFS),其傅里葉正變換如式(2-a)至(2-d)所示。

(2-a)

(2-b)

(2-c)

(2-d)

在教學(xué)過程中,學(xué)生需要掌握四種變換的數(shù)學(xué)概念、物理意義,才能深刻理解傅里葉變換,從數(shù)學(xué)角度考慮,傅里葉變換實現(xiàn)信號的時域表達函數(shù)x(·)和頻域表達函數(shù)X(·)的相互轉(zhuǎn)換,時域函數(shù)自變量為t(連續(xù)信號)或n(離散信號),頻域函數(shù)自變量為Ω(模擬角頻率)、ω(數(shù)字角頻率)或k(代表離散,實際為kΩ0或kω0);每個變換的核函數(shù)不同,核函數(shù)Φ(·)是一個包含了時域函數(shù)和頻域函數(shù)兩個函數(shù)自變量的函數(shù);傅里葉變換表達式為積分或求和運算,如果公式右邊為連續(xù)函數(shù),則傅里葉變換表達式為積分運算,如果公式右邊為離散函數(shù),則傅里葉變換表達式為求和運算。

在物理意義上,四種傅里葉變換不是并行關(guān)系,而是依次出現(xiàn)的,其解析過程如圖2所示,傅里葉變換的大部分知識點都可以聯(lián)系該圖進行講解。

圖2 五種傅里葉變換關(guān)系圖

(ii)CTFT到DTFT,信號的CTFT變換中,信號的時域函數(shù)和頻域函數(shù)皆為連續(xù)函數(shù),要用計算機對信號的時域和頻域函數(shù)進行分析,必須對連續(xù)函數(shù)離散化,時域函數(shù)抽樣前要求進行抗混疊低通濾波,這里進行時域抽樣定理分析;信號時域離散化導(dǎo)致其頻域函數(shù)周期化,所以頻域函數(shù)X(ejω)為一個周期為2π的連續(xù)函數(shù),其傅里葉變換為DTFT,如式(2-c)所示。

3 離散傅里葉變換仿真案例教學(xué)

3.1 離散傅里葉變換

(3)

DFT運算的時域和頻域函數(shù)皆為長度為N點的序列,為傅里葉變換的計算機實現(xiàn)提供了極大的方便和可能。自然信號傅里葉變換的計算機實現(xiàn)的本質(zhì)是DTFS,實現(xiàn)的方法是DFT,DFT運算隱含周期性,這在DFT運算原理和性質(zhì)教學(xué)中是尤其要強調(diào)的內(nèi)容。

3.2 公式表達與代碼實現(xiàn)

DFT運算中時域函數(shù)x(n)和頻域函數(shù)X(k)皆為長度為N的離散序列,DFT運算如式(3)所示,可以通過矩陣運算實現(xiàn),如式(4)所示。

X=WN·x,

(4)

其中,WN為N點DFT運算的變換矩陣。

由DFT變換表達式到變換矩陣有一系列的推導(dǎo)過程,教學(xué)過程中要進行嚴(yán)格的理論推導(dǎo),讓學(xué)生了解由DFT運算表達式到程序?qū)崿F(xiàn)的過程,這里直接給出其變換矩陣,如式(5-a)和(5-b)所示。

(5-a)

(5-b)

根據(jù)式(5-a)(5-b),進行DFT運算的MATLAB編程,矩陣運算效率高,代碼非常簡潔,示例代碼及其注釋如下:

function [xk] = dft(xn) %DFT函數(shù)定義

N = length(xn); %獲取x長度

N = 0∶1∶N-1; %獲取變換矩陣WNnk

K = 0∶1∶N-1;

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

Nk = n’*k;

WNnk = WN.^nk;

Xk = WNnk*xn; %矩陣乘法實現(xiàn)DFT運算

4 快速傅里葉變換仿真案例教學(xué)

4.1 快速傅里葉變換

(a)x(n)波形

圖3中序列取前200點進行顯示,圖3(a)為時域的隨機序列,圖3(b)為通過矩陣乘法實現(xiàn)DFT獲得的幅度譜,圖3(c)為FFT運算獲得的幅度譜,可見隨機信號頻譜在整個頻率軸上為均勻頻譜。通過對比實驗可以發(fā)現(xiàn),圖3(b)和(c)的波形完全相同,說明DFT和FFT運算結(jié)果完全相同。對運算進行計時結(jié)果顯示,DFT運算時間為34.397 1 s,而FFT運算時間為0.000 268 5 s,運算速度大大提升。通過運算速度對比實驗,使學(xué)生對FFT如何實現(xiàn)快速運算產(chǎn)生了極大的興趣。

4.2 運算流圖到代碼實現(xiàn)

FFT運算理論推導(dǎo)復(fù)雜,利用核函數(shù)的周期性、可約性、共軛對稱性和一些特殊點的值,實現(xiàn)將長序列變?yōu)槎绦蛄械腄FT運算,以降低運算次數(shù),核心運算是蝶形運算,在FFT知識點教學(xué)中重點介紹DIT-FFT算法推理及其代碼實現(xiàn)過程,N=8點序列的DIT-FFT運算的信號流如圖4所示。

圖4 DIT-FFT蝶形運算流圖

function xk = ditfft(xn) %DIT-FFT函數(shù)定義

%初始化:補零,倒位序,定義WN

for m = 1∶L %第m級運算

for n =1∶ (2.^(L-m)) %蝶形塊運算

for o = 1∶2.^(m-1)%蝶形運算

循環(huán)體:計算參數(shù)r

循環(huán)體:蝶形運算

end

end

end

課堂教學(xué)中,可以設(shè)計如上文3.1節(jié)的相似實驗,進行DFT、DIT-FFT、FFT的對比實驗,實驗結(jié)果相同,計算時間分別為34.397 1 s、0.000 865 s、0.000 268 s,說明MATLAB內(nèi)置函數(shù)FFT運算速度最快,借此引導(dǎo)學(xué)生學(xué)習(xí)基4FFT算法、混合基FFT算法等更加快速高效的傅里葉變換算法。

5 應(yīng)用教學(xué)案例設(shè)計

5.1 聲音信號及其頻譜分析

傅里葉變換在信號及其頻譜分析中應(yīng)用廣泛,綜合性仿真案例教學(xué)中,對于一維信號可采用語音信號及其頻譜分析的仿真案例教學(xué)。運用MATLAB函數(shù)audiorecorder()驅(qū)動計算機聲卡進行語音信號采集[8],分別采集了男聲和女聲“好了”兩個字的語音信號,采樣頻率為11 025 Hz,截取采樣點數(shù)為7 500點,利用FFT求其頻譜,信號時域波形及其頻譜如圖5所示,圖中顯示0～11 025/2(折疊頻率)之間的頻譜。

(a)女聲時域波形

從圖5可以看出,男女聲的時域波形沒有什么區(qū)別,但是從其頻譜圖可知,男女聲信號的頻譜主要分布于50～2 600 Hz之間,女聲信號的頻譜主要分布在1 000 Hz以上,分布范圍較寬,而男聲信號的頻譜主要分布在1 000 Hz以下,分布范圍較窄,3 000 Hz以上能量較低的部分頻譜應(yīng)為語音錄制時的噪聲,可以通過低通濾波濾除。

5.2 圖像信號頻譜分析

灰度圖像是典型的二維離散信號,二維信號的傅里葉變換的實質(zhì)是分別在每一個維度上進行一維FFT運算,二維傅里葉變換經(jīng)常用于圖像譜分析,圖像譜是將圖像灰度值分解為不同頻率和方向的正弦平面波的加權(quán)疊加。圖像的空間域增強是實現(xiàn)圖像處理、特征提取和模式識別的基礎(chǔ),傅里葉變換實現(xiàn)了從頻率域進行圖像增強的重要手段。圖像頻譜的低頻部分反映了圖像的背景,即灰度變化較為緩慢的部分,高頻部分反映了圖像的邊緣和細節(jié),即灰度變換較快的部分;圖像幅度譜反映了圖像的灰度信息,而相位譜反映了圖像的結(jié)構(gòu)信息[9],這些都可以通過如圖6所示的教學(xué)案例進行說明。

圖6 圖像幅度譜和相位譜關(guān)系

通過交換人和貓的圖像的相位譜合成新的圖像,其案例實現(xiàn)過程為:①用二維傅里葉變換分別獲取人和貓圖像的頻譜;②提取兩個圖像的幅度譜和相位譜;③交換相位譜,合成兩個新的圖像頻譜矩陣;④對兩個圖像頻譜矩陣進行二維傅里葉反變換,取絕對值;⑤顯示變換后的圖像如圖6(c)和(d)所示。由圖6可以明顯看出,相位譜決定了圖像的結(jié)構(gòu)信息。這也可以讓我們聯(lián)想到,在深度學(xué)習(xí)用卷積神經(jīng)網(wǎng)絡(luò)進行貓狗分類的實驗中,可以通過圖像的相位譜矩陣作為輸入信號進行貓狗的分類識別。

6 結(jié)語

離散傅里葉變換是數(shù)字信號處理類課程的重要教學(xué)內(nèi)容,需要從數(shù)學(xué)概念、物理意義和工程應(yīng)用等方面深刻理解傅里葉變換及其在信號頻譜分析中的應(yīng)用。本文通過理論分析并結(jié)合代碼實現(xiàn),運用MATLAB軟件進行一系列可視化教學(xué)案例設(shè)計。通過仿真案例教學(xué),使學(xué)生掌握傅里葉變換嚴(yán)謹(jǐn)?shù)臄?shù)學(xué)推導(dǎo)及其代碼實現(xiàn)過程,真正理解信號頻譜的概念及其在信號分析和處理中的應(yīng)用,激發(fā)了學(xué)生的學(xué)習(xí)興趣,提升了學(xué)生的理論分析能力和工程應(yīng)用開發(fā)能力。