Wave (ℝ)

日往則月來,月往則日來,日月相推而明生焉。《易傳》

振動、振盪

這個世界天天都在振動。地面、空氣、海水、機械、人體等等,都是不斷振動。

振動、振盪是物理學名詞,振動(Vibration)是來回運動,振盪(Oscillation)是來回變化。

震動、震盪是自古以來就有的詞彙。

振動可以用函數表示

每個時間點的振動高低,可以描繪成函數圖形,橫向是時間軸,縱向是每個時刻的振動高低位置。

平穩的振動

最平穩的振動,就是高中物理教的簡諧運動,是等速圓周運動在座標軸上的投影,呈現sin函數。sin函數和cos函數長相一樣,只是起點不同而已。

舉例來說,敲打音叉產生的振動,就非常接近平穩的振動。

振動的快慢:頻率

單位時間振動的次數,稱作「頻率Frequency」。

一秒振動的次數,單位是赫茲Hz。

人類能感知頻率:耳朵能聽到20Hz至20000Hz的空氣振動,低頻低沉、高頻尖銳;眼睛能看到4×10^14Hz至8×10^14Hz的電磁振盪,低頻至高頻分別呈現紅橙黃綠藍靛紫。

振動的高低:振幅

振動的最高(低)距離,稱作「振幅Amplitude」。

人類能感知振幅,受頻率大小影響。就聽覺而言,振幅高則大聲、振幅低則小聲;就視覺而言,振幅高則亮、振幅低則暗。

題外話,人類對於頻率與振幅的區分能力,大略等於取log。

振動的起點:相位

振動的起點位置,稱作「相位Phase」。

注意到,相位是圓周運動cos函數的角度,而不是振動高低位置。物理學家喜歡用相位。

人類幾乎分辨不出相位的差異。

振動有疊加效果

現實世界當中,多個振動時常融合成一個振動,等於各個振動高低位置相加。相同方向則增益、相反方向則抵銷。

寫成數學式子,就是多個函數相加。

振動有傳遞效果:波

一個粒子振動,就會牽引隔壁粒子振動,一傳十、十傳百。宏觀之下,形成「波Wave」,亦譯「波動Wave」。

觀察任意一個粒子,都是在振動。傳遞速度取決於粒子之間的作用力、粒子的質量。作用力強、質量小,則傳遞速度快。

均勻分布的粒子之中,某個粒子振動所產生的波,剛好也呈現sin函數,英文稱作sine wave或者sinusoid。

水的高低起伏,就是水波。空氣的疏密,就是聲波。地的高低起伏與左右晃動,就是地震波。電場與磁場的交互作用,就是電磁波。光波經實驗證明是電磁波。原子的振動,也許是熱。有人覺得氣功也許是波,就叫做氣功波。

Wave (ℂ)

有物混成,先天地生,寂兮寥兮,獨立而不改,周行而不殆,可以為天下母。《老子》

Complex Number

快速複習一下複數吧。實數,再額外考慮i = √-1,就是複數。例如2 + 3i、(1 - √2) + (1/3)i、1 / (-2i - 4)、∛i、sin(i)。

只要是複數,都可以重新整理成左邊實數不乘i、右邊實數有乘i,兩個部分相加的格式,證明省略之;其中不乘i的部分叫做實部(real part),有乘i的部分叫做虛部(imaginary part)。例如1 / (-2i - 4)可以重新整理成-0.2 + 0.1i,其中實部是-0.2,虛部是0.1i。

複數亦可以用圖形表示。

複數平面、二維平面、極座標平面是不同的事情,不要搞混了。

兩個複數相加,就是實部加實部、虛部加虛部。在複數平面上,外觀宛如向量相加。

兩個複數相乘,就是實乘實、虛乘虛、實乘虛、虛乘實,再累加這四個乘積。在複數平面上,外觀宛如長度相乘、角度相加。

一個複數可以重新表示成一個長度與一個角度,叫做極座標表示法。長度可以用畢氏定理求得,角度可以用arctan函數求得。

