1. 程式人生 > >全面解析傅立葉變換(非常詳細)

全面解析傅立葉變換(非常詳細)

前言

第一部分、  DFT
第一章、傅立葉變換的由來
第二章、實數形式離散傅立葉變換(Real DFT)

第三章、複數

第四章、複數形式離散傅立葉變換

前言
“關於傅立葉變換,無論是書本還是在網上可以很容易找到關於傅立葉變換的描述,但是大都是些故弄玄虛的文章,太過抽象,盡是一些讓人看了就望而生畏的公式的羅列,讓人很難能夠從感性上得到理解”---dznlong,

那麼,到底什麼是傅立葉變換演算法列?傅立葉變換所涉及到的公式具體有多複雜列?
傅立葉變換(Fourier transform)是一種線性的積分變換。因其基本思想首先由法國學者傅立葉系統地提出,所以以其名字來命名以示紀念。

   哦,傅立葉變換原來就是一種變換而已,只是這種變換是從時間轉換為頻率的變化。這下,你就知道了,傅立葉就是一種變換,一種什麼變換列?就是一種從時間到頻率的變化或其相互轉化。

  ok,咱們再來總體瞭解下傅立葉變換,讓各位對其有個總體大概的印象,也順便看看傅立葉變換所涉及到的公式,究竟有多複雜:
以下就是傅立葉變換的4種變體(摘自,維基百科)
連續傅立葉變換
   一般情況下,若“傅立葉變換”一詞不加任何限定語,則指的是“連續傅立葉變換”。連續傅立葉變換將平方可積的函式f(t)表示成復指數函式的積分或級數形式。

這是將頻率域的函式F(ω)表示為時間域的函式f(t)的積分形式。

連續傅立葉變換的逆變換 (inverse Fourier transform)為:

即將時間域的函式f(t)表示為頻率域的函式F(ω)的積分。

一般可稱函式f(t)為原函式,而稱函式F(ω)

為傅立葉變換的像函式,原函式和像函式構成一個傅立葉變換對(transform pair)。

除此之外,還有其它型式的變換對,以下兩種型式亦常被使用。在通訊或是訊號處理方面,常以來代換,而形成新的變換對:

 或者是因係數重分配而得到新的變換對:

 一種對連續傅立葉變換的推廣稱為分數傅立葉變換(Fractional Fourier Transform)。分數傅立葉變換(fractional Fourier transform,FRFT)指的就是傅立葉變換(Fourier transform,FT)的廣義化。
分數傅立葉變換的物理意義即做傅立葉變換 a 次,其中 a 不一定要為整數;而做了分數傅立葉變換之後,訊號或輸入函式便會出現在介於時域(time domain)與頻域(frequency domain)之間的分數域(fractional domain)。

當f(t)為偶函式(或奇函式)時,其正弦(或餘弦)分量將消亡,而可以稱這時的變換為餘弦變換(cosine transform)或正弦變換(sine transform).

另一個值得注意的性質是,當f(t)為純實函式時,F(−ω) = F*(ω)成立.

傅立葉級數
   連續形式的傅立葉變換其實是傅立葉級數 (Fourier series)的推廣,因為積分其實是一種極限形式的求和運算元而已。對於周期函式,其傅立葉級數是存在的:

其中Fn為復幅度。對於實值函式,函式的傅立葉級數可以寫成:


其中an和bn是實頻率分量的幅度。

離散時域傅立葉變換
   離散傅立葉變換是離散時間傅立葉變換(DTFT)的特例(有時作為後者的近似)。DTFT在時域上離散,在頻域上則是週期的。DTFT可以被看作是傅立葉級數的逆變換。

離散傅立葉變換
   離散傅立葉變換(DFT),是連續傅立葉變換在時域和頻域上都離散的形式,將時域訊號的取樣變換為在離散時間傅立葉變換(DTFT)頻域的取樣。在形式上,變換兩端(時域和頻域上)的序列是有限長的,而實際上這兩組序列都應當被認為是離散週期訊號的主值序列。即使對有限長的離散訊號作DFT,也應當將其看作經過週期延拓成為週期訊號再作變換。在實際應用中通常採用快速傅立葉變換以高效計算DFT。

   為了在科學計算和數字訊號處理等領域使用計算機進行傅立葉變換,必須將函式xn定義在離散點而非連續域內,且須滿足有限性或週期性條件。這種情況下,使用離散傅立葉變換(DFT),將函式xn表示為下面的求和形式:

其中Xk是傅立葉幅度。直接使用這個公式計算的計算複雜度為O(n*n),而快速傅立葉變換(FFT)可以將複雜度改進為O(n*lgn)。(後面會具體闡述FFT是如何將複雜度降為O(n*lgn)的。)計算複雜度的降低以及數位電路計算能力的發展使得DFT成為在訊號處理領域十分實用且重要的方法。

   下面,比較下上述傅立葉變換的4種變體,

   如上,容易發現:函式在時(頻)域的離散對應於其像函式在頻(時)域的週期性。反之連續則意味著在對應域的訊號的非週期性。也就是說,時間上的離散性對應著頻率上的週期性。同時,注意,離散時間傅立葉變換,時間離散,頻率不離散,它在頻域依然是連續的。
   如果,讀到此,你不甚明白,大沒關係,不必糾結於以上4種變體,繼續往下看,你自會豁然開朗。(有什麼問題,也懇請提出,或者批評指正)

   ok, 本文,接下來,由傅立葉變換入手,後重點闡述離散傅立葉變換、快速傅立葉演算法,到最後徹底實現FFT演算法,全篇力求通俗易懂、閱讀順暢,教你從頭到尾徹底理解傅立葉變換演算法。由於傅立葉變換,也稱傅立葉變換,下文所稱為傅立葉變換,同一個變換,不同叫法,讀者不必感到奇怪。

第一部分、DFT
第一章、傅立葉變換的由來
    要理解傅立葉變換,先得知道傅立葉變換是怎麼變換的,當然,也需要一定的高等數學基礎,最基本的是級數變換,其中傅立葉級數變換是傅立葉變換的基礎公式。
 
一、傅立葉變換的提出

傅立葉是一位法國數學家和物理學家,原名是Jean Baptiste Joseph Fourier(1768-1830), Fourier於1807年在法國科學學會上發表了一篇論文,論文裡描述運用正弦曲線來描述溫度分佈,論文裡有個在當時具有爭議性的決斷:任何連續週期訊號都可以由一組適當的正弦曲線組合而成。

    當時審查這個論文拉格朗日堅決反對此論文的發表,而後在近50年的時間裡,拉格朗日堅持認為傅立葉的方法無法表示帶有稜角的訊號,如在方波中出現非連續變化斜率。直到拉格朗日死後15年這個論文才被髮表出來。
    誰是對的呢?拉格朗日是對的:正弦曲線無法組合成一個帶有稜角的訊號。但是,我們可以用正弦曲線來非常逼近地表示它,逼近到兩種表示方法不存在能量差別,基於此,傅立葉是對的。

    為什麼我們要用正弦曲線來代替原來的曲線呢?如我們也還可以用方波或三角波來代替呀,分解訊號的方法是無窮多的,但分解訊號的目的是為了更加簡單地處理原來的訊號。
    用正餘弦來表示原訊號會更加簡單,因為正餘弦擁有原訊號所不具有的性質:正弦曲線保真度。一個正餘弦曲線訊號輸入後,輸出的仍是正餘弦曲線,只有幅度和相位可能發生變化,但是頻率和波的形狀仍是一樣的。且只有正餘弦曲線才擁有這樣的性質,正因如此我們才不用方波或三角波來表示。

二、傅立葉變換分類
    根據原訊號的不同型別,我們可以把傅立葉變換分為四種類別:
1、非週期性連續訊號        傅立葉變換(Fourier Transform) 
2、週期性連續訊號           傅立葉級數(Fourier Series) 
3、非週期性離散訊號        離散時域傅立葉變換(Discrete Time Fourier Transform) 
4、週期性離散訊號           離散傅立葉變換(Discrete Fourier Transform) 
       下圖是四種原訊號圖例(從上到下,依次是FT,FS,DTFT,DFT): 

 

    這四種傅立葉變換都是針對正無窮大和負無窮大的訊號,即訊號的的長度是無窮大的,我們知道這對於計算機處理來說是不可能的,那麼有沒有針對長度有限的傅立葉變換呢?沒有。因為正餘弦波被定義成從負無窮小到正無窮大,我們無法把一個長度無限的訊號組合成長度有限的訊號。
    面對這種困難,方法是:把長度有限的訊號表示成長度無限的訊號。如,可以把訊號無限地從左右進行延伸,延伸的部分用零來表示,這樣,這個訊號就可以被看成是非週期性離散訊號我們可以用到離散時域傅立葉變換(DTFT)的方法。也可以把訊號用複製的方法進行延伸,這樣訊號就變成了週期性離散訊號,這時我們就可以用離散傅立葉變換方法(DFT)進行變換。本章我們要講的是離散訊號,對於連續訊號我們不作討論,因為計算機只能處理離散的數值訊號,我們的最終目的是運用計算機來處理訊號的。
 
    但是對於非週期性的訊號,我們需要用無窮多不同頻率的正弦曲線來表示,這對於計算機來說是不可能實現的。所以對於離散訊號的變換隻有離散傅立葉變換(DFT)才能被適用,對於計算機來說只有離散的和有限長度的資料才能被處理,對於其它的變換型別只有在數學演算中才能用到,在計算機面前我們只能用DFT方法,後面我們要理解的也正是DFT方法。
    這裡要理解的是我們使用週期性的訊號目的是為了能夠用數學方法來解決問題,至於考慮週期性訊號是從哪裡得到或怎樣得到是無意義的。
 
    每種傅立葉變換都分成實數和複數兩種方法,對於實數方法是最好理解的,但是複數方法就相對複雜許多了,需要懂得有關複數的理論知識,不過,如果理解了實數離散傅立葉變換(real DFT),再去理解複數傅立葉變換就更容易了,所以我們先把複數的傅立葉變換放到一邊去,先來理解實數傅立葉變換,在後面我們會先講講關於複數的基本理論,然後在理解了實數傅立葉變換的基礎上再來理解複數傅立葉變換。
 
    還有,這裡我們所要說的變換(transform)雖然是數學意義上的變換,但跟函式變換是不同的,函式變換是符合一一映射準則的,對於離散數字訊號處理(DSP),有許多的變換:傅立葉變換、拉普拉斯變換、Z變換、希爾伯特變換、離散餘弦變換等,這些都擴充套件了函式變換的定義,允許輸入和輸出有多種的值,簡單地說變換就是把一堆的資料變成另一堆的資料的方法。
 
