演算法筆記 - Wave

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

Fourier Cosine Transform

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

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

離散到離散的餘弦轉換,輸入和輸出都是一串數列。

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

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

連續版本是離散版本的推廣:輸入輸出無限密無限長。

Discrete Cosine Transform物理意義

N個波,頻率是0倍、0.5倍、1倍、1.5倍、……,分別是cos((2π/N)⋅0⋅t)、cos((2π/N)⋅0.5⋅t)、cos((2π/N)⋅1⋅t)、……。寫成代數是cos((2π/N)⋅(f/2)⋅t)。

輸入數列與一個波,置中對齊。N個對應位置,相乘後求和(點積),得到一個輸出數值。

輸入數列,分別投影至N個波,得到N個輸出數值,形成輸出數列。這就是餘弦轉換。

正向餘弦轉換:一個複雜的波,拆解成N個平穩的波,頻率是0倍開始漸增0.5倍,振幅是N個輸出數值,相位都是0。

逆向餘弦轉換:N個平穩的波,頻率是0倍開始漸增0.5倍,分別乘上振幅,疊加成一個複雜的波。

Discrete Cosine Transform有許多版本

N個平穩的波,微調頻率振幅相位。維基百科列出了許多版本,大家習慣以第2型做為正向轉換、以第3型做為逆向轉換。

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

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

餘弦轉換必須指定相位。傅立葉轉換可以自動得到相位。小波轉換可以自由設計波形,不必是平穩的波。

餘弦轉換與其他轉換相比顯得大費周折。不過有一種說法是:大量資料,分段處理,而微調有助於銜接段落。圖片壓縮JPEG、影片壓縮H.265和VP9、聲音壓縮MP3和AAC都使用餘弦轉換。

Discrete Cosine Transform數學公式

正向餘弦轉換
       N-1
y[f] =  ∑ { x[t] ⋅ cos((2π/N)⋅(f/2)⋅(t+0.5)) } ⋅ √2/N
       t=0
       y[0]最後再除以√2

逆向餘弦轉換(反函數)
       N-1
x[t] =  ∑ { y[f] ⋅ cos((2π/N)⋅(f/2)⋅(t+0.5)) } ⋅ √2/N
       f=0
       y[0]事先要除以√2

符號意義:輸入數列x、輸出數列y、數列長度N、時刻t、頻率倍數f/2。時刻加上0.5以置中對齊。

Discrete Cosine Transform是線性函數

餘弦轉換是線性函數!可以寫成矩陣形式!

正向餘弦轉換,視角是橫條點積投影,可以畫成這樣子:

逆向餘弦轉換,視角是直條加權總和,可以畫成這樣子:

演算法(公式解)

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

演算法(Arai-Agui-Nakajima Algorithm)

針對輸入輸出只有8點的餘弦轉換,進行細部加速。

Integer Discrete Cosine Transform

實數運算既複雜又緩慢。改弦易轍,以整數運算趨近正確答案:

https://www.vcodex.com/h264avc-4x4-transform-and-quantization/

整數運算比實數運算簡單,整數餘弦轉換比餘弦轉換迅速。

影像處理,訊號都是整數,習慣採用整數餘弦轉換。

2D Discrete Cosine Transform

餘弦轉換可以推廣到高維度。

二維餘弦轉換,輸入和輸出都是一個N×N方陣。輸入方陣,分別除以N×N種波,得到N×N個輸出數值,形成輸出方陣。