一個長度與一個角度也可以還原成一個複數。實部可以用cos函數求得,虛部可以用sin函數求得。

附帶一提,長度也有人稱作強度(magnitude),角度也有人稱作相位(phase)。

Euler's Formula

強者歐拉發現這世界上有一個神奇數字e,e的純虛數次方竟然在複數平面上繞圈兒。這真是一個超乎常理的發現!

寫成數學公式是:e = cosθ + i * sinθ,複數的長度是常數1,複數的角度是變數θ。等式右邊,是將長度1與角度θ,還原成一個複數cosθ + i * sinθ,外觀很複雜但是本質很簡單。

有了歐拉公式,一個複數也可以重新表示成e的次方、另乘上倍率。次方值即是角度乘i,倍率即是長度。

歐拉公式,定量增加θ,在複數平面上,外觀宛如「等速圓周運動」,逆時針繞圈;只看實部或者只看虛部,外觀宛如「簡諧運動」,先上後下。

繞360°是一圈,剛好回到+1;繞180°是半圈,剛好是-1。因此有了e + 1 = 0這條著名等式,π就是180°。

e運算簡單,考慮長度與角度即可。e性質優美,每轉90°剛好是±1與±i。也許你會漸漸愛上它。

這個e,大約是2.71828183,是自然對數的底數e,是1/x積分後所出現的e。離題了。

Wave

因為e跟波很像,所以可以用e來表示波。實數系統的波,推廣成為複數系統的波。

複數系統的e波,欲降格為實數系統的sin波、cos波,有兩種方式。觀察e = cosθ + i * sinθ這道式子:

一、取實部得到cos波、取虛部得到sin波。即是俯瞰和側視。

二、用e與e-iθ,相加除以二得到cos波,相減除以二得到sin波。

Fourier Transform

Fourier Transform輸入輸出

傅立葉轉換是一對一函數,輸入輸出都是一串複數,可以是連續函數或者離散數列,各有對應名稱。混淆視聽罷了。

輸入 輸出 名稱
連續 連續 Fourier Transform
連續 離散 Fourier Series
離散 連續 Discrete-time Fourier Transform
離散 離散 Discrete Fourier Transform

連續到連續的傅立葉轉換,輸入和輸出都是一個ℝ ⇨ ℂ函數。(ℝ ⇨ ℂ函數很難畫,故圖例為實數函數。)

離散到離散的傅立葉轉換,輸入和輸出都是一串複數數列。離散版本是連續版本的特例:輸入輸出是週期函數、等距取樣。(複數數列很難畫,故圖例為實數數列。)

電腦做運算,數值皆離散。本篇只介紹離散版本。

Fourier Transform物理意義

N個複數波,頻率是0倍到N-1倍,分別是ei⋅(2π/N)⋅0⋅θ、ei⋅(2π/N)⋅1⋅θ、……、ei⋅(2π/N)⋅(N-1)⋅θ。寫成代數是ei⋅(2π/N)⋅n⋅θ。(複數波很難畫,故圖例為實數波。)

輸入數列與一個波,N個對應位置點對點相除、再求總和,得到一個輸出數值。可以簡單想做:輸入數列除以波,求比例。

(學過線性代數的讀者,也可以想做點積、投影。)

輸入數列,分別除以N個波,得到N個輸出數值,形成輸出數列。這就是傅立葉轉換。

正向傅立葉轉換,是把一個複雜的波,拆解成N個平穩的波,頻率是0倍到N-1倍,強度與相位是N個輸出數值的強度與相位。

逆向傅立葉轉換,是把N個平穩的波,頻率是0倍到N-1倍,分別乘上強度、加上相位,再疊加成一個複雜的波。

Fourier Transform數學公式

傅立葉轉換
       N-1
y[n] =  ∑ { x[k] ÷ ei⋅(2π/N)⋅n⋅k } ÷ sqrt(N)
       k=0
       N-1
     =  ∑ { x[k] ⋅ e-i⋅(2π/N)⋅n⋅k } ÷ sqrt(N)
       k=0