三、一個關於實數離散傅立葉變換(Real DFT)的例子

       先來看一個變換例項,下圖是一個原始訊號影象:


       這個訊號的長度是16,於是可以把這個訊號分解9個餘弦波和9個正弦波(一個長度為N的訊號可以分解成N/2+1個正餘弦訊號,這是為什麼呢?結合下面的18個正餘弦圖,我想從計算機處理精度上就不難理解,一個長度為N的訊號,最多隻能有N/2+1個不同頻率,再多的頻率就超過了計算機所能所處理的精度範圍),如下圖:

        9個餘弦訊號:

        9個正弦訊號:

       把以上所有訊號相加即可得到原始訊號,至於是怎麼分別變換出9種不同頻率訊號的,我們先不急,先看看對於以上的變換結果,在程式中又是該怎麼表示的,我們可以看看下面這個示例圖:


 
    上圖中左邊表示時域中的訊號,右邊是頻域訊號表示方法,
從左向右,-->,表示正向轉換(Forward DFT),從右向左,<--,表示逆向轉換(Inverse DFT),
用小寫x[]表示訊號在每個時間點上的幅度值陣列, 用大寫X[]表示每種頻率的副度值陣列(即時間x-->頻率X), 
因為有N/2+1種頻率,所以該陣列長度為N/2+1,
    X[]陣列又分兩種,一種是表示餘弦波的不同頻率幅度值:Re X[],
另一種是表示正弦波的不同頻率幅度值:Im X[],
    Re是實數(Real)的意思,Im是虛數(Imagine)的意思,採用複數的表示方法把正餘弦波組合起來進行表示,但這裡我們不考慮複數的其它作用,只記住是一種組合方法而已,目的是為了便於表達(在後面我們會知道,複數形式的傅立葉變換長度是N,而不是N/2+1)。如此,再回過頭去,看上面的正餘弦各9種頻率的變化,相信,問題不大了。

第二章、實數形式離散傅立葉變換(Real DFT)
       上一章,我們看到了一個實數形式離散傅立葉變換的例子,通過這個例子能夠讓我們先對傅立葉變換有一個較為形象的感性認識,現在就讓我們來看看實數形式離散傅立葉變換的正向和逆向是怎麼進行變換的。在此,我們先來看一下頻率的多種表示方法。
 
一、   頻域中關於頻率的四種表示方法
 
1、序號表示方法,根據時域中訊號的樣本數取0 ~ N/2,用這種方法在程式中使用起來可以更直接地取得每種頻率的幅度值,因為頻率值跟陣列的序號是一一對應的: X[k],取值範圍是0 ~ N/2;
2、分數表示方法,根據時域中訊號的樣本數的比例值取0 ~ 0.5: X[ƒ],ƒ = k/N,取值範圍是0 ~ 1/2;
3、用弧度值來表示,把ƒ乘以一個2π得到一個弧度值,這種表示方法叫做自然頻率(natural frequency):X[ω],ω = 2πƒ = 2πk/N,取值範圍是0 ~ π;
4、以赫茲(Hz)為單位來表示,這個一般是應用於一些特殊應用,如取樣率為10 kHz表示每秒有10,000個樣本數:取值範圍是0到取樣率的一半。
 
