wave [ℝ]

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

振動、振盪

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

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

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

振動可以用函數表示

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

平穩的振動

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

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

振動的快慢:頻率

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

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

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

振動的高低:振幅

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

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

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

振動的起點:相位

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

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

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

振動有疊加效果

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

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

振動有傳遞效果:波

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

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

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

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

振動通常有規律:週期函數

水波、聲波、地震波、電磁波,擷取一段極短時間,放大尺度仔細觀察,可以發現同一種振動模式連接不斷重複出現,振動模式大致固定不變。畫成函數圖形,形成週期函數。

週期函數習慣稱作波。此波非彼波。

「振動vibration」是物理術語。一個粒子來回運動。

「波wave」是物理術語。一群粒子傳遞振動、各自振動。

振動與波的頻率振幅相位通通不同,千萬不要混為一談。

振動與波可能有規律、可能沒有規律,千萬不要混為一談。

然而,波通常有規律,形成週期函數。而且波隨處可見。以致於大家將經典的週期函數通通稱作波。這裡列出幾個例子:

mathematics                  | engineering
-----------------------------| --------------------------
periodic rectanglar function | square wave
periodic triangular function | triangle wave
periodic linear function     | sawtooth wave
sinusoidal function          | sinusoidal wave / sinusoid
sine function                | sine wave
cosine function              | cosine wave
數學     丨工程
一一一一一一一十一一一
週期長方形函數丨方波
週期三角形函數丨三角波
週期一次函數 丨鋸齒波
弦型函數   丨弦波
正弦函數   丨正弦波
餘弦函數   丨餘弦波

弦波

平穩的振動,形成弦波。正弦波、餘弦波都是弦波。

這裡深入介紹弦波的數學式子與函數圖形,教科書經常用到。

首先以拋物線y = x²為例。這不是弦波,只是前情提要。

調整函數圖形,數學式子寫成反運算。

例如x值增加1、y值增加2,數學式子寫成減法,形成(y-2) = (x-1)²。移項整理成函數的模樣,形成y = (x-1)² + 2。

例如x值乘上3倍、y值乘上4倍,數學式子寫成除法,形成(y/4) = (x/3)²。移項整理成函數的模樣,形成y = 4(x/3)²。

接著以餘弦函數y = cos(x)為例。這是弦波,最常用到。

x是角度,y是長度。一律改成長度。x從角度改成長度(圓周長度改成半徑長度),x值除以2π,形成y = cos(2πx)。

或者這樣說,x重複範圍[0,2π],y範圍[-1,+1]。一律改成範圍1。x值除以2π,範圍調整成[0,1],形成y = cos(2πx)。

考慮振幅與週期。x值乘上週期T倍、y值乘上振幅A倍,形成(y/A) = cos(2πx/T)。移項變成y = A cos(2πx/T)。

除法不好算。那就改成乘法。物理學家規定週期的倒數是頻率1/T = f,形成y = A cos(2πfx)。這就是頻率的由來。

對於振動而言,輸入應是現在時刻,輸出應是振動高低。把長度x改成時間t,把長度y改成函數f(t),形成f(t) = A cos(2πft)。

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型作為逆向轉換。

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

餘弦轉換與其他轉換相比顯得大費周折。不過有一種說法是:大量資料,分段處理,而微調有助於銜接段落。圖片壓縮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

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

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

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

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

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

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

複數亦可以畫成圖形。

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

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

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

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

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

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

Euler's formula

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

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

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

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

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

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

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

complex sinusoid

因為e𝑖θ長得像波,所以用e𝑖θ將實弦波推廣成複弦波。

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

複弦波e𝑖θ,俯瞰和側視,即是實弦波cosθ、sinθ。

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

最後引入振幅A、頻率f,大功告成。

實弦波real sinusoid:f(t) = A cos(2πft)。

複弦波complex sinusoid:f(t) = A e𝑖2πft

Fourier transform

Fourier transform

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

輸入丨輸出丨名稱
一一十一一十一一一一一一一一一一一一一一一一一一一
離散丨離散丨discrete Fourier transform
離散丨連續丨discrete-time Fourier transform
連續丨離散丨Fourier series
連續丨連續丨Fourier transform

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

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

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

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

discrete Fourier transform物理意義

N個複弦波,頻率是0倍到N-1倍,分別是e𝑖⋅(2π/N)⋅0⋅t、e𝑖⋅(2π/N)⋅1⋅t、……、e𝑖⋅(2π/N)⋅(N-1)⋅t。寫成代數是e𝑖⋅(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] ÷ e𝑖⋅(2π/N)⋅f⋅t } ÷ √N
       t=0
       N-1
     =  ∑ { x[t] ⋅ e-𝑖⋅(2π/N)⋅f⋅t } ÷ √N
       t=0

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

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

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

discrete Fourier transform是線性函數

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

ω = e𝑖⋅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

物理教科書規定角速度ω = 2πf。數學教科書規定第一個N次單位根ω = e𝑖⋅2π/N。符號恰巧一樣。此處是後者。

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

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

演算法(公式解)

依照公式實作,時間複雜度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-𝑖⋅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(N)。假設高低位數顛倒是O(1)。

Hartley transform

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

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

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

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

  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):

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

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

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」,呼應傅立葉轉換的功能:把波(時間軸)表示成頻譜(頻率軸)。

時域乘法等於頻域循環卷積。請見本站文件「filter」。數列與窗函數相乘,等於數列與窗函數在頻域的循環卷積。

窗函數,多由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

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

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-𝑖⋅(2π/N)⋅f⋅t的次方值,改成任意複數。

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

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

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

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

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

wavelet

wavelet

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

教科書習慣介紹Harr wavelet,一個方波。

其他小波:

wavelet transform

wavelet transform

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

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

順帶一提,餘弦轉換、傅立葉轉換,N個弦波線性獨立,是雙射函數。

哪些小波可以用於小波轉換呢?以線性代數的觀點來看,最簡潔的方式是正規正交基底,反矩陣就是轉置矩陣。正規是指個個範數為1(向量長度是1)。正交是指兩兩內積為0(向量互相垂直)。

順帶一提,傅立葉轉換是正規正交基底,N個弦波互相垂直,再除以√N使得長度(能量)皆是1。

演算法

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

應用

上世紀末曾經流行一陣子。現在乏人問津。

《Application of Wavelet Transform and its Advantages Compared to Fourier Transform》