周期信號雙邊頻譜圖的Matlab實現

王波

摘要：選擇一個特定的周期鋸齒波信號，求出其離散頻譜函數及其幅度頻譜和相位頻譜，編制程序繪制其雙邊頻譜圖，探討利用Matlab輔助教學手段，使周期信號復傅里葉級數及相關理論教學過程直觀、可視化的一般方法。

關鍵詞：周期信號;傅里葉級數;離散頻譜函數;幅度頻譜;相位頻譜;Matlab

中圖分類號：G642? ? ? ? 文獻標識碼：A

文章編號：1009-3044（2021）33-0118-02

開放科學（資源服務）標識碼（OSID）：

1 引言

非正弦周期信號分解為傅里葉級數的理論是具有驚人洞察力的法國科學家傅里葉（J.Fourier）在1807年提出，并于1822年公之于世的[1]。其后，泊松（Poisson）、高斯（Gauss）等人將這一成果應用到電學中，經過近兩個世紀時間的發(fā)展，如今傅里葉分析法已廣泛成熟地應用于眾多的科學技術領域[2]。

在諸如信號與系統(tǒng)、數字信號處理等眾多的高等學校理工科專業(yè)基礎課的課程體系中，傅里葉級數理論是系統(tǒng)分析由時域轉向頻域的重要關節(jié)點，起著承前啟后的重要作用。在依托傳統(tǒng)教學手段的基礎上，借助Matlab軟件使傅里葉級數及相關理論脫離抽象，可在一定程度上降低教學難度，取得更好的教學效果。以下，將以一個特定的周期鋸齒波信號的雙邊頻譜圖的Matlab實現過程為例對此進行說明。

2周期信號的復傅里葉級數及相關理論

2.1復傅里葉級數

狄里赫利（P.Dirichlet）嚴格證明了傅里葉級數理論，當周期信號滿足狄里赫利條件時，可將其展開成傅里葉級數，傅里葉級數分為三角形式傅里葉級數和復指數形式傅里葉級數兩種[3]。三角形式的傅里葉級數表明周期信號由直流成分及眾多的正弦、余弦函數形式的諧波分量疊加組成，而復指數級數形式的傅里葉級數簡稱為復傅里葉級數，其理論的核心內容體現在以下一對數學變換式中[4]。

[ft=n=-∞∞Fnejnω1t? ? ? ? ? ? ? ? ? ? ? ? ? ?（1）Fn=1T-T2T2f（t）e-jnω1tdt? ? ? ? ? ? ? ? ?（2）]

式（1）給出周期為T的連續(xù)時間信號f（t）的復傅里葉級數表示，其物理意義非常明確，意味著f（t）可分解為直流及無窮多個指數形式的諧波分量[Fnejnω1t]，其中的[ω1=2πT]稱為基波角頻率，而[nω1]則為f（t）中的n次諧波角頻率。

2.2離散頻譜函數及周期信號的雙邊頻譜

式（1）中各個諧波分量的復傅里葉系數[Fn]是實現周期信號復傅里葉級數表示的關鍵，可依據式（2）所示的積分變換求出[Fn]。

一般情況下，[Fn]為以各次諧波角頻率[nω1]為自變量的復函數，其極坐標形式為[Fn=Fnejφn]，[Fn]的模[|Fn|]稱為幅度頻譜，[Fn]的輻角[φn]則稱為相位頻譜，它們分別反映了組成周期信號的各個諧波分量的幅度和相位隨[nω1]的變化關系。因為整數n的取值范圍為[-∞～+∞]，所以幅度頻譜和相位頻譜均為雙邊頻譜，在對應的頻譜圖上，譜線都是僅出現在基波角頻率[ω1]整數倍即[nω1]處，故[Fn]亦稱為信號f（t）的離散頻譜函數。

3周期信號雙邊頻譜圖的實現

3.1實例時間信號

選擇一個周期T=2秒，峰值為1且位于每個周期最左側的周期鋸齒波信號f（t），在以t=0為奇函數對稱中心的周期內，其可以用函數[ft=-2Tt=-t]加以描述，以此實例信號探討利用Matlab實現周期信號雙邊頻譜圖的方法。

3.2實例信號的頻譜

實例信號f（t）的基波角頻率為[ω1=2π2=π（rad/s）]，利用式（2）求對應的離散頻譜函數。

[Fn=1T-T2T2f（t）e-jnω1tdt=12-11（-t）e-jnπtdt]

可知[F0=0]，即f（t）中不存在直流分量，而對于f（t）中的各交流諧波分量，則有

[Fn=-j1nπcosnπ=j（-1）n+11nπ（n=±1，±2，±3…）]

其中的幅度頻譜為

[Fn=1nπ=1nπ? ? ? ?（n=+1，+2，+3，…）-1nπ? ?（n=-1，-2，-3，…）? ? ?（3）]

相位頻譜為[φn=π2? ? ?（n=+1，+3，+5，…）-π2? ?（n=-1，-3，-5，…）-π2? ?（n=+2，+4，+6，…）π2? ? ?（n=-2，-4，-6，…）? ? ? （4）]

