|
經(jīng)典算法研究系列：十、從頭到尾徹底理解傅里葉變換算法、上 作者：July、dznlong 推薦閱讀：The Scientist and Engineer's Guide to Digital Signal Processing，By Steven W. Smith, Ph.D。此書地址：http://www./pdfbook.htm。 博主說明：I、本文中闡述離散傅里葉變換方法，是根據(jù)此書：The Scientist and Engineer's Guide to Digital Signal Processing，By Steven W. Smith, Ph.D.而翻譯而成的，此書地址：http://www./pdfbook.htm。II、同時，有相當(dāng)一部分內(nèi)容編輯整理自dznlong的博客，也貼出其博客地址，向原創(chuàng)的作者表示致敬：http://blog.csdn.net/dznlong 。這年頭，真正靜下心寫來原創(chuàng)文章的人，很少了。 從頭到尾徹底理解傅里葉變換算法、上 前言 第一部分、 DFT 第一章、傅立葉變換的由來 第二章、實數(shù)形式離散傅立葉變換（Real DFT） 從頭到尾徹底理解傅里葉變換算法、下 第三章、復(fù)數(shù) 第四章、復(fù)數(shù)形式離散傅立葉變換 前言： “關(guān)于傅立葉變換，無論是書本還是在網(wǎng)上可以很容易找到關(guān)于傅立葉變換的描述，但是大都是些故弄玄虛的文章，太過抽象，盡是一些讓人看了就望而生畏的公式的羅列，讓人很難能夠從感性上得到理解”---dznlong， 那么，到底什么是傅里葉變換算法列?傅里葉變換所涉及到的公式具體有多復(fù)雜列? 傅里葉變換（Fourier transform）是一種線性的積分變換。因其基本思想首先由法國學(xué)者傅里葉系統(tǒng)地提出，所以以其名字來命名以示紀(jì)念。 哦，傅里葉變換原來就是一種變換而已，只是這種變換是從時間轉(zhuǎn)換為頻率的變化。這下，你就知道了，傅里葉就是一種變換，一種什么變換列?就是一種從時間到頻率的變化或其相互轉(zhuǎn)化。 ok，咱們再來總體了解下傅里葉變換，讓各位對其有個總體大概的印象，也順便看看傅里葉變換所涉及到的公式，究竟有多復(fù)雜： 以下就是傅里葉變換的4種變體（摘自，維基百科） 連續(xù)傅里葉變換 一般情況下，若“傅里葉變換”一詞不加任何限定語，則指的是“連續(xù)傅里葉變換”。連續(xù)傅里葉變換將平方可積的函數(shù)f（t）表示成復(fù)指數(shù)函數(shù)的積分或級數(shù)形式。  這是將頻率域的函數(shù)F(ω)表示為時間域的函數(shù)f（t）的積分形式。 連續(xù)傅里葉變換的逆變換 (inverse Fourier transform)為：  即將時間域的函數(shù)f（t）表示為頻率域的函數(shù)F(ω)的積分。 一般可稱函數(shù)f（t）為原函數(shù)，而稱函數(shù)F(ω)為傅里葉變換的像函數(shù)，原函數(shù)和像函數(shù)構(gòu)成一個傅里葉變換對（transform pair）。 除此之外，還有其它型式的變換對，以下兩種型式亦常被使用。在通信或是信號處理方面，常以  來代換，而形成新的變換對:  或者是因系數(shù)重分配而得到新的變換對：  一種對連續(xù)傅里葉變換的推廣稱為分?jǐn)?shù)傅里葉變換（Fractional Fourier Transform）。分?jǐn)?shù)傅里葉變換(fractional Fourier transform,FRFT)指的就是傅里葉變換(Fourier transform,FT)的廣義化。 分?jǐn)?shù)傅里葉變換的物理意義即做傅里葉變換 a 次，其中 a 不一定要為整數(shù)；而做了分?jǐn)?shù)傅里葉變換之后，信號或輸入函數(shù)便會出現(xiàn)在介于時域(time domain)與頻域(frequency domain)之間的分?jǐn)?shù)域(fractional domain)。 當(dāng)f（t）為偶函數(shù)（或奇函數(shù)）時，其正弦（或余弦）分量將消亡，而可以稱這時的變換為余弦變換（cosine transform）或正弦變換（sine transform）. 另一個值得注意的性質(zhì)是，當(dāng)f（t）為純實函數(shù)時，F(xiàn)(?ω) = F*(ω)成立. 傅里葉級數(shù) 連續(xù)形式的傅里葉變換其實是傅里葉級數(shù) (Fourier series)的推廣，因為積分其實是一種極限形式的求和算子而已。對于周期函數(shù)，其傅里葉級數(shù)是存在的：  其中Fn為復(fù)幅度。對于實值函數(shù)，函數(shù)的傅里葉級數(shù)可以寫成：  其中an和bn是實頻率分量的幅度。 離散時域傅里葉變換 離散傅里葉變換是離散時間傅里葉變換（DTFT）的特例（有時作為后者的近似）。DTFT在時域上離散，在頻域上則是周期的。DTFT可以被看作是傅里葉級數(shù)的逆變換。 離散傅里葉變換 離散傅里葉變換（DFT），是連續(xù)傅里葉變換在時域和頻域上都離散的形式，將時域信號的采樣變換為在離散時間傅里葉變換（DTFT）頻域的采樣。在形式上，變換兩端（時域和頻域上）的序列是有限長的，而實際上這兩組序列都應(yīng)當(dāng)被認(rèn)為是離散周期信號的主值序列。即使對有限長的離散信號作DFT，也應(yīng)當(dāng)將其看作經(jīng)過周期延拓成為周期信號再作變換。在實際應(yīng)用中通常采用快速傅里葉變換以高效計算DFT。 為了在科學(xué)計算和數(shù)字信號處理等領(lǐng)域使用計算機進(jìn)行傅里葉變換，必須將函數(shù)xn定義在離散點而非連續(xù)域內(nèi)，且須滿足有限性或周期性條件。這種情況下，使用離散傅里葉變換（DFT），將函數(shù)xn表示為下面的求和形式：  其中Xk是傅里葉幅度。直接使用這個公式計算的計算復(fù)雜度為O（n*n），而快速傅里葉變換（FFT）可以將復(fù)雜度改進(jìn)為O（n*lgn）。（后面會具體闡述FFT是如何將復(fù)雜度降為O（n*lgn）的。）計算復(fù)雜度的降低以及數(shù)字電路計算能力的發(fā)展使得DFT成為在信號處理領(lǐng)域十分實用且重要的方法。 下面，比較下上述傅立葉變換的4種變體，  如上，容易發(fā)現(xiàn)：函數(shù)在時（頻）域的離散對應(yīng)于其像函數(shù)在頻（時）域的周期性。反之連續(xù)則意味著在對應(yīng)域的信號的非周期性。也就是說，時間上的離散性對應(yīng)著頻率上的周期性。同時，注意，離散時間傅里葉變換，時間離散，頻率不離散，它在頻域依然是連續(xù)的。 如果，讀到此，你不甚明白，大沒關(guān)系，不必糾結(jié)于以上4種變體，繼續(xù)往下看，你自會豁然開朗。（有什么問題，也懇請?zhí)岢觯蛘吲u指正） ok， 本文，接下來，由傅里葉變換入手，后重點闡述離散傅里葉變換、快速傅里葉算法，到最后徹底實現(xiàn)FFT算法，全篇力求通俗易懂、閱讀順暢，教你從頭到尾徹底理解傅里葉變換算法。由于傅里葉變換，也稱傅立葉變換，下文所稱為傅立葉變換，同一個變換，不同叫法，讀者不必感到奇怪。 第一部分、DFT 第一章、傅立葉變換的由來 要理解傅立葉變換，先得知道傅立葉變換是怎么變換的，當(dāng)然，也需要一定的高等數(shù)學(xué)基礎(chǔ)，最基本的是級數(shù)變換，其中傅立葉級數(shù)變換是傅立葉變換的基礎(chǔ)公式。 一、傅立葉變換的提出 傅立葉是一位法國數(shù)學(xué)家和物理學(xué)家，原名是Jean Baptiste Joseph Fourier(1768-1830), Fourier于1807年在法國科學(xué)學(xué)會上發(fā)表了一篇論文，論文里描述運用正弦曲線來描述溫度分布，論文里有個在當(dāng)時具有爭議性的決斷：任何連續(xù)周期信號都可以由一組適當(dāng)?shù)恼仪€組合而成。 當(dāng)時審查這個論文拉格朗日堅決反對此論文的發(fā)表，而后在近50年的時間里，拉格朗日堅持認(rèn)為傅立葉的方法無法表示帶有棱角的信號，如在方波中出現(xiàn)非連續(xù)變化斜率。直到拉格朗日死后15年這個論文才被發(fā)表出來。 誰是對的呢？拉格朗日是對的：正弦曲線無法組合成一個帶有棱角的信號。但是，我們可以用正弦曲線來非常逼近地表示它，逼近到兩種表示方法不存在能量差別，基于此，傅立葉是對的。 為什么我們要用正弦曲線來代替原來的曲線呢？如我們也還可以用方波或三角波來代替呀，分解信號的方法是無窮多的，但分解信號的目的是為了更加簡單地處理原來的信號。 用正余弦來表示原信號會更加簡單，因為正余弦擁有原信號所不具有的性質(zhì)：正弦曲線保真度。一個正余弦曲線信號輸入后，輸出的仍是正余弦曲線，只有幅度和相位可能發(fā)生變化，但是頻率和波的形狀仍是一樣的。且只有正余弦曲線才擁有這樣的性質(zhì)，正因如此我們才不用方波或三角波來表示。 二、傅立葉變換分類 根據(jù)原信號的不同類型，我們可以把傅立葉變換分為四種類別： 1、非周期性連續(xù)信號 傅立葉變換（Fourier Transform） 2、周期性連續(xù)信號 傅立葉級數(shù)(Fourier Series) 3、非周期性離散信號 離散時域傅立葉變換（Discrete Time Fourier Transform） 4、周期性離散信號 離散傅立葉變換(Discrete Fourier Transform) 下圖是四種原信號圖例（從上到下，依次是FT，F(xiàn)S，DTFT，DFT）：  這四種傅立葉變換都是針對正無窮大和負(fù)無窮大的信號，即信號的的長度是無窮大的，我們知道這對于計算機處理來說是不可能的，那么有沒有針對長度有限的傅立葉變換呢？沒有。因為正余弦波被定義成從負(fù)無窮小到正無窮大，我們無法把一個長度無限的信號組合成長度有限的信號。 面對這種困難，方法是：把長度有限的信號表示成長度無限的信號。如，可以把信號無限地從左右進(jìn)行延伸，延伸的部分用零來表示，這樣，這個信號就可以被看成是非周期性離散信號，我們可以用到離散時域傅立葉變換（DTFT）的方法。也可以把信號用復(fù)制的方法進(jìn)行延伸，這樣信號就變成了周期性離散信號，這時我們就可以用離散傅立葉變換方法（DFT）進(jìn)行變換。本章我們要講的是離散信號，對于連續(xù)信號我們不作討論，因為計算機只能處理離散的數(shù)值信號，我們的最終目的是運用計算機來處理信號的。 但是對于非周期性的信號，我們需要用無窮多不同頻率的正弦曲線來表示，這對于計算機來說是不可能實現(xiàn)的。所以對于離散信號的變換只有離散傅立葉變換（DFT）才能被適用，對于計算機來說只有離散的和有限長度的數(shù)據(jù)才能被處理，對于其它的變換類型只有在數(shù)學(xué)演算中才能用到，在計算機面前我們只能用DFT方法，后面我們要理解的也正是DFT方法。 這里要理解的是我們使用周期性的信號目的是為了能夠用數(shù)學(xué)方法來解決問題，至于考慮周期性信號是從哪里得到或怎樣得到是無意義的。 每種傅立葉變換都分成實數(shù)和復(fù)數(shù)兩種方法，對于實數(shù)方法是最好理解的，但是復(fù)數(shù)方法就相對復(fù)雜許多了，需要懂得有關(guān)復(fù)數(shù)的理論知識，不過，如果理解了實數(shù)離散傅立葉變換(real DFT)，再去理解復(fù)數(shù)傅立葉變換就更容易了，所以我們先把復(fù)數(shù)的傅立葉變換放到一邊去，先來理解實數(shù)傅立葉變換，在后面我們會先講講關(guān)于復(fù)數(shù)的基本理論，然后在理解了實數(shù)傅立葉變換的基礎(chǔ)上再來理解復(fù)數(shù)傅立葉變換。 還有，這里我們所要說的變換(transform)雖然是數(shù)學(xué)意義上的變換，但跟函數(shù)變換是不同的，函數(shù)變換是符合一一映射準(zhǔn)則的，對于離散數(shù)字信號處理（DSP），有許多的變換：傅立葉變換、拉普拉斯變換、Z變換、希爾伯特變換、離散余弦變換等，這些都擴展了函數(shù)變換的定義，允許輸入和輸出有多種的值，簡單地說變換就是把一堆的數(shù)據(jù)變成另一堆的數(shù)據(jù)的方法。 三、一個關(guān)于實數(shù)離散傅立葉變換(Real DFT)的例子 先來看一個變換實例，下圖是一個原始信號圖像：  這個信號的長度是16，于是可以把這個信號分解9個余弦波和9個正弦波（一個長度為N的信號可以分解成N/2+1個正余弦信號，這是為什么呢？結(jié)合下面的18個正余弦圖,我想從計算機處理精度上就不難理解，一個長度為N的信號，最多只能有N/2+1個不同頻率，再多的頻率就超過了計算機所能所處理的精度范圍），如下圖： 9個余弦信號： 9個正弦信號： 把以上所有信號相加即可得到原始信號，至于是怎么分別變換出9種不同頻率信號的，我們先不急，先看看對于以上的變換結(jié)果，在程序中又是該怎么表示的，我們可以看看下面這個示例圖： 上圖中左邊表示時域中的信號，右邊是頻域信號表示方法， 從左向右，-->，表示正向轉(zhuǎn)換(Forward DFT)，從右向左，<--，表示逆向轉(zhuǎn)換(Inverse DFT)， 用小寫x[]表示信號在每個時間點上的幅度值數(shù)組, 用大寫X[]表示每種頻率的副度值數(shù)組（即時間x-->頻率X）, 因為有N/2+1種頻率，所以該數(shù)組長度為N/2+1， X[]數(shù)組又分兩種，一種是表示余弦波的不同頻率幅度值：Re X[]， 另一種是表示正弦波的不同頻率幅度值：Im X[]， Re是實數(shù)(Real)的意思，Im是虛數(shù)(Imagine)的意思，采用復(fù)數(shù)的表示方法把正余弦波組合起來進(jìn)行表示，但這里我們不考慮復(fù)數(shù)的其它作用，只記住是一種組合方法而已，目的是為了便于表達(dá)（在后面我們會知道，復(fù)數(shù)形式的傅立葉變換長度是N，而不是N/2+1）。如此，再回過頭去，看上面的正余弦各9種頻率的變化，相信，問題不大了。 第二章、實數(shù)形式離散傅立葉變換（Real DFT） 上一章，我們看到了一個實數(shù)形式離散傅立葉變換的例子，通過這個例子能夠讓我們先對傅立葉變換有一個較為形象的感性認(rèn)識，現(xiàn)在就讓我們來看看實數(shù)形式離散傅立葉變換的正向和逆向是怎么進(jìn)行變換的。在此，我們先來看一下頻率的多種表示方法。 一、 頻域中關(guān)于頻率的四種表示方法 1、序號表示方法，根據(jù)時域中信號的樣本數(shù)取0 ~ N/2，用這種方法在程序中使用起來可以更直接地取得每種頻率的幅度值，因為頻率值跟數(shù)組的序號是一一對應(yīng)的: X[k]，取值范圍是0 ~ N/2； 2、分?jǐn)?shù)表示方法，根據(jù)時域中信號的樣本數(shù)的比例值取0 ~ 0.5: X[?]，? = k/N，取值范圍是0 ~ 1/2； 3、用弧度值來表示，把?乘以一個2π得到一個弧度值，這種表示方法叫做自然頻率(natural frequency)：X[ω]，ω = 2π? = 2πk/N，取值范圍是0 ~ π； 4、以赫茲(Hz)為單位來表示，這個一般是應(yīng)用于一些特殊應(yīng)用，如取樣率為10 kHz表示每秒有10,000個樣本數(shù)：取值范圍是0到取樣率的一半。 二、 DFT基本函數(shù) ck[i] = cos(2πki/N) sk[i] = sin(2πki/N) 其中k表示每個正余弦波的頻率，如為2表示在0到N長度中存在兩個完整的周期，10即有10個周期，如下圖： 上圖中至于每個波的振幅(amplitude)值(Re X[k],Im X[k])是怎么算出來的,這個是DFT的核心，也是最難理解的部分，我們先來看看如何把分解出來的正余弦波合成原始信號(Inverse DFT)。 三、 合成運算方法(Real Inverse DFT) DFT合成等式（合成原始時間信號，頻率-->時間，逆向變換）： 如果有學(xué)過傅立葉級數(shù)，對這個等式就會有似曾相識的感覺，不錯！這個等式跟傅立葉級數(shù)是非常相似的： 當(dāng)然，差別是肯定是存在的，因為這兩個等式是在兩個不同條件下運用的，至于怎么證明DFT合成公式，這個我想需要非常強的高等數(shù)學(xué)理論知識了，這是研究數(shù)學(xué)的人的工作，對于普通應(yīng)用者就不需要如此的追根究底了，但是傅立葉級數(shù)是好理解的，我們起碼可以從傅立葉級數(shù)公式中看出DFT合成公式的合理性。 _ _ DFT合成等式中的Im X[k]和Re X[k]跟之前提到的Im X[k]和Re X[k]是不一樣的，下面是轉(zhuǎn)換方法（關(guān)于此公式的解釋，見下文）: 但k等于0和N/2時,實數(shù)部分的計算要用下面的等式: 上面四個式中的N是時域中點的總數(shù)，k是從0到N/2的序號。 為什么要這樣進(jìn)行轉(zhuǎn)換呢？這個可以從頻譜密度(spectral density)得到理解，如下圖就是個頻譜圖： 這是一個頻譜圖，橫坐標(biāo)表示頻率大小，縱坐標(biāo)表示振幅大小，原始信號長度為N（這里是32），經(jīng)DFT轉(zhuǎn)換后得到的17個頻率的頻譜，頻譜密度表示每單位帶寬中為多大的振幅，那么帶寬是怎么計算出來的呢？看上圖，除了頭尾兩個，其余點的所占的寬度是2/N，這個寬度便是每個點的帶寬，頭尾兩個點的帶寬是1/N,而Im X[k]和Re X[k]表示的是頻譜密度，即每一個單位帶寬的振幅大小，但 表示2/N（或1/N）帶寬的振幅大小，所以 分別應(yīng)當(dāng)是Im X[k]和Re X[k]的2/N（或1/N）。 頻譜密度就象物理中物質(zhì)密度，原始信號中的每一個點就象是一個混合物，這個混合物是由不同密度的物質(zhì)組成的，混合物中含有的每種物質(zhì)的質(zhì)量是一樣的，除了最大和最小兩個密度的物質(zhì)外，這樣我們只要把每種物質(zhì)的密度加起來就可以得到該混合物的密度了，又該混合物的質(zhì)量是單位質(zhì)量，所以得到的密度值跟該混合物的質(zhì)量值是一樣的。 至于為什么虛數(shù)部分是負(fù)數(shù)，這是為了跟復(fù)數(shù)DFT保持一致，這個我們將在后面會知道這是數(shù)學(xué)計算上的需要（Im X[k]在計算時就已經(jīng)加上了一個負(fù)號（稍后，由下文，便可知）， 再加上負(fù)號，結(jié)果便是正的，等于沒有變化）。 如果已經(jīng)得到了DFT結(jié)果，這時要進(jìn)行逆轉(zhuǎn)換，即合成原始信號，則可按如下步驟進(jìn)行轉(zhuǎn)換： 1、先根據(jù)上面四個式子計算得出 的值； 2、再根據(jù)DFT合成等式得到原始信號數(shù)據(jù)。 下面是用BASIC語言來實現(xiàn)的轉(zhuǎn)換源代碼： 100 ‘DFT逆轉(zhuǎn)換方法 110 ‘/XX[]數(shù)組存儲計算結(jié)果（時域中的原始信號） 120 ‘/REX[]數(shù)組存儲頻域中的實數(shù)分量，IMX[]為虛分量 130 ‘ 140 DIM XX[511] 150 DIM REX[256] 160 DIM IMX[256] 170 ‘ 180 PI = 3.14159265 190 N% = 512 200 ‘ 210 GOSUB XXXX ‘轉(zhuǎn)到子函數(shù)去獲取REX[]和IMX[]數(shù)據(jù) 220 ‘ 230 ‘ 240 ‘ 250 FOR K% = 0 TO 256 260 REX[K%] = REX[K%] / (N%/2) 270 IMX[K%] = -IMX[K%] / (N%/2) 280 NEXT k% 290 ‘ 300 REX[0] = REX[0] / N 310 REX[256] = REX[256] / N 320 ‘ 330 ‘ 初始化XX[]數(shù)組 340 FOR I% = 0 TO 511 350 XX[I%] = 0 360 NEXT I% 370 ‘ 380 ‘ 390 ‘ 400 ‘ 410 ‘ 420 FOR K% =0 TO 256 430 FOR I%=0 TO 511 440 ‘ 450 XX[I%] = XX[I%] + REX[K%] * COS(2 * PI * K% * I% / N%) 460 XX[I%] = XX[I%] + IMX[K%] * SIN(2 * PI * K% * I% / N%) 470 ‘ 480 NEXT I% 490 NEXT K% 500 ‘ 510 END 上面代碼中420至490換成如下形式也許更好理解，但結(jié)果都是一樣的： 420 FOR I% =0 TO 511 430 FOR K%=0 TO 256 440 ‘ 450 XX[I%] = XX[I%] + REX[K%] * COS(2 * PI * K% * I% / N%) 460 XX[I%] = XX[I%] + IMX[K%] * SIN(2 * PI * K% * I% / N%) 470 ‘ 480 NEXT I% 490 NEXT K% 四、 分解運算方法（DFT） 有三種完全不同的方法進(jìn)行DFT：一種方法是通過聯(lián)立方程進(jìn)行求解, 從代數(shù)的角度看，要從N個已知值求N個未知值，需要N個聯(lián)立方程，且N個聯(lián)立方程必須是線性獨立的，但這是這種方法計算量非常的大且極其復(fù)雜，所以很少被采用；第二種方法是利用信號的相關(guān)性（correlation）進(jìn)行計算，這個是我們后面將要介紹的方法；第三種方法是快速傅立葉變換（FFT），這是一個非常具有創(chuàng)造性和革命性的的方法，因為它大大提高了運算速度，使得傅立葉變換能夠在計算機中被廣泛應(yīng)用，但這種算法是根據(jù)復(fù)數(shù)形式的傅立葉變換來實現(xiàn)的，它把N個點的信號分解成長度為N的頻域，這個跟我們現(xiàn)在所進(jìn)行的實域DFT變換不一樣，而且這種方法也較難理解，這里我們先不去理解，等先理解了復(fù)數(shù)DFT后，再來看一下FFT。有一點很重要，那就是這三種方法所得的變換結(jié)果是一樣的，經(jīng)過實踐證明，當(dāng)頻域長度為32時，利用相關(guān)性方法進(jìn)行計算效率最好，否則FFT算法效率較高?，F(xiàn)在就讓我們來看一下相關(guān)性算法。 利用第一種方法、信號的相關(guān)性(correlation)可以從噪聲背景中檢測出已知的信號，我們也可以利用這個方法檢測信號波中是否含有某個頻率的信號波：把一個待檢測信號波乘以另一個信號波，得到一個新的信號波，再把這個新的信號波所有的點進(jìn)行相加，從相加的結(jié)果就可以判斷出這兩個信號的相似程度。如下圖： 上面a和 b兩個圖是待檢測信號波，圖a很明顯可以看出是個3個周期的正弦信號波，圖b的信號波則看不出是否含有正弦或余弦信號，圖c和d都是個3個周期的正弦信號波，圖e和f分別是a、b兩圖跟c、d兩圖相乘后的結(jié)果，圖e所有點的平均值是0.5，說明信號a含有振幅為1的正弦信號c，但圖f所有點的平均值是0，則說明信號b不含有信號d。這個就是通過信號相關(guān)性來檢測是否含有某個信號的方法。 第二種方法：相應(yīng)地，我也可以通過把輸入信號和每一種頻率的正余弦信號進(jìn)行相乘（關(guān)聯(lián)操作），從而得到原始信號與每種頻率的關(guān)聯(lián)程度（即總和大?。@個結(jié)果便是我們所要的傅立葉變換結(jié)果，下面兩個等式便是我們所要的計算方法： 第二個式子中加了個負(fù)號，是為了保持復(fù)數(shù)形式的一致，前面我們知道在計算 時又加了個負(fù)號，所以這只是個形式的問題，并沒有實際意義，你也可以把負(fù)號去掉，并在計算 時也不加負(fù)號。 這里有一點必須明白一個正交的概念：兩個函數(shù)相乘，如果結(jié)果中的每個點的總和為0，則可認(rèn)為這兩個函數(shù)為正交函數(shù)。要確保關(guān)聯(lián)性算法是正確的，則必須使得跟原始信號相乘的信號的函數(shù)形式是正交的，我們知道所有的正弦或余弦函數(shù)是正交的，這一點我們可以通過簡單的高數(shù)知識就可以證明它，所以我們可以通過關(guān)聯(lián)的方法把原始信號分離出正余弦信號。當(dāng)然，其它的正交函數(shù)也是存在的，如：方波、三角波等形式的脈沖信號，所以原始信號也可被分解成這些信號，但這只是說可以這樣做，卻是沒有用的。 下面是實域傅立葉變換的BASIC語言代碼： 到此為止，我們對傅立葉變換便有了感性的認(rèn)識了吧。但要記住，這只是在實域上的離散傅立葉變換，其中雖然也用到了復(fù)數(shù)的形式，但那只是個替代的形式，并無實際意義，現(xiàn)實中一般使用的是復(fù)數(shù)形式的離散傅立葉變換，且快速傅立葉變換是根據(jù)復(fù)數(shù)離散傅立葉變換來設(shè)計算法的，在后面我們先來復(fù)習(xí)一下有關(guān)復(fù)數(shù)的內(nèi)容，然后再在理解實域離散傅立葉變換的基礎(chǔ)上來理解復(fù)數(shù)形式的離散傅立葉變換。 前期回顧，在上一篇：十、從頭到尾徹底理解傅里葉變換算法、上里，我們講了傅立葉變換的由來、和實數(shù)形式離散傅立葉變換（Real DFT）倆個問題， 本文接上文，著重講下復(fù)數(shù)、和復(fù)數(shù)形式離散傅立葉變換等倆個問題。 第三章、復(fù)數(shù) 復(fù)數(shù)擴展了我們一般所能理解的數(shù)的概念，復(fù)數(shù)包含了實數(shù)和虛數(shù)兩部分，利用復(fù)數(shù)的形式可以把由兩個變量表示的表達(dá)式變成由一個變量(復(fù)變量)來表達(dá)，使得處理起來更加自然和方便。 我們知道傅立葉變換的結(jié)果是由兩部分組成的，使用復(fù)數(shù)形式可以縮短變換表達(dá)式，使得我們可以單獨處理一個變量（這個在后面的描述中我們就可以更加確切地知道），而且快速傅立葉變換正是基于復(fù)數(shù)形式的，所以幾乎所有描述的傅立葉變換形式都是復(fù)數(shù)的形式。 但是復(fù)數(shù)的概念超過了我們?nèi)粘Ｉ钪兴芾斫獾母拍睿斫鈴?fù)數(shù)是較難的，所以我們在理解復(fù)數(shù)傅立葉變換之前，先來專門復(fù)習(xí)一下有關(guān)復(fù)數(shù)的知識，這對后面的理解非常重要。 一、 復(fù)數(shù)的提出 在此，先讓我們看一個物理實驗：把一個球從某點向上拋出，然后根據(jù)初速度和時間來計算球所在高度，這個方法可以根據(jù)下面的式子計算得出： 其中h表示高度，g表示重力加速度(9.8m/s2)，v表示初速度，t表示時間?，F(xiàn)在反過來，假如知道了高度，要求計算到這個高度所需要的時間，這時我們又可以通過下式來計算： （多謝JERRY_PRI提出： 1、根據(jù)公式h=-(gt2/2)+Vt（gt后面的2表示t的平方），我們可以討論最終情況，也就是說小球運動到最高點時，v=gt，所以，可以得到t=sqt(2h/g) 且在您給的公式中，根號下為1-(2h)/g，化成分?jǐn)?shù)形式為(g-2h)/g，g和h不能直接做加減運算。 2、g是重力加速度，單位是m/s2，h的單位是m，他們兩個相減的話在物理上沒有意義，而且使用您給的那個公式反向回去的話推出的是h=-(gt2/2)+gt啊（gt后面的2表示t的平方）。 3、直接推到可以得出t=v/g±sqt((v2-2hg)/g2)（v和g后面的2都表示平方），那么也就是說當(dāng)v2<2hg時會產(chǎn)生復(fù)數(shù)，但是如果從實際的v2是不可能小于2hg的，所以我感覺復(fù)數(shù)不能從實際出發(fā)去推到，只能從抽象的角度說明一下。 ） 經(jīng)過計算我們可以知道，當(dāng)高度是3米時，有兩個時間點到達(dá)該高度：球向上運動時的時間是0.38秒，球向下運動時的時間是1.62秒。但是如果高度等于10時，結(jié)果又是什么呢？根據(jù)上面的式子可以發(fā)現(xiàn)存在對負(fù)數(shù)進(jìn)行開平方運算，我們知道這肯定是不現(xiàn)實的。 第一次使用這個不一般的式子的人是意大利數(shù)學(xué)家Girolamo Cardano（1501-1576），兩個世紀(jì)后，德國偉大數(shù)學(xué)家Carl Friedrich Gause（1777-1855）提出了復(fù)數(shù)的概念，為后來的應(yīng)用鋪平了道路，他對復(fù)數(shù)進(jìn)行這樣表示：復(fù)數(shù)由實數(shù)（real）和虛數(shù)(imaginary)兩部分組成，虛數(shù)中的根號負(fù)1用i來表示（在這里我們用j來表示，因為i在電力學(xué)中表示電流的意思）。 我們可以把橫坐標(biāo)表示成實數(shù)，縱坐標(biāo)表示成虛數(shù)，則坐標(biāo)中的每個點的向量就可以用復(fù)數(shù)來表示，如下圖： 上圖中的ABC三個向量可以表示成如下的式子： A = 2 + 6j B = -4 – 1.5j C = 3 – 7j 這樣子來表達(dá)方便之處在于運用一個符號就能把兩個原來難以聯(lián)系起來的數(shù)組合起來了，不方便的是我們要分辨哪個是實數(shù)和哪個是虛數(shù)，我們一般是用Re( )和Im( )來表示實數(shù)和虛數(shù)兩部分，如： Re A = 2 Im A = 6 Re B = -4 Im B = -1.5 Re C = 3 Im C = -7 復(fù)數(shù)之間也可以進(jìn)行加減乘除運算： 這里有個特殊的地方是j2等于-1，上面第四個式子的計算方法是把分子和分母同時乘以c – dj，這樣就可消去分母中的j了。 復(fù)數(shù)也符合代數(shù)運算中的交換律、結(jié)合律、分配律： A B = B A (A + B) + C = A + (B + C) A(B + C) = AB + AC 二、 復(fù)數(shù)的極坐標(biāo)表示形式 前面提到的是運用直角坐標(biāo)來表示復(fù)數(shù)，其實更為普遍應(yīng)用的是極坐標(biāo)的表示方法，如下圖： 上圖中的M即是數(shù)量積(magnitude)，表示從原點到坐標(biāo)點的距離，θ是相位角(phase angle)，表示從X軸正方向到某個向量的夾角，下面四個式子是計算方法： 我們還可以通過下面的式子進(jìn)行極坐標(biāo)到直角坐標(biāo)的轉(zhuǎn)換： a + jb = M (cosθ + j sinθ) 上面這個等式中左邊是直角坐標(biāo)表達(dá)式，右邊是極坐標(biāo)表達(dá)式。 還有一個更為重要的等式——歐拉等式（歐拉，瑞士的著名數(shù)學(xué)家，Leonhard Euler，1707-1783）： ejx = cos x + j sin x 這個等式可以從下面的級數(shù)變換中得到證明： 上面中右邊的兩個式子分別是cos(x)和sin(x)的泰勒(Taylor)級數(shù)。 這樣子我們又可以把復(fù)數(shù)的表達(dá)式表示成指數(shù)的形式了： a + jb = M ejθ （這便是復(fù)數(shù)的兩個表達(dá)式） 指數(shù)形式是數(shù)字信號處理中數(shù)學(xué)方法的支柱，也許是因為用指數(shù)形式進(jìn)行復(fù)數(shù)的乘除運算極為簡單的緣故吧： 三、復(fù)數(shù)是數(shù)學(xué)分析中的一個工具 為什么要使用復(fù)數(shù)呢？其實它只是個工具而已，就如釘子和錘子的關(guān)系，復(fù)數(shù)就象那錘子，作為一種使用的工具。我們把要解決的問題表達(dá)成復(fù)數(shù)的形式（因為有些問題用復(fù)數(shù)的形式進(jìn)行運算更加方便），然后對復(fù)數(shù)進(jìn)行運算，最后再轉(zhuǎn)換回來得到我們所需要的結(jié)果。 有兩種方法使用復(fù)數(shù)，一種是用復(fù)數(shù)進(jìn)行簡單的替換，如前面所說的向量表達(dá)式方法和前一節(jié)中我們所討論的實域DFT，另一種是更高級的方法：數(shù)學(xué)等價(mathematical equivalence)，復(fù)數(shù)形式的傅立葉變換用的便是數(shù)學(xué)等價的方法，但在這里我們先不討論這種方法，這里我們先來看一下用復(fù)數(shù)進(jìn)行替換中的問題。 用復(fù)數(shù)進(jìn)行替換的基本思想是：把所要分析的物理問題轉(zhuǎn)換成復(fù)數(shù)的形式，其中只是簡單地添加一個復(fù)數(shù)的符號j，當(dāng)返回到原來的物理問題時，則只是把符號j去掉就可以了。 有一點要明白的是并不是所有問題都可以用復(fù)數(shù)來表示，必須看用復(fù)數(shù)進(jìn)行分析是否適用，有個例子可以看出用復(fù)數(shù)來替換原來問題的表達(dá)方式明顯是謬誤的：假設(shè)一箱的蘋果是5美元，一箱的桔子是10美元，于是我們把它表示成 5 + 10j，有一個星期你買了6箱蘋果和2箱桔子，我們又把它表示成6 + 2j，最后計算總共花的錢是(5 + 10j)(6 + 2j) = 10 + 70j，結(jié)果是買蘋果花了10美元的，買桔子花了70美元，這樣的結(jié)果明顯是錯了，所以復(fù)數(shù)的形式不適合運用于對這種問題的解決。 四、用復(fù)數(shù)來表示正余弦函數(shù)表達(dá)式 對于象M cos (ωt + φ)和A cos(ωt ) + B sin(ωt )表達(dá)式，用復(fù)數(shù)來表示，可以變得非常簡潔，對于直角坐標(biāo)形式可以按如下形式進(jìn)行轉(zhuǎn)換： 上式中余弦幅值A(chǔ)經(jīng)變換生成a，正弦幅值B的相反數(shù)經(jīng)變換生成b：A <=> a，B<=> -b，但要注意的是，這不是個等式，只是個替換形式而已。 對于極坐標(biāo)形式可以按如下形式進(jìn)行轉(zhuǎn)換： 上式中，M <=> M，θ<=>φ。 這里虛數(shù)部分采用負(fù)數(shù)的形式主要是為了跟復(fù)數(shù)傅立葉變換表達(dá)式保持一致，對于這種替換的方法來表示正余弦，符號的變換沒有什么好處，但替換時總會被改變掉符號以跟更高級的等價變換保持形式上的一致。 在離散信號處理中，運用復(fù)數(shù)形式來表示正余弦波是個常用的技術(shù)，這是因為利用復(fù)數(shù)進(jìn)行各種運算得到的結(jié)果跟原來的正余弦運算結(jié)果是一致的，但是，我們要小心使用復(fù)數(shù)操作，如加、減、乘、除，有些操作是不能用的，如兩個正弦信號相加，采用復(fù)數(shù)形式進(jìn)行相加，得到的結(jié)果跟替換前的直接相加的結(jié)果是一樣的，但是如果兩個正弦信號相乘，則采用復(fù)數(shù)形式來相乘結(jié)果是不一樣的。幸運的是，我們已嚴(yán)格定義了正余弦復(fù)數(shù)形式的運算操作條件： 1、參加運算的所有正余弦的頻率必須是一樣的； 2、運算操作必須是線性的，如兩個正弦信號可以進(jìn)行相加減，但不能進(jìn)行乘除，象信號的放大、衰減、高低通濾波等系統(tǒng)都是線性的，象平方、縮短、取限等則不是線性的。要記住的是卷積和傅立葉分析也只有線性操作才可以進(jìn)行。 下圖是一個相量變換(我們把正弦或余弦波變成復(fù)數(shù)的形式稱為相量變換，Phasor transform)的例子，一個連續(xù)信號波經(jīng)過一個線性處理系統(tǒng)生成另一個信號波，從計算過程我們可以看出采用復(fù)數(shù)的形式使得計算變化十分的簡潔： 在第二章中我們描述的實數(shù)形式傅立葉變換也是一種替換形式的復(fù)數(shù)變換，但要注意的是那還不是復(fù)數(shù)傅立葉變換，只是一種代替方式而已。下一章、即，第四章，我們就會知道復(fù)數(shù)傅立葉變換是一種更高級的變換，而不是這種簡單的替換形式。 第四章、復(fù)數(shù)形式離散傅立葉變換 復(fù)數(shù)形式的離散傅立葉變換非常巧妙地運用了復(fù)數(shù)的方法，使得傅立葉變換變換更加自然和簡潔，它并不是只是簡單地運用替換的方法來運用復(fù)數(shù)，而是完全從復(fù)數(shù)的角度來分析問題，這一點跟實數(shù)DFT是完全不一樣的。 一、 把正余弦函數(shù)表示成復(fù)數(shù)的形式 通過歐拉等式可以把正余弦函數(shù)表示成復(fù)數(shù)的形式： cos( x ) = 1/2 e j(-x) + 1/2 ejx sin( x ) = j (1/2 e j(-x) - 1/2 ejx) 從這個等式可以看出，如果把正余弦函數(shù)表示成復(fù)數(shù)后，它們變成了由正負(fù)頻率組成的正余弦波，相反地，一個由正負(fù)頻率組成的正余弦波，可以通過復(fù)數(shù)的形式來表示。 我們知道，在實數(shù)傅立葉變換中，它的頻譜是0 ~ π(0 ~ N/2),但無法表示-π~ 0的頻譜，可以預(yù)見，如果把正余弦表示成復(fù)數(shù)形式，則能夠把負(fù)頻率包含進(jìn)來。 二、 把變換前后的變量都看成復(fù)數(shù)的形式 復(fù)數(shù)形式傅立葉變換把原始信號x[n]當(dāng)成是一個用復(fù)數(shù)來表示的信號，其中實數(shù)部分表示原始信號值，虛數(shù)部分為0，變換結(jié)果X[k]也是個復(fù)數(shù)的形式，但這里的虛數(shù)部分是有值的。 在這里要用復(fù)數(shù)的觀點來看原始信號，是理解復(fù)數(shù)形式傅立葉變換的關(guān)鍵（如果有學(xué)過復(fù)變函數(shù)則可能更好理解，即把x[n]看成是一個復(fù)數(shù)變量，然后象對待實數(shù)那樣對這個復(fù)數(shù)變量進(jìn)行相同的變換）。 三、 對復(fù)數(shù)進(jìn)行相關(guān)性算法（正向傅立葉變換） 從實數(shù)傅立葉變換中可以知道，我們可以通過原始信號乘以一個正交函數(shù)形式的信號，然后進(jìn)行求總和，最后就能得到這個原始信號所包含的正交函數(shù)信號的分量。 現(xiàn)在我們的原始信號變成了復(fù)數(shù)，我們要得到的當(dāng)然是復(fù)數(shù)的信號分量，我們是不是可以把它乘以一個復(fù)數(shù)形式的正交函數(shù)呢？答案是肯定的，正余弦函數(shù)都是正交函數(shù)，變成如下形式的復(fù)數(shù)后，仍舊還是正交函數(shù)（這個從正交函數(shù)的定義可以很容易得到證明）： cos x + j sin x, cos x – j sin x，…… 這里我們采用上面的第二個式子進(jìn)行相關(guān)性求和，為什么用第二個式子呢?，我們在后面會知道，正弦函數(shù)在虛數(shù)中變換后得到的是負(fù)的正弦函數(shù)，這里我們再加上一個負(fù)號，使得最后的得到的是正的正弦波，根據(jù)這個于是我們很容易就可以得到了復(fù)數(shù)形式的DFT正向變換等式： 這個式子很容易可以得到歐拉變換式子： 其實我們是為了表達(dá)上的方便才用到歐拉變換式，在解決問題時我們還是較多地用到正余弦表達(dá)式。 對于上面的等式，我們要清楚如下幾個方面（也是區(qū)別于實數(shù)DFT的地方）： 1、X[k]、x[n]都是復(fù)數(shù)，但x[n]的虛數(shù)部分都是由0組成的，實數(shù)部分表示原始信號； 2、k的取值范圍是0 ~ N-1 (也可以表達(dá)成0 ~ 2π)，其中0 ~ N/2（或0 ~ π）是正頻部分， N/2 ~ N-1（π~ 2π）是負(fù)頻部分，由于正余弦函數(shù)的對稱性，所以我們把 –π~ 0表示成π~ 2π，這是出于計算上方便的考慮。 3、其中的j是一個不可分離的組成部分，就象一個等式中的變量一樣，不能隨便去掉，去掉之后意義就完全不一樣了，但我們知道在實數(shù)DFT中，j只是個符號而已，把j去掉，整個等式的意義不變； 4、下圖是個連續(xù)信號的頻譜，但離散頻譜也是與此類似的，所以不影響我們對問題的分析： 上面的頻譜圖把負(fù)頻率放到了左邊，是為了迎合我們的思維習(xí)慣，但在實際實 現(xiàn)中我們一般是把它移到正的頻譜后面的。 從上圖可以看出，時域中的正余弦波（用來組成原始信號的正余弦波）在復(fù)數(shù)DFT的頻譜中被分成了正、負(fù)頻率的兩個組成部分，基于此等式中前面的比例系數(shù)是1/N（或1/2π），而不是2/N，這是因為現(xiàn)在把頻譜延伸到了2π,但把正負(fù)兩個頻率相加即又得到了2/N,又還原到了實數(shù)DFT的形式，這個在后面的描述中可以更清楚地看到。 由于復(fù)數(shù)DFT生成的是一個完整的頻譜，原始信號中的每一個點都是由正、負(fù)兩個頻率組合而成的，所以頻譜中每一個點的帶寬是一樣的，都是1/N，相對實數(shù)DFT，兩端帶寬比其它點的帶寬少了一半；復(fù)數(shù)DFT的頻譜特征具有周期性：-N/2 ~ 0與N/2 ~ N-1是一樣的，實域頻譜呈偶對稱性（表示余弦波頻譜），虛域頻譜呈奇對稱性（表示正弦波頻譜）。 四、 逆向傅立葉變換 假設(shè)我們已經(jīng)得到了復(fù)數(shù)形式的頻譜X[k]，現(xiàn)在要把它還原到復(fù)數(shù)形式的原始信號x[n]，當(dāng)然應(yīng)該是把X[k]乘以一個復(fù)數(shù)，然后再進(jìn)行求和，最后得到原始信號x[n]，這個跟X[k]相乘的復(fù)數(shù)首先讓我們想到的應(yīng)該是上面進(jìn)行相關(guān)性計算的復(fù)數(shù)： cos(2πkn/N) – j si(2πkn/N)， 但其中的負(fù)號其實是為了使得進(jìn)行逆向傅立葉變換時把正弦函數(shù)變?yōu)檎姆枺驗樘摂?shù)j的運算特殊性，使得原來應(yīng)該是正的正弦函數(shù)變?yōu)榱素?fù)的正弦函數(shù)（我們從后面的推導(dǎo)會看到這一點），所以這里的負(fù)號只是為了糾正符號的作用，在進(jìn)行逆向DFT時，我們可以把負(fù)號去掉，于是我們便得到了這樣的逆向DFT變換等式： x[n] = X[k] (cos(2πkn/N) + j sin(2πkn/N)) 我們現(xiàn)在來分析這個式子，會發(fā)現(xiàn)這個式其實跟實數(shù)傅立葉變換是可以得到一樣結(jié)果的。我們先把X[k]變換一下： X[k] = Re X[k] + j Im X[k] 這樣我們就可以對x[n]再次進(jìn)行變換，如： x[n] = (Re X[k] + j Im X[k]) (cos(2πkn/N) + j sin(2πkn/N)) = ( Re X[k] cos(2πkn/N) + j Im X[k] cos(2πkn/N) +j Re X[k] sin(2πkn/N) - Im X[k] sin(2πkn/N) ) = ( Re X[k] （cos(2πkn/N) + j sin(2πkn/N)） + ---------------------(1) Im X[k] （ - sin(2πkn/N) + j cos(2πkn/N))） ---------------(2) 這時我們就把原來的等式分成了兩個部分，第一個部分是跟實域中的頻譜相乘，第二個部分是跟虛域中的頻譜相乘，根據(jù)頻譜圖我們可以知道，Re X[k]是個偶對稱的變量，Im X[k]是個奇對稱的變量，即 Re X[k] = Re X[- k] Im X[k] = - Im X[-k] 但k的范圍是0 ~ N-1，0~N/2表示正頻率，N/2~N-1表示負(fù)頻率，為了表達(dá)方便我們把N/2~N-1用-k來表示，這樣在從0到N-1的求和過程中對于(1)和(2)式分別有N/2對的k和-k的和，對于（1）式有： Re X[k] (cos(2πkn/N) + j sin(2πkn/N)) + Re X[- k] (cos( - 2πkn/N) + j sin( -2πkn/N)) 根據(jù)偶對稱性和三角函數(shù)的性質(zhì)，把上式化簡得到： Re X[k] (cos(2πkn/N) + j sin(2πkn/N)) + Re X[ k] (cos( 2πkn/N) - j sin( 2πkn/N)) 這個式子最后的結(jié)果是： 2 Re X[ k] cos(2πkn/N)。 再考慮到求Re X[ k]等式中有個比例系數(shù)1/N，把1/N乘以2，這樣的結(jié)果不就是跟實數(shù)DFT中的式子一樣了嗎？ 對于(2)式，用同樣的方法，我們也可以得到這樣的結(jié)果： -2 Im X[k] sin(2πkn/N) 注意上式前面多了個負(fù)符號，這是由于虛數(shù)變換的特殊性造成的，當(dāng)然我們肯定不能把負(fù)符號的正弦函數(shù)跟余弦來相加，還好，我們前面是用cos(2πkn/N) – j sin(2πkn/N)進(jìn)行相關(guān)性計算，得到的Im X[k]中有個負(fù)的符號，這樣最后的結(jié)果中正弦函數(shù)就沒有負(fù)的符號了，這就是為什么在進(jìn)行相關(guān)性計算時虛數(shù)部分要用到負(fù)符號的原因（我覺得這也許是復(fù)數(shù)形式DFT美中不足的地方，讓人有一種拼湊的感覺）。 從上面的分析中可以看出，實數(shù)傅立葉變換跟復(fù)數(shù)傅立葉變換，在進(jìn)行逆變換時得到的結(jié)果是一樣的，只不過是殊途同歸吧。本文完。（July、dznlong） 本人July對本博客所有任何文章、內(nèi)容和資料享有版權(quán)。 轉(zhuǎn)載務(wù)必注明作者本人及出處，并通知本人。 |
|
|
