傅立葉變換的實際應用解析

如果你對信號處理感興趣，說一天征服傅里葉有點夸張。當然，沒有反覆實踐和鉆研數(shù)學，您無法在一天里學會傅立葉變換的方方面面。無論如何, 這個在線課程將提供給您怎樣進行傅立葉變換運算的基本知識。

能有效和能非常簡單地領會的原因是我們使用了一種不太傳統(tǒng)的逼近。重要的是你將學習傅立葉變換的要素而完全不用超過加法和乘法的數(shù)學計算! 我將設法在不超過以下六節(jié)里解釋在對音像信號處理中傅立葉變換的實際應用。

步驟 1: 一些簡單的前提

在下面，您需要理解以下四件最基本的事情: 加法，乘、除法。什么是正弦，余弦和正弦信號。明顯地，我將跳第一二件事和將解釋位最后一個。您大概還記得您在學校學過的“ 三角函數(shù)”[1] ，它神秘地用于 與角度一起從它們的內角計算它們的邊長，反之亦然。我們這里不需要所有這些事，我們只需要知道二個 最重要的三角函數(shù)，"正弦" 和"余弦" 的外表特征。這相當簡單: 他們看起來象是以峰頂和谷組成的從觀察點向左右無限伸展的非常簡單的波浪。

（附圖一）

如同你所知道的，這兩種波形是周期性的，這意味著在一定的時間、周期之后，它們看起來再次一樣。兩種波形看起來也很象，但當正弦波在零點開始時余弦波開始出現(xiàn)在最大值。在實踐中，我們如何判定我們在一個給定時間所觀測到的波形是開始在它的最大值或在零? 問的好: 我們不能。實踐上沒有辦法區(qū)分正弦波和余弦波，因此看起來象正弦或余弦波的我們統(tǒng)稱為正弦波，在希臘語中譯作"正弦類"。正弦波的一個重要性質是"頻率"。它告訴我們在一個給定的時間內有多少個波峰和波谷。高頻意味許多波峰和波谷， 低頻率意味少量波峰和波谷:

（附圖二）

步驟 2: 了解傅立葉定理

Jean-Baptiste Joseph Fourier 是孩子們中讓父母感到驕傲和慚愧的的一個，因為他十四歲時就開始對他們說非常復雜的數(shù)學用語。他的一生中做了很多重要工作，但最重大的發(fā)現(xiàn)可能是解決了材料熱傳導問題。

他推導出了描述熱在某一媒介中如何傳導的公式，即用三角函數(shù)的無窮級數(shù)來解決這個問題（就是我們在上面討論過的正弦、余弦函數(shù)）。主要和我們話題有關的是：傅里葉的發(fā)現(xiàn)總結成一般規(guī)律就是任意復雜的信號都能由一個個混合在一起的正弦函數(shù)的和來表示。

這是一個例子：

（附圖三）

在這里你看到的是一個原始的信號，以及如何按某一確定的關系（“ 配方” ）混合在一起的正弦函數(shù)混合物（我們稱它們?yōu)榉至浚┧平?。我們將簡略地談論一下那份配方。如你所知，我們用的正弦函?shù)愈多 其結果就愈精確地接近我們的原始信號波形。在“ 現(xiàn)實” 世界中，在信號連續(xù)的地方，即你能以無窮小的間 隔來測量它們，精度僅受你的測試設備限制，你需要無限多的正弦函數(shù)才能完美地建立任意一個給定的信 號。幸運地是，和數(shù)字信號處理者們一樣，我們不是生活在那樣的世界。相反，我們將處理僅以有限精度 每隔一定間隔被測量的現(xiàn)實世界的采樣信號。因而，我們不需要無限多地正弦函數(shù)，我們只需要非常多。稍后我們也將討論這個“ 非常多” 是多少。目前重要的一點是你能夠想象，任意一個在你計算機上的信號， 都能用簡單正弦波按配方組成。

步驟 3: “ 非常多” 是多少