傅立葉轉換的反函數,稱作逆向傅立葉轉換
       N-1
x[n] =  ∑ { y[k] ⋅ ei⋅(2π/N)⋅n⋅k } ÷ sqrt(N)
       k=0

為了加快計算速度,正向傅立葉轉換經常改成不除以sqrt(N),逆向傅立葉轉換經常改成多除以sqrt(N)。

       N-1
y[n] =  ∑ { x[k] ÷ ei⋅(2π/N)⋅n⋅k }
       k=0
       N-1
x[n] =  ∑ { y[k] ⋅ ei⋅(2π/N)⋅n⋅k } ÷ N
       k=0

傅立葉轉換是線性變換,其矩陣恰是正規正交矩陣。

ω = ei⋅2π/N

[y0  ]   [ ω-0⋅0     ω-0⋅1     ω-0⋅2     ..  ω-0⋅(N-1)     ]  [x0  ]
[y1  ]   [ ω-1⋅0     ω-1⋅1     ω-1⋅2     ..  ω-1⋅(N-1)     ]  [x1  ]
[y2  ] = [ ω-2⋅0     ω-2⋅1     ω-2⋅2     ..  ω-2⋅(N-1)     ]  [x2  ]
[:   ]   [ :         :         :             :            ]  [:   ]
[yN-1]   [ ω-(N-1)⋅0 ω-(N-1)⋅1 ω-(N-1)⋅2 ..  ω-(N-1)⋅(N-1) ]  [xN-1]

[x0  ]       [ ω0⋅0     ω0⋅1     ω0⋅2     ..  ω0⋅(N-1)     ]  [y0  ]
[x1  ]    1  [ ω1⋅0     ω1⋅1     ω1⋅2     ..  ω1⋅(N-1)     ]  [y1  ]
[x2  ] = --- [ ω2⋅0     ω2⋅1     ω2⋅2     ..  ω2⋅(N-1)     ]  [y2  ]
[:   ]    N  [ :        :        :            :           ]  [:   ]
[xN-1]       [ ω(N-1)⋅0 ω(N-1)⋅1 ω(N-1)⋅2 ..  ω(N-1)⋅(N-1) ]  [yN-1]

複數波,變成離散數列,可以畫成這樣子:

傅立葉轉換的矩陣,可以畫成這樣子:

演算法(公式解)

依照公式實作,時間複雜度是O(N^2)。

演算法(Cooley-Tukey Algorithm)

時間複雜度優於O(N^2)的傅立葉轉換演算法,老人家就直接稱作「快速傅立葉轉換Fast Fourier Transform, FFT」。

這裡介紹最經典的快速傅立葉轉換。公式的偶數項與奇數項分開整理,採用Dynamic Programming,時間複雜度是O(NlogN)。由於必須剛好對半分,所以N必須剛好是2的次方。當N不是2的次方,可在輸入數列末端補零,理由容後介紹。

逆向快速傅立葉轉換的計算原理也是一樣的,此處省略。

【所有文獻皆歸類為Divide and Conquer,不太準確。】

                           FFT
(x0 x1 x2 x3 x4 x5 x6 x7) ----> (y0 y1 y2 y3 y4 y5 y6 y7)

N = 8, ω = e-i⋅2π/N   注意到ω放入了負號,讓下面的數學式子比較簡潔

y0 = x0ω0 + x1ω0 + x2ω0 + x3ω0 + x4ω0 + x5ω0 + x6ω0 + x7ω0
   = (x0ω0 + x2ω0 + x4ω0 + x6ω0) + (x1ω0 + x3ω0 + x5ω0 + x7ω0)
   = (x0ω0 + x2ω0 + x4ω0 + x6ω0) + ω0 ⋅ (x1ω0 + x3ω0 + x5ω0 + x7ω0)
   = (x0 x2 x4 x6)轉換結果第0項 + ω0 ⋅ (x1 x3 x5 x7)轉換結果第0項
   = y偶0 + ω0 ⋅ y奇0