二、   DFT基本函式
 
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合成等式(合成原始時間訊號,頻率-->時間,逆向變換):

如果有學過傅立葉級數,對這個等式就會有似曾相識的感覺,不錯!這個等式跟傅立葉級數是非常相似的:

           當然,差別是肯定是存在的,因為這兩個等式是在兩個不同條件下運用的,至於怎麼證明DFT合成公式,這個我想需要非常強的高等數學理論知識了,這是研究數學的人的工作,對於普通應用者就不需要如此的追根究底了,但是傅立葉級數是好理解的,我們起碼可以從傅立葉級數公式中看出DFT合成公式的合理性。
                                  _            _
       DFT合成等式中的Im X[k]和Re X[k]跟之前提到的Im X[k]和Re X[k]是不一樣的,下面是轉換方法(關於此公式的解釋,見下文):

       
       但k等於0和N/2時,實數部分的計算要用下面的等式:

              
       上面四個式中的N是時域中點的總數,k是從0到N/2的序號。
       為什麼要這樣進行轉換呢?這個可以從頻譜密度(spectral density)得到理解,如下圖就是個頻譜圖:

       
       這是一個頻譜圖,橫座標表示頻率大小,縱座標表示振幅大小,原始訊號長度為N(這裡是32),經DFT轉換後得到的17個頻率的頻譜,頻譜密度表示每單位頻寬中為多大的振幅,那麼頻寬是怎麼計算出來的呢?看上圖,除了頭尾兩個,其餘點的所佔的寬度是2/N,這個寬度便是每個點的頻寬,頭尾兩個點的頻寬是1/N,而Im X[k]和Re X[k]表示的是頻譜密度,即每一個單位頻寬的振幅大小,但表示2/N(或1/N)頻寬的振幅大小,所以分別應當是Im X[k]和Re X[k]的2/N(或1/N)。
 
頻譜密度就象物理中物質密度,原始訊號中的每一個點就象是一個混合物,這個混合物是由不同密度的物質組成的,混合物中含有的每種物質的質量是一樣的,除了最大和最小兩個密度的物質外,這樣我們只要把每種物質的密度加起來就可以得到該混合物的密度了,又該混合物的質量是單位質量,所以得到的密度值跟該混合物的質量值是一樣的。
 
       至於為什麼虛數部分是負數,這是為了跟複數DFT保持一致,這個我們將在後面會知道這是數學計算上的需要(Im X[k]在計算時就已經加上了一個負號(稍後,由下文,便可知),再加上負號,結果便是正的,等於沒有變化)。
 
       如果已經得到了DFT結果,這時要進行逆轉換,即合成原始訊號,則可按如下步驟進行轉換:
1、先根據上面四個式子計算得出的值;
2、再根據DFT合成等式得到原始訊號資料。
下面是用BASIC語言來實現的轉換原始碼:
100 ‘DFT逆轉換方法
110 ‘/XX[]陣列儲存計算結果(時域中的原始訊號)
120 ‘/REX[]陣列儲存頻域中的實數分量,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 ‘轉到子函式去獲取REX[]和IMX[]資料
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[]陣列
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換成如下形式也許更好理解,但結果都是一樣的:
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)
 
      有三種完全不同的方法進行DFT:一種方法是通過聯立方程進行求解, 從代數的角度看,要從N個已知值求N個未知值,需要N個聯立方程,且N個聯立方程必須是線性獨立的,但這是這種方法計算量非常的大且極其複雜,所以很少被採用;第二種方法是利用訊號的相關性(correlation)進行計算,這個是我們後面將要介紹的方法;第三種方法是快速傅立葉變換(FFT),這是一個非常具有創造性和革命性的的方法,因為它大大提高了運算速度,使得傅立葉變換能夠在計算機中被廣泛應用,但這種演算法是根據複數形式的傅立葉變換來實現的,它把N個點的訊號分解成長度為N的頻域,這個跟我們現在所進行的實域DFT變換不一樣,而且這種方法也較難理解,這裡我們先不去理解,等先理解了複數DFT後,再來看一下FFT。有一點很重要,那就是這三種方法所得的變換結果是一樣的,經過實踐證明,當頻域長度為32時,利用相關性方法進行計算效率最好,否則FFT演算法效率較高。現在就讓我們來看一下相關性演算法。
 