正如我們所知道的，復雜形狀的波形能由混合在一起的正弦波所建立。我們也許要問需要多少正弦波來 構造任意一個在計算機上給定的信號。當然，倘若我們知道正在處理的信號是如何組成的，這可能至少是 一個單個正弦波。在許多情況下，我們處理的現(xiàn)實世界的信號可能有非常復雜的結構，以至于我們不能深 入知道實際上有多少“ 分量” 波存在。在這種情況下，即使我們無法知道原始的信號是由多少個正弦波來構 成的，肯定存在一個我們將需要多少正弦波的上限。盡管如此，這實際上沒解決有多少的問題。讓我們試 著來直觀地逼近它: 假設一個信號我們有 1000 個樣采，可能存在的最短周期正弦波（即多數(shù)波峰波谷在其 中）以交替的波峰波谷分布在每個采樣內。因此，最高頻率的正弦波將有 500 個波峰和 500 個波谷在我 們的 1000 個采樣中，且每隔一個采樣是波峰。下圖中的黑點表示我們的采樣，所以，最高頻率的正弦波 以看起來象這樣:

The highest frequency sine wave（附圖四）

現(xiàn)在讓我們來看一下最低頻率正弦波可能多么低。如果我們只給一個單獨的采樣點，我們將如何能測量穿過這點的正弦波的峰頂和谷? 我們做不到，因為有許多不同周期正弦波穿過這點。

（附圖五）

所以，一個單獨數(shù)據(jù)點不足以告訴我們關于頻率的任何事。現(xiàn)在，如果我們有兩個采樣，那么穿過這 兩點的正弦波的最低頻率是什么？在這種情況下它很常簡單。只有一個穿過這兩點的非常低頻率的正弦 波。它看起來向這樣：

The lowest frequency sine wave（附圖六）

想象最左面的兩個點是二個釘子和一個跨越它們之間的弦（圖六描述三個數(shù)據(jù)點是因為正弦波的周期性，但我們實際上只需要最左面的兩點說明它的頻率）。我們能領會的最低的頻率是來回地擺動在二個釘 子之間的弦，象是圖六中左邊兩點之間的正弦波所做的。如果我們有 1000 個采樣，那么兩個“ 釘子” 相當 于第一個或最后的采樣，比如 1 號采樣和 1000 號采樣。從對樂器的體驗我們知道，當長度增加時弦的頻 率將下降。所以我們可以想象，當我們將兩個釘子向彼此遠離的方向移開時，最小正弦波的頻率將變得更 低。例如，如果我們選擇 2000 個采樣，因為我們的“ 釘子” ，在 1 號或 2000 號采樣間，所以最低正弦波將 更低。事實上，它將低兩倍，因為因為我們的釘子比 1000 個采樣時遠兩倍。這樣，如果我們有更多的采 樣，我們將能辨別出一個更低頻率的正弦波，因為它們的零交叉點（我們的“ 釘子” ）將移動得更遠。這對 了解下面的解釋是非常重要的。

象我們看到的那樣，在兩個“ 釘子” 之后，我們的波形開始重復上升斜坡（第一個釘子和第三個釘子同樣）。這意味著任意兩個相鄰的釘子準確地包含完整正弦波的一半，換句話說一個峰或一個谷或半個周期。

概括一下我們剛學過的東西，我們知道，一個采樣正弦波的上限頻率是所有其它有一個波峰和波谷的 采樣，并且，低頻下限是我們看到的正好匹配采樣數(shù)的正弦波的周期的一半。但等一下，這難到不意味著， 當上限頻率保持固定時，當有很多采樣時最低頻率可以降低？確實如此! 結果我們將從一個低頻開始增加 更多正弦波來組成一個較長的未知內容的信號。

一切都清楚了，但我們仍然不知道我們最終需要多少正弦波。就象我們現(xiàn)在知道任意正弦波分量所能 具的有上下限頻率一樣，我們能計算出適合這兩個極限之間的量是多少。既然我們沿著最左到最右的采樣 固定了最低正弦波分量，我們要所有其它正弦波最好也使用這些釘子(為什么我們要不同地對待他們? 所有的正弦波被同等的創(chuàng)建！)。假設正弦波束是系在吉他上二個固定點的弦。它們能只搖擺在二個釘子之間(除非他們斷了)，就象我們下圖的正弦波。這導致如下關系，我們的最低分量（1）以 1/2 周期裝配，第二分量（2）以 1 個周期裝配，第三分量以 1 1/2 周期裝配以此類推直到我們看到的 1000 個采樣。形象地, 看起來象這:

（附圖七）

在，如果我們算一下要多少正弦波以那種方法裝配我們的 1000 個采樣，就會發(fā)現(xiàn)我們精確地需要1000 個正弦波疊加起來表示 1000 個采樣。實際上，我們總是發(fā)現(xiàn)我們需要和采樣一樣多的正弦波。

步驟 4 關于烹飪食譜

在前面的段落我們看到, 任一個給定的在計算機上的信號能被正弦波混合物來構造。我們考慮了他們的頻率，并且考慮了需要多大的最低和最高頻率的正弦波來完美地重建任一個我們所分析信號。我們明白 了為確定所需最低的正弦波分量，我們考察的采樣的數(shù)量是重要，但我們還未論述實際正弦波如何必需被 混合產(chǎn)生某一確定的結果。由疊加正弦波組成任何指定的信號，我們需要測量他們的另外一個方面。實際 上，頻率不是我們需要知道的唯一的事。我們還需要知道正弦波的幅度，也就是說每個正弦波幅度有多高 才能混合在一起產(chǎn)生我們需要的輸入信號。高度是正弦波的峰頂?shù)母叨?，意即峰頂和零線之間距離。幅度 赿高，我們聽到的聲音也就赿大。所以，如果您有一個含有許多低音的信號，無疑可以預期混合體中的低 頻率正弦波的分量比例比高頻正弦波分量更大。因此一般情況下，低音中的低頻正弦波有一個比高頻正弦 波更高的幅度。在我們的分析中，我們將需要確定各個分量正弦波的幅度以完成我們的配方。

步驟 5：關于蘋果和桔子

如果你一直跟著我，我們幾乎完成了通向傅里葉變換的旅程。我們學了需要多少正弦波，它的數(shù)量依 賴于我們查看的采樣的數(shù)量有一個頻率上下限界，并且不知道怎么確定單個分量的幅度以完成我們的配 方。我們一直不清楚究竟如何從我們的采樣來確定實際的配方。直觀上我們可以斷定能找到正弦波的幅度， 設法把一個已知頻率正弦波和采樣作對比，我們測量找出它們有多么接近。如果它們精確地相等，我們知 道該正弦波存在著相同的幅度，如果我們發(fā)現(xiàn)我們的信號與參考正弦波一點也不匹配，我們將認為這個不 存在。盡管如此，我們如何高效地把一個已知的正弦波同采樣信號進行比較？幸運地是，數(shù)字信號處理工 作者早已解決了如何作這些。事實上，這象加法和乘法一樣容易－我們取一個已知頻率的單位正弦波（這意味著它的振幅是 1，可從我們的計算器或計算機中精確地獲得）和我們的信號采樣相乘。累加乘積之后， 我們將得到我們正在觀測的這個頻率上正弦波分量的幅度。

這是個舉例，一個簡單的完成這些工作的 C 代碼片段：

Listing 1.1: The direct realization of the Discrete Sine Transform (DST):

#define M_PI 3.14159265358979323846

long bin,k;

double arg;

for (bin = 0; bin < transformLength; bin++) {

transformData[bin] = 0.;

for (k = 0; k < transformLength; k++) {

arg = (float)bin * M_PI *(float)k / (float)transformLength;

transformData[bin] += inputData[k] * sin(arg);
}
}

這段代碼變換存儲在 inputData[0...transformLength-1]中我們測量的采樣點成為一個正弦波分量的幅度隊列 transformData[0...transformLength-1]。根據(jù)通用術語，我們稱參考正弦波的頻率步長為盒（bin），

這意味著它們被認為象是一個我們放置我們估計的任意分量波的幅度的容器。離散正弦變換(DST)是一個普通程序，它假設我們無法想象我們的信號看起來象什么樣，否則我們能使用一個更加高效率的方法來確 定正弦波分量的幅度（例如，我們予先知道, 我們的信號是一個已知頻率的正弦波。我們能直接地找出它 的高度而不用計算正弦波的整個范圍。實現(xiàn)這個有效的逼近是基于傅里葉原理，它能在文獻的戈策爾(Goe rtzel)算法條目下找到）。