y1 = x0ω0 + x1ω1 + x2ω2 + x3ω3 + x4ω4 + x5ω5 + x6ω6 + x7ω7
   = (x0ω0 + x2ω2 + x4ω4 + x6ω6) + (x1ω1 + x3ω3 + x5ω5 + x7ω7)
   = (x0ω0 + x2ω2 + x4ω4 + x6ω6) + ω1 ⋅ (x1ω0 + x3ω2 + x5ω4 + x7ω6)
   = (x0υ0 + x2υ1 + x4υ2 + x6υ3) + ω1 ⋅ (x1υ0 + x3υ1 + x5υ2 + x7υ3)
   = (x0 x2 x4 x6)轉換結果第1項 + ω1 ⋅ (x1 x3 x5 x7)轉換結果第2項
   = y偶1 + ω1 ⋅ y奇1

y2 = x0ω0 + x1ω2 + x2ω4 + x3ω6 + x4ω8 + x5ω10 + x6ω12 + x7ω14
   = (x0ω0 + x2ω4 + x4ω8 + x6ω12) + (x1ω2 + x3ω6 + x5ω10 + x7ω14)
   = (x0ω0 + x2ω4 + x4ω8 + x6ω12) + ω2 ⋅ (x1ω0 + x3ω4 + x5ω8 + x7ω12)
   = (x0υ0 + x2υ2 + x4υ4 + x6υ6 ) + ω2 ⋅ (x1υ0 + x3υ2 + x5υ4 + x7υ6 )
   = (x0 x2 x4 x6)轉換結果第2項 + ω2 ⋅ (x1 x3 x5 x7)轉換結果第2項
   = y偶2 + ω2 ⋅ y奇2

y3 = x0ω0 + x1ω3 + x2ω6 + x3ω9 + x4ω12 + x5ω15 + x6ω18 + x7ω21
   = (x0ω0 + x2ω6 + x4ω12 + x6ω18) + (x1ω3 + x3ω9 + x5ω15 + x7ω21)
   = (x0ω0 + x2ω6 + x4ω12 + x6ω18) + ω3 ⋅ (x1ω0 + x3ω6 + x5ω12 + x7ω18)
   = (x0υ0 + x2υ3 + x4υ6  + x6υ9 ) + ω3 ⋅ (x1υ0 + x3υ3 + x5υ6  + x7υ9 )
   = (x0 x2 x4 x6)轉換結果第3項 + ω3 ⋅ (x1 x3 x5 x7)轉換結果第3項
   = y偶3 + ω3 ⋅ y奇3

注意到 ω8 = 1
y4 = x0ω0 + x1ω4 + x2ω8 + x3ω12 + x4ω16 + x5ω20 + x6ω24 + x7ω28
   = (x0ω0 + x2ω8 + x4ω16 + x6ω24) + (x1ω4 + x3ω12 + x5ω20 + x7ω28)
   = (x0ω0 + x2ω8 + x4ω16 + x6ω24) + ω4 ⋅ (x1ω0 + x3ω8 + x5ω16 + x7ω24)
   = (x0ω0 + x2ω0 + x4ω0  + x6ω0 ) + ω4 ⋅ (x1ω0 + x3ω0 + x5ω0  + x7ω0 )
   = (x0 x2 x4 x6)轉換結果第0項 + ω4 ⋅ (x1 x3 x5 x7)轉換結果第0項
   = y偶0 + ω4 ⋅ y奇0