Plot3D[Cos[1.5x]Cos[1.5y], {x, 0, 2Pi}, {y, 0, 2Pi}, PlotRange -> {-1, 1}, Axes -> False, ColorFunction -> (ColorData["CherryTones"][Rescale[#3, {-2, 2}]] &)]

依照公式實作,時間複雜度為O(N⁴)。快速演算法是每一橫條各自餘弦轉換,然後每一直條各自餘弦轉換,時間複雜度為O(NN² + NN²) = O(N³)。

Gibbs Phenomenon

連續到離散的餘弦轉換,有個缺點:先拆解再疊加,產生針刺。函數曲線劇烈起伏之處(斜率很大或者很小)尤其明顯。

彷彿多項式內插的Runge Phenomenon。

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,俯瞰和側視,即是實數波cosθ、sinθ。

換句話說,觀察e = cosθ + i⋅sinθ這道式子:取實部得到實數波cosθ、取虛部得到實數波sinθ。

波有兩種繪圖方式:三維空間螺旋線、複數平面繞圓圈。

Fourier Transform

Fourier Transform

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

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

離散到離散的傅立葉轉換,輸入和輸出都是一串複數數列。

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

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

連續版本是離散版本的推廣:輸入輸出無限密無限長。

Discrete Fourier Transform物理意義

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

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

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

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

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

Discrete Fourier Transform數學公式

正向傅立葉轉換
       N-1
y[f] =  ∑ { x[t] ÷ ei⋅(2π/N)⋅f⋅t } ÷ √N
       t=0
       N-1
     =  ∑ { x[t] ⋅ e-i⋅(2π/N)⋅f⋅t } ÷ √N
       t=0

逆向傅立葉轉換(反函數)
       N-1
x[t] =  ∑ { y[f] ⋅ ei⋅(2π/N)⋅f⋅t } ÷ √N
       f=0

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

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

Discrete Fourier Transform是線性函數

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

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

演算法(Horner's Rule)

依照公式實作,點積視作多項式,時間複雜度是O(N²)。

演算法(Cooley-Tukey Algorithm)

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

這裡介紹最經典的快速傅立葉轉換。公式的偶數項與奇數項分開整理,採用Dynamic Programming,時間複雜度是O(NlogN)。

由於必須剛好對半分,所以N必須剛好是2的次方。當N不是2的次方,可在輸入數列末端補零,理由容後介紹。

逆向轉換的演算法也是一樣的,此處省略。

                           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₇)轉換結果第1項
   = 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)。

Hartley Transform

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

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

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

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

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

= y[f] ⋅ cos((2π/N)⋅f⋅1) + y[-f] ⋅ sin((2π/N)⋅f⋅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θ

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

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

實數運算比複數運算簡單,哈特利轉換比傅立葉轉換迅速。

聲音處理,訊號都是實數,習慣採用哈特利轉換,再把結果調整成傅立葉轉換。

2D Discrete Fourier Transform

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

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

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

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 Transform的性質

輸入輸出對應

連續到連續的傅立葉轉換,輸入輸出有著特殊對應關係。

因為正向轉換幾乎等於逆向轉換,所以輸入輸出對調之後,對應依然成立。

運算對應

加法-加法         a + b = aft + bft
倍率-倍率         a ⋅ k = aft ⋅ k
乘法-摺積         a × b = aft ∗ bft
摺積-乘法         a ∗ b = aft × bft
微分-角加速度     a′ =  aft ⋅ 2πif     d/dt a(t) = 2πif ⋅ aft(f) 
角加速度-微分     a ⋅ 2πit = aft′
平方和(能量)守恆  ‖a‖ = ‖aft‖          ∑ a(t)² = ∑ aft(f)² 

因為正向轉換幾乎等於逆向轉換,所以運算對調之後,對應依然成立。

Laplace Transform

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

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

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

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

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

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

計算學家不使用拉普拉斯轉換。但是因為上述性質通通可以推廣到拉普拉斯轉換,所以訊號處理教科書喜歡採用拉普拉斯轉換。

Wavelet(Under Construction!)

Wavelet

「小波」,自訂特殊造型的波。

http://www.faculty.jacobs-university.de/llinsen/teaching/320491/Lecture05.pdf

Wavelet Transform

「小波轉換」是一對一函數,輸入和輸出都是一串實數,可以是離散數列或者連續函數。

哪些小波可以用於小波轉換呢?以線性代數的觀點來看,N個波必須構成N維空間,才有反矩陣。換句話說,線性獨立導致一對一函數。

餘弦轉換、傅立葉轉換,N個波線性獨立,是一對一函數。

哪些小波可以用於小波轉換呢?以線性代數的觀點來看,最簡潔的方式是正規正交基底,反矩陣就是轉置矩陣。正規是指個個範數為1。正交是指兩兩內積為0。換句話說,N個小波能量皆是1(一個小波的N個位置的平方的總和是1),N個小波互相垂直。

傅立葉轉換是正規正交基底,N個波互相垂直,再除以√N使得能量皆是1。

演算法

時間複雜度優於O(N²)的小波轉換演算法,老人家稱作「快速小波轉換Fast Wavelet Transform, FWT」。

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