這些就是你堅持想要的我們?yōu)槭裁从眠@樣的方法計算正弦變換的一個解釋：對我們用一個已知頻率正弦 波的乘積來作一種非常直觀逼近的理由，可以設想，這大致相當于一個固有頻率的“ 共振” 在系統(tǒng)內發(fā)生時 物理世界發(fā)生的事情。sin(arg)項本質上是一個獲得由輸入信號波形激勵的諧振器。如果輸入（信號）有在 我們正觀測的頻率上的分量，它的輸出將是參考正弦波諧振的幅度。因為我們的參考波是單位幅度的，輸出是一個在那個頻率上的分量的實際幅度的一個直接測量。因為諧振器只是簡單的濾波器，變換（不可否 認是在稍微寬松條件下）被認為有極窄的帶通濾波器組的特征，它位于我們估值的頻率中心的周圍。這有 助于解釋一個事實，為什么傅立葉變換提供了對信號進行過濾的一個高效工具。

只是為了完備性: 當然，上述程序是可逆的，當我們知道它的正弦波分量時，我們的信號(在數(shù)字精確度 極限內)能完全被重建，通過簡單地把正弦波加起來。這留下給讀者做為一個練習。同樣程序能改變使用余 弦波做為基本函數(shù)工作-我們只需簡單地改變 sin(arg)條件到 cos(arg)來獲得離散余弦變換的直接實現(xiàn)(DC T)。

現(xiàn)在，就象在這篇文章較前面的段落中我們討論過的那樣，我們在實踐中沒有辦法區(qū)分一個被測量的正 弦類函數(shù)象是正弦波還是余弦波。做為代替我們總是測量正弦信號，且正弦和余弦變換在實踐中沒有太大 的用途，除了一些特殊情況（象圖象壓縮的地方，即毎塊圖象具有能用一個基本的余弦或正弦函數(shù)較好模 擬特性，例如能用余弦基本函數(shù)較好表現(xiàn)的相同顏色的大區(qū)域）。正弦信號是一個比正弦或余弦波更一般 的片斷，因為它可以開始在在它的周期中的一個任意位置。我們記得，當余弦波開始于 1 時，正弦波總開 始在 0。當我們采取正弦波作為參考，余弦波開始在它的周期的最后 1/4 之處。一般用度或弧度測量它們 的偏移量，這是兩個一般與三角函數(shù)相關的單位。一個完整的周期等于 360°（代表 "度") 或 2π 個弧度(代表"2π" ，"π" 發(fā)音象“pie” 。π 是希臘字表示數(shù) 3.14159265358979323846... 在三角學方面有重要意義) 。

余弦波因而有一個 90°或 π/2 的偏移。這偏移叫正弦信號的相位，因此余弦波相對正弦信號有 90°或 π/2 相位。

相位的事情就有這些內容。因為我們一直不能限定信號在 0°或 90°相位開始（因為我們正觀測一個我們 可能無法控制的信號），它對同時直接唯一的描述信號的頻率、振幅、相位至關重要。以正弦或余弦做變換相位限制在 0°或 90°，一個具有任意相位的正弦信號將引起相鄰頻率出現(xiàn)假峰（因為它們試圖“ 幫助” 分析，強制給被測信號加上一個 0°或 90°的相位作用）。它有些象用一圓石頭去填滿一個方孔：你需要小一 些的圓石頭去填充剩余的空間，并且更小的石頭填好依然留出空的空間，等等。我們所需要的是能處理一 般信號的變換，它能處理任意相位正弦波構成的信號。

步驟 6：離散傅葉變換

從正弦變換到傅里葉變換的步驟是簡單的，只需用更一般的方法。在正弦變換中對每個頻率上的測度 使用正弦波，在傅里葉變換中正弦、余弦波二者都使用。就是說，對任意的當前頻率，我們以同一頻率的 正弦和余弦波來“ 比較” （或“ 共振” ）被測信號。如果我們的信號看起來很象正弦波，變換的正弦部份將有 一個大的幅值。如果它看起來象余弦波，變換的余弦部份將有一個大的幅值。如果看起來象反相的正弦波， 也就是說，它開始于 0 但下降至－1 取代上升至 1，它的正弦部份將有一個大的負幅值。這表明用＋、－ 符號和正弦、余弦相位能表示任意給定頻率的正弦信號[2]。

Listing 1.2: The direct realization of the Discrete Fourier Transform[3]:

#define M_PI 3.14159265358979323846

long bin, k;

double arg, sign = -1.; /* sign = -1 -> FFT, 1 -> iFFT */

for (bin = 0; bin <= transformLength/2; bin++) {

cosPart[bin] = (sinPart[bin] = 0.);

for (k = 0; k < transformLength; k++) {

arg = 2.*(float)bin*M_PI*(float)k/(float)transformLength;

sinPart[bin] += inputData[k] * sign * sin(arg);

cosPart[bin] += inputData[k] * cos(arg);
}
}

我們仍遺留著一個問題，就是如何獲得傅里葉變換所缺乏的那些有用的東西。我說過傳里葉變換的優(yōu)赿性超過正弦和余弦變換是因為用正弦信號工作。但至今我們還未看到任何正弦信號，仍只有正弦和余弦。好，這需要一點附加處理步聚：

#define M_PI 3.14159265358979323846

long bin;

for (bin = 0; bin <= transformLength/2; bin++) {

/* frequency */

frequency[bin] = (float)bin * sampleRate / (float)transformLength;

/* magnitude */

magnitude[bin] = 20. * log10( 2. * sqrt(sinPart[bin]*sinPart[bin] + cosPart[bin]*c

osPart[bin]) / (float)transformLength);

/* phase */

phase[bin] = 180.*atan2(sinPart[bin], cosPart[bin]) / M_PI - 90.;

}

在運行清單 1.3 所示的關于 DFT 輸出的代碼段之后，我們結束被看作以正弦信號波的和的輸入信號表示。

K 序正弦信號是用 frequency[k],magnitude[k]和 phase[k]來描述的。單位是 Hz(Hertz,周/秒)，dB(Decibel),和°(Degree)。請注意在經(jīng)過清單 1.3 的后加工（處理）即把正弦和余弦函數(shù)部份轉換成一個單一的正弦信 號之后，我們命名 K 序正弦信號的振幅－DFT 存貯為幅度，且它總是取相對值。我們可以說一個-1.0 的振 幅對應于 1.0 的幅度，對應于相位+或-180°。在文獻中，做做傅里葉變換的場合，隊列 magnitude[] 被稱 作被測信號的幅度譜，隊列 phase[]被稱作被測信號的相位譜。

如用分貝測量存貯幅度的參考，輸入波也期望有一個在[-1.0,1.0]之間的采樣值，相對于 0dB 幅度滿刻 度數(shù)字。做為一個 DFT 的有趣應用，比如清單 1.3 就可被用于寫一個基于離散傅里葉變換的譜分析。

結論

象我們已知那樣，傅里葉變換和其系列的離散正弦和余弦變換，提供了把一個信號分解成一束分波的 便利工具。結果有正弦或余弦之一，或正弦信號（用正弦和余弦波的組合來描述）。在傅里葉變換中同時 使用正弦和余弦波的好處是我們因而能引入相位的概念，它使變換更一般化，因而我們能用它有效清楚地 分析既不是純正弦也不是純余弦的正弦信號，當然其它信號也一樣。

傅里葉變換與被考察信號無關，因而無論我們正分析的信號是一個正弦信號或是一些其它的更復雜的， 變換需要相同的操作數(shù)。這就是為什么傅里葉變換被稱做無參數(shù)變換的原因，這意味著它對需要的信號“ 智能的” 分析沒有直接的幫助（在考察一個我們已知是一個信號是正弦曲線的情況下，我們更喜歡精確地獲 得關于相位，頻率，幅度的信息以代替一串在一些預定頻率上的正弦和余弦波）。

現(xiàn)在我們也知道了我們是在求輸入信號在一組固定頻率柵格上的值，輸入信號實際存在的頻率組在這 組柵格上可能不起作用。我們在分析中利用的柵格是人為的，因為我們幾乎按照關于它們的頻率的嘗試來 選擇參考正弦、余弦波。說到了這些，馬上清楚了一個將要很容易遇到的要點，即被測信號的頻率位于變 換柵格的頻率之間。因此，有一個頻率發(fā)生在位于兩個頻率柵格之間的正弦曲線，在變換中將不好被描述。包圍著與輸入信號頻率最接近的柵格的相鄰的柵格將試圖‘改正’ 頻率的背離。因而，輸入信號的能量將拖 尾至數(shù)個相鄰的柵格。這也是傅里葉變換不能迅速地分析聲音返回它的基波和諧波（并且，這也是為什么 我們稱正弦和余弦波為分波而不諧波和泛音）。