利用第一種方法、訊號的相關性(correlation)可以從噪聲背景中檢測出已知的訊號,我們也可以利用這個方法檢測訊號波中是否含有某個頻率的訊號波:把一個待檢測訊號波乘以另一個訊號波,得到一個新的訊號波,再把這個新的訊號波所有的點進行相加,從相加的結果就可以判斷出這兩個訊號的相似程度。如下圖:

        上面a和 b兩個圖是待檢測訊號波,圖a很明顯可以看出是個3個週期的正弦訊號波,圖b的訊號波則看不出是否含有正弦或餘弦訊號,圖c和d都是個3個週期的正弦訊號波,圖e和f分別是a、b兩圖跟c、d兩圖相乘後的結果,圖e所有點的平均值是0.5,說明訊號a含有振幅為1的正弦訊號c,但圖f所有點的平均值是0,則說明訊號b不含有訊號d。這個就是通過訊號相關性來檢測是否含有某個訊號的方法。
 
       第二種方法:相應地,我也可以通過把輸入訊號和每一種頻率的正餘弦訊號進行相乘(關聯操作),從而得到原始訊號與每種頻率的關聯程度(即總和大小),這個結果便是我們所要的傅立葉變換結果,下面兩個等式便是我們所要的計算方法:

       第二個式子中加了個負號,是為了保持複數形式的一致,前面我們知道在計算時又加了個負號,所以這只是個形式的問題,並沒有實際意義,你也可以把負號去掉,並在計算時也不加負號。

       這裡有一點必須明白一個正交的概念:兩個函式相乘,如果結果中的每個點的總和為0,則可認為這兩個函式為正交函式。要確保關聯性演算法是正確的,則必須使得跟原始訊號相乘的訊號的函式形式是正交的,我們知道所有的正弦或餘弦函式是正交的,這一點我們可以通過簡單的高數知識就可以證明它,所以我們可以通過關聯的方法把原始訊號分離出正餘弦訊號。當然,其它的正交函式也是存在的,如:方波、三角波等形式的脈衝訊號,所以原始訊號也可被分解成這些訊號,但這只是說可以這樣做,卻是沒有用的。
       下面是實域傅立葉變換的BASIC語言程式碼:


 
       到此為止,我們對傅立葉變換便有了感性的認識了吧。但要記住,這只是在實域上的離散傅立葉變換,其中雖然也用到了複數的形式,但那只是個替代的形式,並無實際意義,現實中一般使用的是複數形式的離散傅立葉變換,且快速傅立葉變換是根據複數離散傅立葉變換來設計演算法的,在後面我們先來複習一下有關複數的內容,然後再在理解實域離散傅立葉變換的基礎上來理解複數形式的離散傅立葉變換。

第三章、複數

        複數擴充套件了我們一般所能理解的數的概念,複數包含了實數和虛數兩部分,利用複數的形式可以把由兩個變量表示的表示式變成由一個變數(復變數)來表達,使得處理起來更加自然和方便。
        我們知道傅立葉變換的結果是由兩部分組成的,使用複數形式可以縮短變換表示式,使得我們可以單獨處理一個變數(這個在後面的描述中我們就可以更加確切地知道),而且快速傅立葉變換正是基於複數形式的,所以幾乎所有描述的傅立葉變換形式都是複數的形式。
       但是複數的概念超過了我們日常生活中所能理解的概念,要理解複數是較難的,所以我們在理解複數傅立葉變換之前,先來專門複習一下有關複數的知識,這對後面的理解非常重要。
 
一、 複數的提出
 
      在此,先讓我們看一個物理實驗:把一個球從某點向上丟擲,然後根據初速度和時間來計算球所在高度,這個方法可以根據下面的式子計算得出:

其中h表示高度,g表示重力加速度(9.8m/s2),v表示初速度,t表示時間。現在反過來,假如知道了高度,要求計算到這個高度所需要的時間,這時我們又可以通過下式來計算:

        

多謝JERRY_PRI提出

    1、根據公式h=-(gt2/2)+Vt(gt後面的2表示t的平方),我們可以討論最終情況,也就是說小球運動到最高點時,v=gt,所以,可以得到t=sqt(2h/g)
