Wave (ℝ)

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

振動、振盪

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

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

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

振動可以用函數表示

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

平穩的振動

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

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

振動的快慢:頻率

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

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

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

振動的高低:振幅

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

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

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

振動的起點:相位

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

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

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

振動有疊加效果

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

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

振動有傳遞效果:波

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

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

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

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

Fourier Cosine Transform(Under Construction!)

Fourier Cosine Transform

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

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

連續到連續的傅立葉餘弦轉換,輸入和輸出都是一個ℝ ⇨ ℝ函數。

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

電腦做運算,數值皆離散。計算學家只使用離散版本。

稍後提到的傅立葉轉換,性質優美,完全可以取代傅立葉餘弦轉換。現在只有圖片壓縮、聲音壓縮使用傅立葉餘弦轉換,原因好像是區塊之間的連續性。

另外也有sin版本,不過沒人用。

Fourier Cosine Transform物理意義

N個cos波,頻率是0倍到N-1倍,分別是cos((2π/N)⋅0⋅θ)、cos((2π/N)⋅1⋅θ)、……、cos((2π/N)⋅(N-1)⋅θ)。寫成代數是cos((2π/N)⋅n⋅θ)。

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

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

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

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

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

Gibbs Phenomenon

Lapped Transform

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

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²)。

演算法(Cooley-Tukey Algorithm)

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

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

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

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

                           FFT
(x₀ x₁ x₂ x₃ x₄ x₅ x₆ x₇) ----> (y₀ y₁ y₂ y₃ y₄ y₅ y₆ y₇)

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

y₀ = x₀ω⁰ + x₁ω⁰ + x₂ω⁰ + x₃ω⁰ + x₄ω⁰ + x₅ω⁰ + x₆ω⁰ + x₇ω⁰
   = (x₀ω⁰ + x₂ω⁰ + x₄ω⁰ + x₆ω⁰) + (x₁ω⁰ + x₃ω⁰ + x₅ω⁰ + x₇ω⁰)
   = (x₀ω⁰ + x₂ω⁰ + x₄ω⁰ + x₆ω⁰) + ω⁰ ⋅ (x₁ω⁰ + x₃ω⁰ + x₅ω⁰ + x₇ω⁰)
   = (x₀ x₂ x₄ x₆)轉換結果第0項 + ω⁰ ⋅ (x₁ x₃ x₅ x₇)轉換結果第0項
   = y偶0 + ω⁰ ⋅ y奇0

y₁ = x₀ω⁰ + x₁ω¹ + x₂ω² + x₃ω³ + x₄ω⁴ + x₅ω⁵ + x₆ω⁶ + x₇ω⁷
   = (x₀ω⁰ + x₂ω² + x₄ω⁴ + x₆ω⁶) + (x₁ω¹ + x₃ω³ + x₅ω⁵ + x₇ω⁷)
   = (x₀ω⁰ + x₂ω² + x₄ω⁴ + x₆ω⁶) + ω¹ ⋅ (x₁ω⁰ + x₃ω² + x₅ω⁴ + x₇ω⁶)
   = (x₀υ⁰ + x₂υ¹ + x₄υ² + x₆υ³) + ω¹ ⋅ (x₁υ⁰ + x₃υ¹ + x₅υ² + x₇υ³)
   = (x₀ x₂ x₄ x₆)轉換結果第1項 + ω¹ ⋅ (x₁ x₃ x₅ x₇)轉換結果第2項
   = y偶1 + ω¹ ⋅ y奇1

y₂ = x₀ω⁰ + x₁ω² + x₂ω⁴ + x₃ω⁶ + x₄ω⁸ + x₅ω¹⁰ + x₆ω¹² + x₇ω¹⁴
   = (x₀ω⁰ + x₂ω⁴ + x₄ω⁸ + x₆ω¹²) + (x₁ω² + x₃ω⁶ + x₅ω¹⁰ + x₇ω¹⁴)
   = (x₀ω⁰ + x₂ω⁴ + x₄ω⁸ + x₆ω¹²) + ω² ⋅ (x₁ω⁰ + x₃ω⁴ + x₅ω⁸ + x₇ω¹²)
   = (x₀υ⁰ + x₂υ² + x₄υ⁴ + x₆υ⁶ ) + ω² ⋅ (x₁υ⁰ + x₃υ² + x₅υ⁴ + x₇υ⁶ )
   = (x₀ x₂ x₄ x₆)轉換結果第2項 + ω² ⋅ (x₁ x₃ x₅ x₇)轉換結果第2項
   = y偶2 + ω² ⋅ y奇2