注意到 ω8 = 1
y5 = x0ω0 + x1ω5 + x2ω10 + x3ω15 + x4ω20 + x5ω20 + x6ω25 + x7ω30
   = (x0ω0 + x2ω10 + x4ω20 + x6ω30) + (x1ω5 + x3ω10 + x5ω25 + x7ω30)
   = (x0ω0 + x2ω10 + x4ω20 + x6ω30) + ω5 ⋅ (x1ω0 + x3ω10 + x5ω20 + x7ω30)
   = (x0ω0 + x2ω2  + x4ω4  + x6ω6 ) + ω5 ⋅ (x1ω0 + x3ω4  + x5ω6  + x7ω8 )
   = (x0υ0 + x2υ1  + x4υ2  + x6υ3 ) + ω5 ⋅ (x1υ0 + x3υ2  + x5υ3  + x7υ4 )
   = (x0 x2 x4 x6)轉換結果第1項 + ω5 ⋅ (x1 x3 x5 x7)轉換結果第1項
   = y偶1 + ω5 ⋅ y奇1

y6與y7當做作業,請讀者自己推導吧!
y0 = y偶0 + y奇0 ⋅ ω0
y1 = y偶1 + y奇1 ⋅ ω1
y2 = y偶2 + y奇2 ⋅ ω2
y3 = y偶3 + y奇3 ⋅ ω3
y4 = y偶0 + y奇0 ⋅ ω4
y5 = y偶1 + y奇1 ⋅ ω5
y6 = y偶2 + y奇2 ⋅ ω6
y7 = y偶3 + y奇3 ⋅ ω7

觀察DP的遞推過程,偶數項與奇數項分開處理之後,索引值不是連續的,不好取值。預先重新排列陣列元素,符合遞推過程,就可以減少cache miss,還可以重複使用記憶體、節省空間。

如何重新排列呢?有趣的是,索引值的二進位表示法,高低位數顛倒之後,恰是正確結果!

重新排列的時間複雜度為O(NlogN)。似乎可以加速為O(N)。

2D Fourier Transform

傅立葉轉換可以推廣到高維度。

離散版本的二維傅立葉轉換,輸入輸出都是一個N×N複數方陣。輸入方陣,分別除以N×N種二維複數波,得到N×N個輸出數值,形成輸出方陣。(由於二維複數波很難畫,以下改畫二維實數波。)

依照公式實作,時間複雜度為O(N^4)。快速的演算法,是每一橫條各自傅立葉轉換,然後每一直條各自傅立葉轉換,時間複雜度為O(NNlogN + NNlogN) = O(NNlogN)。

Hartley Transform

哈特利轉換是一對一函數,輸入和輸出都是一串實數。

哈特利轉換與傅立葉轉換如出一轍,只少了虛數i而已。

傅立葉轉換:
    2πnk           2πnk    -i2πnk/N
cos ———— - i ⋅ sin ———— = e
     N              N
哈特利轉換:
    2πnk           2πnk       2πnk
cos ———— +     sin ———— = cas ————
     N              N          N  
另一個哈特利轉換,比較沒人用:
    2πnk           2πnk
cos ———— -     sin ————
     N              N  
傅立葉轉換:
       N-1
y[n] =  ∑ { x[k] ÷ ei⋅(2π/N)⋅n⋅k } ÷ sqrt(N)
       k=0
哈特利轉換:
       N-1
y[n] =  ∑ { x[k] ⋅ cas((2π/N)⋅n⋅k) } ÷ sqrt(N)
       k=0

哈特利轉換的輸出,可以調整成傅立葉轉換的輸出,O(N):

http://mathworld.wolfram.com/HartleyTransform.html

實數運算比複數運算還要簡單,所以哈特利轉換比傅立葉轉換還要快速。聲音處理、影像處理時,訊號都是實數、甚至是整數,剛好也符合哈特利轉換的輸入格式。因此一般都是套用哈特利轉換進行計算,再把結果調整成傅立葉轉換。

由於哈特利轉換與傅立葉轉換的公式幾乎相同,所以兩者的演算法也是一一對應。這裡介紹的也是運用Divide and Conquer的方法。不一樣的是奇數項的處理方式,提出常數的步驟變複雜了。

  N-1
   ∑ { x[k] ⋅ cas((2π/N)⋅n⋅k) }
  k=1,3,5,...
  N/2-1