且在您給的公式中,根號下為1-(2h)/g,化成分數形式為(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都表示平方),那麼也就是說當v2<2hg時會產生複數,但是如果從實際的v2是不可能小於2hg的,所以我感覺複數不能從實際出發去推到,只能從抽象的角度說明一下。

)      
      經過計算我們可以知道,當高度是3米時,有兩個時間點到達該高度:球向上運動時的時間是0.38秒,球向下運動時的時間是1.62秒。但是如果高度等於10時,結果又是什麼呢?根據上面的式子可以發現存在對負數進行開平方運算,我們知道這肯定是不現實的。
      第一次使用這個不一般的式子的人是義大利數學家Girolamo Cardano(1501-1576),兩個世紀後,德國偉大數學家Carl Friedrich Gause(1777-1855)提出了複數的概念,為後來的應用鋪平了道路,他對複數進行這樣表示:複數由實數(real)和虛數(imaginary)兩部分組成,虛數中的根號負1用i來表示(在這裡我們用j來表示,因為i在電力學中表示電流的意思)。
 
       我們可以把橫座標表示成實數,縱座標表示成虛數,則座標中的每個點的向量就可以用複數來表示,如下圖:
               
       上圖中的ABC三個向量可以表示成如下的式子:
 
            A = 2 + 6j
            B = -4 – 1.5j
            C = 3 – 7j
 
       這樣子來表達方便之處在於運用一個符號就能把兩個原來難以聯絡起來的數組合起來了,不方便的是我們要分辨哪個是實數和哪個是虛數,我們一般是用Re( )和Im( )來表示實數和虛數兩部分,如:
 
            Re A = 2      Im A = 6
            Re B = -4     Im B = -1.5
            Re C = 3      Im C = -7 
 
       複數之間也可以進行加減乘除運算:
            
  
       這裡有個特殊的地方是j2等於-1,上面第四個式子的計算方法是把分子和分母同時乘以c – dj,這樣就可消去分母中的j了。
 
       複數也符合代數運算中的交換律、結合律、分配律:
 
              A B = B A
              (A + B) + C = A + (B + C)
              A(B + C) = AB + AC
 
 
二、 複數的極座標表示形式
 
       前面提到的是運用直角座標來表示複數,其實更為普遍應用的是極座標的表示方法,如下圖:


              
       上圖中的M即是數量積(magnitude),表示從原點到座標點的距離,θ是相位角(phase angle),表示從X軸正方向到某個向量的夾角,下面四個式子是計算方法:
                     
     我們還可以通過下面的式子進行極座標到直角座標的轉換:
 
             a + jb = M (cosθ + j sinθ)


     上面這個等式中左邊是直角座標表示式,右邊是極座標表示式。
 
     還有一個更為重要的等式——尤拉等式(尤拉,瑞士的著名數學家,Leonhard Euler,1707-1783):
             ejx = cos x + j sin x 
 
     這個等式可以從下面的級數變換中得到證明:

      上面中右邊的兩個式子分別是cos(x)和sin(x)的泰勒(Taylor)級數。
 

      這樣子我們又可以把複數的表示式表示成指數的形式了:
 
             a + jb = M ejθ (這便是複數的兩個表示式)
 
      指數形式是數字訊號處理中數學方法的支柱,也許是因為用指數形式進行復數的乘除運算極為簡單的緣故吧:

              
三、複數是數學分析中的一個工具
 
       為什麼要使用複數呢?其實它只是個工具而已,就如釘子和錘子的關係,複數就象那錘子,作為一種使用的工具。我們把要解決的問題表達成複數的形式(因為有些問題用複數的形式進行運算更加方便),然後對複數進行運算,最後再轉換回來得到我們所需要的結果。
 
       有兩種方法使用複數,一種是用複數進行簡單的替換,如前面所說的向量表示式方法和前一節中我們所討論的實域DFT,另一種是更高階的方法:數學等價(mathematical equivalence),複數形式的傅立葉變換用的便是數學等價的方法,但在這裡我們先不討論這種方法,這裡我們先來看一下用複數進行替換中的問題。
 
       用複數進行替換的基本思想是:把所要分析的物理問題轉換成複數的形式,其中只是簡單地新增一個複數的符號j,當返回到原來的物理問題時,則只是把符號j去掉就可以了。
 
       有一點要明白的是並不是所有問題都可以用複數來表示,必須看用複數進行分析是否適用,有個例子可以看出用複數來替換原來問題的表達方式明顯是謬誤的:假設一箱的蘋果是5美元,一箱的桔子是10美元,於是我們把它表示成 5 + 10j,有一個星期你買了6箱蘋果和2箱桔子,我們又把它表示成6 + 2j,最後計算總共花的錢是(5 + 10j)(6 + 2j) = 10 + 70j,結果是買蘋果花了10美元的,買桔子花了70美元,這樣的結果明顯是錯了,所以複數的形式不適合運用於對這種問題的解決。
 