y₃ = x₀ω⁰ + x₁ω³ + x₂ω⁶ + x₃ω⁹ + x₄ω¹² + x₅ω¹⁵ + x₆ω¹⁸ + x₇ω²¹
   = (x₀ω⁰ + x₂ω⁶ + x₄ω¹² + x₆ω¹⁸) + (x₁ω³ + x₃ω⁹ + x₅ω¹⁵ + x₇ω²¹)
   = (x₀ω⁰ + x₂ω⁶ + x₄ω¹² + x₆ω¹⁸) + ω³ ⋅ (x₁ω⁰ + x₃ω⁶ + x₅ω¹² + x₇ω¹⁸)
   = (x₀υ⁰ + x₂υ³ + x₄υ⁶  + x₆υ⁹ ) + ω³ ⋅ (x₁υ⁰ + x₃υ³  + x₅υ⁶ + x₇υ⁹ )
   = (x₀ x₂ x₄ x₆)轉換結果第3項 + ω³ ⋅ (x₁ x₃ x₅ x₇)轉換結果第3項
   = y偶3 + ω³ ⋅ y奇3

注意到 ω⁸ = 1
y₄ = x₀ω⁰ + x₁ω⁴ + x₂ω⁸ + x₃ω¹² + x₄ω¹⁶ + x₅ω²⁰ + x₆ω²⁴ + x₇ω²⁸
   = (x₀ω⁰ + x₂ω⁸ + x₄ω¹⁶ + x₆ω²⁴) + (x₁ω⁴ + x₃ω¹² + x₅ω²⁰ + x₇ω²⁸)
   = (x₀ω⁰ + x₂ω⁸ + x₄ω¹⁶ + x₆ω²⁴) + ω⁴ ⋅ (x₁ω⁰ + x₃ω⁸ + x₅ω¹⁶ + x₇ω²⁴)
   = (x₀ω⁰ + x₂ω⁰ + x₄ω⁰  + x₆ω⁰ ) + ω⁴ ⋅ (x₁ω⁰ + x₃ω⁰ + x₅ω⁰  + x₇ω⁰ )
   = (x₀ x₂ x₄ x₆)轉換結果第0項 + ω⁴ ⋅ (x₁ x₃ x₅ x₇)轉換結果第0項
   = y偶0 + ω⁴ ⋅ y奇0

y₅ y₆ y₇以此類推
y₀ = y偶0 + y奇0 ⋅ ω⁰
y₁ = y偶1 + y奇1 ⋅ ω¹
y₂ = y偶2 + y奇2 ⋅ ω²
y₃ = y偶3 + y奇3 ⋅ ω³
y₄ = y偶0 + y奇0 ⋅ ω⁴
y₅ = y偶1 + y奇1 ⋅ ω⁵
y₆ = y偶2 + y奇2 ⋅ ω⁶
y₇ = y偶3 + y奇3 ⋅ ω⁷

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

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

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

2D Fourier Transform

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

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

依照公式實作,時間複雜度為O(N⁴)。快速的演算法,是每一橫條各自傅立葉轉換,然後每一直條各自傅立葉轉換,時間複雜度為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}]

Laplace Transform

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

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

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

傅立葉轉換是振幅為1、相位為0、頻率為定值,平穩振動的波;拉普拉斯轉換是振幅頻率相位隨時變動的波,窮舉所有變動方式。

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

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

據我所知,計算學家不使用拉普拉斯轉換。

Wavelet

Wavelet Transform

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

Linear Canonical Transform

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

Gyrator Transform
https://arxiv.org/abs/1707.03689