=  ∑ { x[2k+1] ⋅ cas((2π/N)⋅n⋅(2k+1)) }
  k=0,1,2,...
  N/2-1
=  ∑ { x[2k+1] ⋅ ( cas((2π/N)⋅n⋅2k) ⋅ cos((2π/N)⋅n⋅1)
  k=0,1,2,...       + cas(-(2π/N)⋅n⋅2k) ⋅ sin((2π/N)⋅n⋅1) ) }
  N/2-1
=  ∑ { x[2k+1] ⋅ ( cas((2π/(N/2))⋅n⋅k) ⋅ cos((2π/N)⋅n⋅1)
  k=0,1,2,...       + cas(-(2π/(N/2))⋅n⋅k) ⋅ sin((2π/N)⋅n⋅1) ) }
  N/2-1
=  ∑ { x[2k+1] ⋅ cas( (2π/(N/2))⋅n⋅k) } ⋅ cos((2π/N)⋅n⋅1)
  k=0,1,2,...
  N/2-1
+  ∑ { x[2k+1] ⋅ cas(-(2π/(N/2))⋅n⋅k) } ⋅ sin((2π/N)⋅n⋅1)
  k=0,1,2,...

= y[n] ⋅ cos((2π/N)⋅n⋅1) + y[-n] ⋅ sin((2π/N)⋅n⋅1)
θ = 2π / N
y0 = y偶0 + y奇0 * cos0θ + y奇0 ⋅ sin0θ
y1 = y偶1 + y奇1 * cos1θ + y奇3 ⋅ sin1θ
y2 = y偶2 + y奇2 * cos2θ + y奇2 ⋅ sin2θ
y3 = y偶3 + y奇3 * cos3θ + y奇1 ⋅ sin3θ
y4 = y偶0 + y奇0 * cos4θ + y奇0 ⋅ sin4θ
y5 = y偶1 + y奇1 * cos5θ + y奇3 ⋅ sin5θ
y6 = y偶2 + y奇2 * cos6θ + y奇2 ⋅ sin6θ
y7 = y偶3 + y奇3 * cos7θ + y奇1 ⋅ sin7θ

下面是據我所知效率最好的實作程式碼:

http://home.iae.nl/users/mhx/fft.c

http://reocities.com/ResearchTriangle/8869/fft_summary.html

Fourier Transform的性質

Frequency Spectrum

傅立葉轉換,輸出數列有N個複數,可以畫成函數。一般不畫實部與虛部,而是畫長度與角度,具備物理意義。

這N個複數的長度(強度)畫成函數,稱為「強度頻譜」。

這N個複數的角度(相位)畫成函數,稱為「相位頻譜」。

兩者合稱為「頻譜」。

附帶一提,當輸入數列皆是實數,則輸出數列將共軛對稱:長度(強度)相等、角度(相位)負號。教科書為了讓圖片美觀,經常循環位移令中央為低頻、畫成折線圖、強度取log10。讀者要注意!

我們得以運用正向傅立葉轉換分解一個波,運用逆向傅立葉轉換合成一個波,運用頻譜解讀波的詳細內容。傅立葉轉換是一對一函數,一種波對應一種頻譜。頻譜的左側到右側,是低頻到高頻。

甚至可以把一個波實施正向傅立葉轉換,將低頻數值或者高頻數值改成零,再實施逆向傅立葉轉換,改造原本的波。這是十分常用的技巧。

頻譜是非常實用的分析工具。凡是學習科學的人,都有必要了解頻譜!各種物質的振動或振盪,皆可求得頻譜,發掘其特性。例如震譜是震波的頻譜,光譜是光波的頻譜,聲譜是聲波的頻譜。世間萬物皆有譜,應用無限廣泛。

解讀頻譜

範例:一條實數數列,16個數字,實施傅立葉轉換,得到頻譜。

起點是1,平穩振動1次,振幅為1,形成cos波:對應傅立葉轉換的一倍頻率波,頻譜第一點的強度是8、相位是0,其餘的強度和相位是0。