3.3雙邊頻譜圖的實現

根據上述討論內容編制本文后續(xù)3.4部分給出的Matlab程序，通過其中代碼的順序執(zhí)行，可完成以下4步操作。

（1）以交互方式輸入一個數值，后續(xù)代碼將依據此值設定頻譜圖中繪制的最高次諧波譜線的諧波次數;

（2）依次以[pi*t-1、0]為兩個實參調用sawtooth函數，生成并存儲實例周期鋸齒波信號f（t）的相關數據;

（3）根據式（3）、式（4）分別計算并存儲幅度頻譜[|Fn|]和相位頻譜[φn]的相關數據;

（4）在圖形窗口中，利用從上到下排列的（a）、（b）和（c）三個子窗口依次繪制f（t）的波形及其對應的雙邊幅度頻譜圖和雙邊相位頻譜圖。

設定最高次諧波次數為15時，運行程序將得到如圖1所示的圖形窗口，周期信號的復傅里葉級數分解的思想得以形象、直觀地圖示，復雜理論的教學過程變得相對簡單、方便。

可以借助圖1說明周期信號幅度頻譜的諧波性（譜線只能在基波角頻率[ω1]的整數倍即諧波角頻率[nω1]處出現）、離散性（頻譜圖由頻率離散的譜線組成，每根譜線代表一個諧波分量，兩根相鄰譜線之間的角頻率間隔為[ω1]）和收斂性（各譜線的高度即各諧波分量的幅度隨著諧波次數增高而逐漸減小并直至趨于無窮?。?。另外，幅度頻譜[|Fn|]是諧波角頻率[nω1]的偶函數，相位頻譜[φn]則是[nω1]的奇函數的特征在圖1中也是顯而易見的。

3.4 程序代碼

display（'請輸入最高次諧波次數Nmax'）;

Nmax=input（'Nmax= '）;

t=-5：0.0001：5;

f=sawtooth（pi*（t-1），0）;

AbsFn=zeros（2*Nmax+1，1）;Phase=AbsFn;

fori=（Nmax+2）：（2*Nmax+1）

AbsFn（i）=1/（i-（Nmax+1））/pi;

AbsFn（2*Nmax+2-i）=AbsFn（i）;

end

fori=（Nmax+2）：2：（2*Nmax+1）

Phase（i）=pi/2;

Phase（2*（Nmax+1）-i）=-pi/2;

end

fori=（Nmax+3）：2：（2*Nmax+1）

Phase（i）=-pi/2;

Phase（2*（Nmax+1）-i）=pi/2;

end

subplot（311）;plot（t，f，'linewidth'，1.5）;

axis（[-4.8 4.8 -1.1 1.1]）;

xlabel（'＼itt? ? ? ＼rm單位：s'，'Fontsize'，9）;

ylabel（'（a）? 周期鋸齒信號 f（t）'，'Fontsize'，9）;

k=-Nmax：Nmax;

subplot（312）;stem（k，AbsFn，'filled'）;

axis（[-NmaxNmax 0 0.34]）;

xlabel（'n＼omega_1? ? ? 單位：rad/s'，'Fontsize'，9）;

ylabel（'（b）? 幅度頻譜 |F_n|'，'Fontsize'，9）;

subplot（313）;stem（k，Phase，'filled'）;

axis（[-NmaxNmax -1.8 1.8]）;

xlabel（'n＼omega_1? ? ? 單位：rad/s'，'Fontsize'，9）;

ylabel（'（c）? 相位頻譜 ＼phi_n'，'Fontsize'，9）;

4結束語

不限于實例鋸齒波信號，只要按照上述說明的方法、步驟，將信號的函數表達式代入式（2）求出其對應的離散頻譜函數[Fn]，并根據[Fn]的具體形式計算出[|Fn|]和[φn]的相關數據，對Matlab程序中的相關代碼加以有針對性地簡單修改，此程序即可以應用于諸如周期性三角波、周期性方波等許多一般周期信號的雙邊頻譜圖的繪制。

本文所敘述的工作，僅僅是Matlab實際應用的一個簡單的嘗試探索。由于在數值計算處理、數據可視化分析和系統(tǒng)直觀描述等方面具有強大的能力，Matlab一定會在高等學校眾多課程的理論和實踐教學過程中得到更廣泛的應用。選擇恰當的知識點，適時使用Matlab輔助教學手段，一定會取得較之傳統(tǒng)教學手段更好的教學效果。

參考文獻：

[1] 燕慶明.信號與系統(tǒng)教程[M].北京：高等教育出版社，2013：86.

[2] 馬金龍.信號與系統(tǒng)[M].北京：科學出版社，2010：114.

[3] 王明泉.信號與系統(tǒng)[M].北京：科學出版社，2008：79.

[4] 沈元隆.信號與系統(tǒng)[M].北京：人民郵電出版社，2003：71.

【通聯(lián)編輯：王力】