四、用複數來表示正餘弦函式表示式
 
       對於象M cos (ωt + φ)和A cos(ωt ) + B sin(ωt )表示式,用複數來表示,可以變得非常簡潔,對於直角座標形式可以按如下形式進行轉換:
       
       上式中餘弦幅值A經變換生成a,正弦幅值B的相反數經變換生成b:A <=> a,B<=> -b,但要注意的是,這不是個等式,只是個替換形式而已。
 
       對於極座標形式可以按如下形式進行轉換:
       
      
       上式中,M <=> M,θ<=>φ。
   這裡虛數部分採用負數的形式主要是為了跟複數傅立葉變換表示式保持一致,對於這種替換的方法來表示正餘弦,符號的變換沒有什麼好處,但替換時總會被改變掉符號以跟更高階的等價變換保持形式上的一致。
 
        在離散訊號處理中,運用複數形式來表示正餘弦波是個常用的技術,這是因為利用複數進行各種運算得到的結果跟原來的正餘弦運算結果是一致的,但是,我們要小心使用複數操作,如加、減、乘、除,有些操作是不能用的,如兩個正弦訊號相加,採用複數形式進行相加,得到的結果跟替換前的直接相加的結果是一樣的,但是如果兩個正弦訊號相乘,則採用複數形式來相乘結果是不一樣的。幸運的是,我們已嚴格定義了正餘弦複數形式的運算操作條件:
1、參加運算的所有正餘弦的頻率必須是一樣的;
2、運算操作必須是線性的,如兩個正弦訊號可以進行相加減,但不能進行乘除,象訊號的放大、衰減、高低通濾波等系統都是線性的,象平方、縮短、取限等則不是線性的。要記住的是卷積和傅立葉分析也只有線性操作才可以進行。
   
       下圖是一個相量變換(我們把正弦或餘弦波變成複數的形式稱為相量變換,Phasor transform)的例子,一個連續訊號波經過一個線性處理系統生成另一個訊號波,從計算過程我們可以看出採用複數的形式使得計算變化十分的簡潔:
 
    在第二章中我們描述的實數形式傅立葉變換也是一種替換形式的複數變換,但要注意的是那還不是複數傅立葉變換,只是一種代替方式而已。下一章、即,第四章,我們就會知道複數傅立葉變換是一種更高階的變換,而不是這種簡單的替換形式。 

第四章、複數形式離散傅立葉變換

    複數形式的離散傅立葉變換非常巧妙地運用了複數的方法,使得傅立葉變換變換更加自然和簡潔,它並不是只是簡單地運用替換的方法來運用複數,而是完全從複數的角度來分析問題,這一點跟實數DFT是完全不一樣的。
 
一、  把正餘弦函式表示成複數的形式

    通過尤拉等式可以把正餘弦函式表示成複數的形式:
 
                    cos( x ) = 1/2 e j(-x) + 1/2 ejx 
                    sin( x ) = j (1/2 e j(-x) - 1/2 ejx)

    從這個等式可以看出,如果把正餘弦函式表示成複數後,它們變成了由正負頻率組成的正餘弦波,相反地,一個由正負頻率組成的正餘弦波,可以通過複數的形式來表示。
 
    我們知道,在實數傅立葉變換中,它的頻譜是0 ~ π(0 ~ N/2),但無法表示-π~ 0的頻譜,可以預見,如果把正餘弦表示成複數形式,則能夠把負頻率包含進來。
 
二、  把變換前後的變數都看成複數的形式
 
    複數形式傅立葉變換把原始訊號x[n]當成是一個用複數來表示的訊號,其中實數部分表示原始訊號值,虛數部分為0,變換結果X[k]也是個複數的形式,但這裡的虛數部分是有值的。
    在這裡要用複數的觀點來看原始訊號,是理解複數形式傅立葉變換的關鍵(如果有學過複變函式則可能更好理解,即把x[n]看成是一個複數變數,然後象對待實數那樣對這個複數變數進行相同的變換)。
 
三、  對複數進行相關性演算法(正向傅立葉變換)
 
     從實數傅立葉變換中可以知道,我們可以通過原始訊號乘以一個正交函式形式的訊號,然後進行求總和,最後就能得到這個原始訊號所包含的正交函式訊號的分量。

     現在我們的原始訊號變成了複數,我們要得到的當然是複數的訊號分量,我們是不是可以把它乘以一個複數形式的正交函式呢?答案是肯定的,正餘弦函式都是正交函式,變成如下形式的複數後,仍舊還是正交函式(這個從正交函式的定義可以很容易得到證明):

                   cos x + j sin x, cos x – j sin x,……
 
     這裡我們採用上面的第二個式子進行相關性求和,為什麼用第二個式子呢?,我們在後面會知道,正弦函式在虛數中變換後得到的是負的正弦函式,這裡我們再加上一個負號,使得最後的得到的是正的正弦波,根據這個於是我們很容易就可以得到了複數形式的DFT正向變換等式

       
     這個式子很容易可以得到尤拉變換式子:
       
 
     其實我們是為了表達上的方便才用到尤拉變換式,在解決問題時我們還是較多地用到正餘弦表示式。
 
       對於上面的等式,我們要清楚如下幾個方面(也是區別於實數DFT的地方):