起點調成0,也就是相位調成-π/2:依然對應傅立葉轉換的一倍頻率波,強度依舊,相位是-π/2。

平穩振動調成2次:對應傅立葉轉換的兩倍頻率波,頻譜第二點的強度是8、相位是-π/2,其餘的強度和相位是0。

振幅調成2:強度變兩倍。

振動基準從0調成1:對應傅立葉轉換的零倍頻率波,其功效是數列總和,頻譜第零點的強度變16。

頻譜的缺點(一)

問題來了。平穩振動1.5次,頻譜如何?

你可能馬上聯想到「加權平均值」的概念,第一點和第二點有強度,強度各半。

但是事實並非如此。1.5倍,頻譜呈現「人」型,所有頻率皆有強度,漏得到處都是。這個現象稱作「spectral leakage」。

這個現象有兩種解讀:

一、1倍和2倍頻率波,疊加(加權平均值)之後,結果是1倍,不是1.5倍。更明確來說是最大公因數。

傅立葉轉換的0倍頻率波到n-1倍頻率波,皆無法組合出1.5倍。只好湊合各種頻率,盡量趨近1.5倍。

二、離散版本的傅立葉轉換,輸入輸出是循環數列。1.5倍,循環之後,其實不是平穩的振動,因而產生許多高頻波。

當振動次數不是整數次、頻率不是整數倍,那麼傅立葉轉換無法精準量測!這是重大缺點!

Window Function

然而數學家尚未發明更好的方式。當今主流仍是傅立葉轉換。

為了克服spectral leakage這個重大缺點,數學家想出了「窗函數」。

原本數列,乘上一個窗函數:中央高、兩端趨近零的數列。如此令原本數列左右兩端連續,抑制頻譜多餘強度。

窗函數非常多種,功效略有差異。請讀者自行研究。

Window Function的快速演算法

傅立葉轉換,有時候輸入稱作「時域Time Domain」、輸出稱作「頻域Frequency Domain」,呼應傅立葉轉換的功能:把波(時間軸)表示成頻譜(頻率軸)。

時域乘法等於頻域循環摺積,請參考本站文件「Convolution」。數列與窗函數相乘,等於數列與窗函數在頻域的循環摺積。

窗函數,多由cos波組成;窗函數在頻域,只有少數幾點有值。例如Hann窗,從時域轉頻域,只有三點有值。

因此,與其在時域套用窗函數,不如在頻域套用窗函數。過程非常簡單:每個值減去兩側的值(相位差不多是π),附帶權重。這揭露了窗函數的真正功效──頻譜之中,平者更平,尖者更尖。

最後額外補充一下。採用連續版本的傅立葉轉換,窗函數頻譜,外觀是一個尖峰。取abs和log,外觀是一個大圓丘(main lobe),附帶連綿小矮丘(side lobe)。很多資工系老師上課只教連續版本,但是我們根本不會用到連續版本!

F = FourierTransform[HannWindow[x], x, w] Plot[F, {w, 0, +70}, PlotRange->{-0.05,+0.2}, Axes->None] F = FourierTransform[HannWindow[x], x, w] Plot[F, {w, 0, +70}, PlotRange->{-0.001,+0.001}, Axes->None] F = Abs[FourierTransform[HannWindow[x], x, w]] LogPlot[F, {w, 0, +70}, Axes->None]

頻譜的缺點(二)

當波形不是完美的sin波,那麼傅立葉轉換無法精準量測!這是重大缺點!

目前無解。自己保重。

頻譜的缺點(三)

聲音波形經常疊加。舉例來說,兩個頻率不同的音叉,同時敲擊,耳膜感受到的振動,差不多就是兩個sin波相加。更明確來說是兩個sin波的加權平均值。

傅立葉轉換是線性變換。換句話說,輸入數列們的加權平均值,經過傅立葉轉換,等於輸出數列們的加權平均值;但是不等於頻譜們的加權平均值!