簡單的說，沒有進一步的快速處理，DFT 和一個狹窄的壩一樣，細小并行的帶通濾波器組（“ 通道” ）和每個通道帶有附加的相位信息。這對分析信號、做濾波器和運用其它的技巧是有益的（改變一個信號的 音調而不改變它的速度是它們其中之一，說明在 DSPdimension.com 上另一篇不同的文章中），但它需 要對少量普通任務附加快速處理。同樣，它能被認為是使用除了正弦和余弦波基本函數(shù)的變換系列的一個 特例。在這個方向上展開概念超出了這篇文章的范圍。

最后，重要的是要提及一個更高效的 DFT 工具，也就是一個被稱做快速傅里葉變換的算法。它最初 是由庫利和圖克在 1969 年構思的（它的根源仍然要追溯到高斯和其它人的工作）。FFT 只是一個高效的 算法，它比上面給出的以直接逼近計算 DFT 所化的時間少，它是結果完全相同的其它方法。無論如何，F(xiàn) FT 是以庫利/圖克算法實施的，它需要變換長度是 2 的冪。在實踐中，對大多數(shù)應用來說這是一個可以接 受的限制。有大量的以不同方法實施 FFT 的可利用的文獻，因而，可以說有足夠多不同的 FFT 實現(xiàn)，其中一些并不需要經(jīng)典 FFT 的 2 的冪的限制。

下面清單 1.4 以程序 smbFft()給出了一個 FFT 的實現(xiàn)。

Listing 1.4: The Discrete Fast Fourier Transform (FFT):

#define M_PI 3.14159265358979323846

void smbFft(float *fftBuffer, long fftFrameSize, long sign)

/*

FFT routine, (C)1996 S.M.Bernsee. Sign = -1 is FFT, 1 is iFFT (inverse)

Fills fftBuffer[0...2*fftFrameSize-1] with the Fourier transform of the time domain data in fftBuffer

[0...2*fftFrameSize-1]. The FFT array takes and returns the cosine and sine parts in an interleaved

manner, ie. fftBuffer[0] = cosPart[0], fftBuffer[1] = sinPart[0], asf. fftFrameSize must be a power of

2. It expects a complex input signal (see footnote 2), ie. when working with ＇common＇ audio sig

nals


our input signal


has to be passed as {in[0],0.,in[1],0.,in[2],0.,...} asf. In that case, the transform

of


the frequencies of


interest is in fftBuffer[0...fftFrameSize].

*/

{

float wr, wi, arg, *p1, *p2, temp;

float tr, ti, ur, ui, ur, ui, *p1r, *p1i, *p2r, *p2i;

long i, bitm, j, le, le2, k;

for (i = 2; i < 2*fftFrameSize-2; i += 2) {

for (bitm = 2, j = 0; bitm < 2*fftFrameSize; bitm <<= 1) {

if (i & bitm) j++;

j <<= 1;

}

if (i < j) {


p1 = fftBuffer+i; p2 = fftBuffer+j;

temp = *p1; *(p1++) = *p2;

*(p2++) = temp; temp = *p1;

*p1 = *p2; *p2 = temp;

}

}

for (k = 0, le = 2; k < (long)(log(fftFrameSize)/log(2.)); k++) {

le <<= 1;

le2 = le>>1;

ur = 1.0;

ui = 0.0;

arg = M_PI / (le2>>1);

wr = cos(arg);

wi = sign*sin(arg);

for (j = 0; j < le2; j += 2) {

p1r = fftBuffer+j; p1i = p1r+1;

p2r = p1r+le2; p2i = p2r+1;

for (i = j; i < 2*fftFrameSize; i += le) {

tr = *p2r * ur - *p2i * ui;

ti = *p2r * ui + *p2i * ur;

*p2r = *p1r - tr; *p2i = *p1i - ti;

*p1r += tr; *p1i += ti;

p1r += le; p1i += le;

p2r += le; p2i += le;

}


tr = ur*wr - ui*wi;

ui = ur*wi + ui*wr;

ur = tr;