1、X[k]、x[n]都是複數,但x[n]的虛數部分都是由0組成的,實數部分表示原始訊號;
2、k的取值範圍是0 ~ N-1 (也可以表達成0 ~ 2π),其中0 ~ N/2(或0 ~ π)是正頻部分,

N/2 ~ N-1(π~ 2π)是負頻部分,由於正餘弦函式的對稱性,所以我們把 –π~ 0表示成π~ 2π,這是出於計算上方便的考慮。
3、其中的j是一個不可分離的組成部分,就象一個等式中的變數一樣,不能隨便去掉,去掉之後意義就完全不一樣了,但我們知道在實數DFT中,j只是個符號而已,把j去掉,整個等式的意義不變;
4、下圖是個連續訊號的頻譜,但離散頻譜也是與此類似的,所以不影響我們對問題的分析:
             

 
     上面的頻譜圖把負頻率放到了左邊,是為了迎合我們的思維習慣,但在實際實

現中我們一般是把它移到正的頻譜後面的。

     從上圖可以看出,時域中的正餘弦波(用來組成原始訊號的正餘弦波)在複數DFT的頻譜中被分成了正、負頻率的兩個組成部分,基於此等式中前面的比例係數是1/N(或1/2π),而不是2/N,這是因為現在把頻譜延伸到了2π,但把正負兩個頻率相加即又得到了2/N,又還原到了實數DFT的形式,這個在後面的描述中可以更清楚地看到。

     由於複數DFT生成的是一個完整的頻譜,原始訊號中的每一個點都是由正、負兩個頻率組合而成的,所以頻譜中每一個點的頻寬是一樣的,都是1/N,相對實數DFT,兩端頻寬比其它點的頻寬少了一半;複數DFT的頻譜特徵具有周期性:-N/2 ~ 0與N/2 ~ N-1是一樣的,實域頻譜呈偶對稱性(表示餘弦波頻譜),虛域頻譜呈奇對稱性(表示正弦波頻譜)。
 

四、  逆向傅立葉變換
 
     假設我們已經得到了複數形式的頻譜X[k],現在要把它還原到複數形式的原始訊號x[n],當然應該是把X[k]乘以一個複數,然後再進行求和,最後得到原始訊號x[n],這個跟X[k]相乘的複數首先讓我們想到的應該是上面進行相關性計算的複數:

                     cos(2πkn/N) – j si(2πkn/N),
 
     但其中的負號其實是為了使得進行逆向傅立葉變換時把正弦函式變為正的符號,因為虛數j的運算特殊性,使得原來應該是正的正弦函式變為了負的正弦函式(我們從後面的推導會看到這一點),所以這裡的負號只是為了糾正符號的作用,在進行逆向DFT時,我們可以把負號去掉,於是我們便得到了這樣的逆向DFT變換等式

                     x[n] = X[k] (cos(2πkn/N) + j sin(2πkn/N))

     我們現在來分析這個式子,會發現這個式其實跟實數傅立葉變換是可以得到一樣結果的。我們先把X[k]變換一下:

                     X[k] = Re X[k] + j Im X[k]


      這樣我們就可以對x[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)

      這時我們就把原來的等式分成了兩個部分,第一個部分是跟實域中的頻譜相乘,第二個部分是跟虛域中的頻譜相乘,根據頻譜圖我們可以知道,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表示負頻率,為了表達方便我們把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))

      根據偶對稱性和三角函式的性質,把上式化簡得到:
 
                    Re X[k] (cos(2πkn/N) + j sin(2πkn/N)) + Re X[ k] (cos( 2πkn/N) - j sin( 2πkn/N))

      這個式子最後的結果是:

                    2 Re X[ k] cos(2πkn/N)。
      
      再考慮到求Re X[ k]等式中有個比例係數1/N,把1/N乘以2,這樣的結果不就是跟實數DFT中的式子一樣了嗎?
 
      對於(2)式,用同樣的方法,我們也可以得到這樣的結果:

                    -2 Im X[k] sin(2πkn/N)

       注意上式前面多了個負符號,這是由於虛數變換的特殊性造成的,當然我們肯定不能把負符號的正弦函式跟餘弦來相加,還好,我們前面是用cos(2πkn/N) – j sin(2πkn/N)進行相關性計算,得到的Im X[k]中有個負的符號,這樣最後的結果中正弦函式就沒有負的符號了,這就是為什麼在進行相關性計算時虛數部分要用到負符號的原因(我覺得這也許是複數形式DFT美中不足的地方,讓人有一種拼湊的感覺)。
 
       從上面的分析中可以看出,實數傅立葉變換跟複數傅立葉變換,在進行逆變換時得到的結果是一樣的,只不過是殊途同歸吧。本文完。(July、dznlong)