輸出數列是複數。複數加法是向量相加,複數倍率是向量伸縮。向量相加不等於長度相加、角度相加。(唯一例外:所有波都是整數次的平穩振動。因為頻譜幾乎都是零。)

多個波形疊加,不會正確反映於頻譜!這是重大缺點!

然而大家仍用頻譜分解頻率,無法可管。自己保重。

Sparse Fourier Transform

http://groups.csail.mit.edu/netmit/sFFT/

http://people.csail.mit.edu/indyk/fourier-gsip.pdf

只計算特定頻率的強度與相位。速度較快。

ListPlot[Table[Sin[x*2*Pi/16], {x, 0, 15}]] ListPlot[Abs[Fourier[Table[Sin[x*2*Pi/16], {x, 0, 15}]]], PlotRange->{0, 2}, Filling->Axis] ListPlot[Arg[Fourier[Table[Sin[x*2*1.5*Pi/16], {x, 0, 15}]]], PlotRange->{-4, +4}, Filling->Axis] ListPlot[Abs[Fourier[Table[HannWindow[(x-16)/32], {x, 0, 31}]]], PlotRange->{0, 2}, Filling->Axis] ListPlot[Arg[Fourier[Table[HannWindow[(x-16)/32], {x, 0, 31}]]], PlotRange->{-4, 4}, Filling->Axis] ListPlot[Table[HannWindow[(x-16)/32], {x, 0, 31}], PlotRange->{0, 1}, Filling->Axis, FillingStyle->Red, PlotStyle->Red, Axes->None] ListPlot[Table[Cos[x*2*Pi/32] * HannWindow[(x-16)/32], {x, 0, 31}]] ListLinePlot[Table[Cos[x*2*Pi/64], {x, 0, 63}]] ListLinePlot[Abs[Fourier[Table[Cos[x*2*Pi/60], {x, 0, 63}]]], PlotRange->{0, 8}] ListLinePlot[Abs[Fourier[Table[Cos[x*2*Pi/60] * HannWindow[(x-32)/64], {x, 0, 63}]]], PlotRange->{0, 8}]

Fourier Cosine Transform

Fourier Cosine Transform / Fourier Sine Transform

複數波太前衛,實數波仍有人用。cos波和sin波都有人用,cos波比較常見。

傅立葉餘弦轉換是一對一函數,輸入和輸出都是一串實數,可以是連續函數或者離散數列,各有對應名稱。

輸入 輸出 名稱
連續 連續 Fourier Cosine Transform
連續 離散 Fourier Cosine Series
離散 連續 似乎沒有名稱
離散 離散 Discrete Cosine Transform

遇到不連續函數容易失真,例如Gibbs Phenomenon。讀者可自行尋找資料。

Lapped Transform

https://en.wikipedia.org/wiki/Modified_discrete_cosine_transform

Fourier Bessel Transform(Under Construction!)

複數波改成Bessel函數。

Zernike Polynomials

Fourier Ramanujan Transform

http://perso.ens-lyon.fr/patrick.flandrin/NewOrl2.pdf

複數波改成Ramanujan Sum。

norm從N改成φ(n)。垂直是互質。逆轉換需要定義無限長數列的內積。

Laplace Transform

拉普拉斯轉換是傅立葉轉換的進化版本。有兩個地方不同:

一、e-i⋅(2π/N)⋅n⋅k的次方值,改成任意複數。

次方值的實部,影響振幅;次方值的虛部,影響相位、頻率。

傅立葉轉換是振幅為1、相位為0、頻率為定值,平穩振動的波;拉普拉斯轉換是振幅頻率相位會變動的波(其中一例是Bessel函數),窮舉所有變動方式。

二、積分起點-oo,改成0。

傅立葉轉換處理負索引值;拉普拉斯轉換不處理負索引值,符合真實世界常見情況。

Wavelet

Wavelet Transform

http://en.wikipedia.org/wiki/Wavelet

Linear Canonical Transform

http://en.wikipedia.org/wiki/Linear_canonical_transformation

Gyrator Transform
http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=7112636