Wave
程度★ 難度★★★
振動、振盪
這個世界天天都在振動。地面、空氣、海水、機械、人體等等，都是不斷振動。震動、震盪是自古以來就有的詞彙。
振動可以用函數表示
每個時間點的振動高低，可以描繪成函數圖形，橫向是時間軸，縱向是每個時刻的振動高低位置。
平穩的振動
最平穩的振動，就是高中物理教的簡諧運動，是等速圓周運動在座標軸上的投影，呈現 sin 函數。 sin 函數和 cos
函數長相一樣，只是起點不同而已。
舉例來說，敲打音叉產生的振動，就非常接近平穩的振動。
振動的快慢：頻率
單位時間振動的次數，稱做「頻率（ Frequency ）」。一秒振動的次數，單位是赫茲 Hz 。
人類能感知頻率：耳朵能聽到 20Hz 至 20000Hz 的空氣振動，低頻低沉、高頻尖銳；眼睛能看到 4×10^14Hz 至 8×10^14Hz
的電磁振盪，低頻至高頻分別呈現紅橙黃綠藍靛紫。
振動的高低：振幅
振動的最高（低）距離，稱做「振幅（ Amplitude ）」。
人類能感知振幅，受頻率大小影響。就聽覺而言，振幅高則大聲、振幅低則小聲；就視覺而言，振幅高則亮、振幅低則暗。題外話，人類對於頻率與振幅的區分能力，大略等於取 log 。
振動的起點：相位
振動的起點，稱做「相位（ Phase ）」。注意到，相位是圓周運動 cos 函數的角度，而不是振動高低位置。物理學家喜歡用相位。
人類幾乎分辨不出相位的差異。
振動有疊加效果
現實世界當中，多個振動時常融合成一個振動，等於各個振動高低位置相加。相同方向則增益、相反方向則抵銷。
寫成數學式子，就是多個函數相加。
振動有傳遞效果：波
一個粒子振動，就會導致隔壁粒子振動，一傳十、十傳百。宏觀之下，形成「波（ Wave ）」。傳遞速度取決於粒子之間的作用力、粒子的質量。作用力強、質量小，則傳遞速度快。
均勻分佈的粒子之中，某個粒子振動所產生的波，剛好也呈現 sin 函數，英文稱做 Sine Wave 或者 Sinusoid 。
http://140.111.56.210/file/ph/phy_6-1-1-1_2/1.html
水的高低起伏，就是水波。空氣的疏密，就是聲波。地的高低起伏與左右晃動，就是地震波。電場與磁場的交互作用，就是電磁波。光波經實驗證明是電磁波。原子的振動，也許是熱。有人覺得氣功也許是波，就叫做氣功波。
觀察波上任意一點，都是在振動。許多情況下，民眾喜歡用波來指稱振動。
Complex Number
程度★ 難度★★★
有物混成，先天地生，寂兮寥兮，
獨立而不改，周行而不殆，可以為天下母。《老子》
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 iθ = cosθ + i * sinθ ，複數的長度是常數 1 ，複數的角度是變數 θ 。等式右邊，是將長度與角度，還原成一個複數，外觀很複雜但是本質很簡單。
有了歐拉公式，一個複數也可以重新表示成 e 的次方、另乘上倍率。次方值即是角度乘 i ，倍率即是長度。
歐拉公式，定量增加 θ ，在複數平面上，外觀宛如「等速圓周運動」，逆時針繞圈；只看實部或者只看虛部，外觀宛如「簡諧運動」，先上後下。
繞 360° 是一圈，剛好回到 +1 ；繞 180° 是半圈，剛好是 -1 。因此有了 e iπ + 1 = 0 這條著名等式， π
就是 180° 。這個 e ，大約是 2.71828183 ，是自然對數的底數 e ，是 1/x 積分後所出現的 e 。離題了。
Wave
e iθ 運算簡單，考慮長度與角度即可。 e iθ 性質優美，每轉 90° 剛好是 ±1 與 ±i
。也許你會漸漸愛上它。現在有了神奇數字 e ，就改用 e iθ 來表示波。
要把 e iθ 還原成 sin 波、 cos 波，有兩種方式。觀察 e iθ = cosθ + i * sinθ 這道式子：
一種是取實部得到 cos 波、取虛部得到 sin 波。
一種是用 e iθ 與 e -iθ ，相加除以二得到 cos 波，相減除以二得到 sin 波。
Fourier Transform
程度★ 難度★★★
Fourier Transform
傅立葉轉換的輸出入都是一串複數，是一對一函數。輸入可以是連續函數或者離散數列，輸出亦如是。根據連續與離散的差異，傅立葉轉換有著不同的名稱。混淆視聽罷了。
輸入 輸出 名稱 連續 連續 Fourier Transform 連續 離散 Discrete Fourier Transform 離散 連續 Continuous Fourier Transform 離散 離散 Discrete-time Fourier Transform電腦做運算，數值皆離散，最常用到離散到離散的傅立葉轉換。公式如下：
傅立葉轉換：
N-1
y[n] = Σ { x[k] ÷ ei*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
傅立葉轉換是線性變換，其矩陣恰是正交基底。W = e^(i*2π*/N) [ y0 ] [ W-0*0 W-0*1 W-0*2 ... W-0*(N-1) ] [ x0 ] [ y1 ] [ W-1*0 W-1*1 W-1*2 ... W-1*(N-1) ] [ x1 ] [ . ] = [ . . . . ] * [ . ] [ . ] [ . . . . ] [ . ] [ yN-1 ] [ W-(N-1)*0 W-(N-1)*1 W-(N-1)*2 ... W-(N-1)*(N-1) ] [ xN-1 ] [ x0 ] [ W0*0 W0*1 W0*2 ... W0*(N-1) ] [ y0 ] [ x1 ] 1 [ W1*0 W1*1 W1*2 ... W1*(N-1) ] [ y1 ] [ . ] = --- [ . . . . ] * [ . ] [ . ] N [ . . . . ] [ . ] [ xN-1 ] [ W(N-1)*0 W(N-1)*1 W(N-1)*2 ... W(N-1)*(N-1) ] [ yN-1 ]傅立葉轉換可以推廣到高維度。舉例來說，二維的傅立葉轉換，輸出入都是一個複數矩陣，轉換過程是：橫向每一條先各自傅立葉轉換，然後直向每一條再各自傅立葉轉換。
Fourier Transform
只給數學公式，讀者應該是霧裡看花。N 個波，頻率是 n=0 倍到 n=N-1 倍，分別是 e i*2π*(n/N) 。
（學過線性代數的讀者，也可以想做內積、投影。）
輸入數列分別除以 N 個波，得到 N 個輸出數值。這就是傅立葉轉換。
正向傅立葉轉換，是把一個複雜的波，拆解成 N 個平穩的波，頻率是 0 倍到 N-1 倍，振幅與相位是 N 個輸出數值。逆向傅立葉轉換，是把 N 個平穩的波，頻率是 0 倍到 N-1 倍，分別乘上振幅、加上相位，再疊加成一個複雜的波，
Frequency Spectrum
傅立葉轉換輸出數列有 N 個複數，可以畫成函數──一般不畫實部與虛部，而是畫長度與角度，具備物理意義。這 N 個複數的長度（振幅）畫成函數，稱為「振幅頻譜」。
這 N 個複數的角度（相位）畫成函數，稱為「相位頻譜」。
兩者合稱為「頻譜」。
輸入數列是一串實數時，輸出數列會前後對稱並且共軛（長度相等、角度負號）。此時頻譜可以只畫一半。我們得以運用傅立葉轉換分解一個波，運用逆向傅立葉轉換合成一個波。運用「頻譜」，就能解讀波的詳細內容。
演算法
按照公式實作，時間複雜度是 O(N^2) 。- const double π = 2.0 * acos(0);
- const int N = 10;
- complex<double> x[N], y[N];
- // cos()和sin()並不是O(1)。
- // 內部迴圈不斷呼叫cos()與sin()，因此整體不是O(N^2)。
- void fourier_transform()
- {
- double ω = 2.0 * π / N;
- for (int i=0; i<N; i++)
- {
- y[i] = 0;
- for (int j=0; j<N; j++)
- y[i] += x[j] / complex<double>(cos(ω*i*j), sin(ω*i*j));
- }
- }
- const int N = 10;
- complex<double> x[N], y[N];
- void fourier_transform()
- {
- double ω = -2.0 * 2.0 * acos(0) / N;
- complex<double> eω(cos(ω), sin(ω)), edθ(1, 0);
- for (int i=0; i<N; i++, edθ *= eω)
- {
- complex<double> eθ(1, 0);
- y[i] = 0;
- for (int j=0; j<N; ++j, eθ *= edθ)
- y[i] += x[j] * eθ;
- }
- }
- const int N = 10;
- double xR[N], xI[N]; // 輸入，實部與虛部。
- double yR[N], yI[N]; // 輸出，實部與虛部。
- // 三角函數和角公式，求得 θ += dθ 之後的三角函數值。
- void add(double& cosθ, double& sinθ, double sindθ, double cosdθ)
- {
- double temp_cosθ = cosθ;
- cosθ = cosθ * cosdθ - sinθ * sindθ;
- sinθ = temp_cosθ * sindθ + sinθ * cosdθ;
- }
- void fourier_transform()
- {
- double ω = -2.0 * 2.0 * acos(0) / N;
- double cosω = cos(ω), sinω = sin(ω);
- double cosdθ = 1, sindθ = 0;
- for (int i=0; i<N; i++)
- {
- double cosθ = 1, sinθ = 0;
- yR[i] = yI[i] = 0;
- for (int j=0; j<N; ++j)
- {
- yR[j] += xR[j] * cosθ - xI[j] * sinθ;
- yI[j] += xR[j] * sinθ + xI[j] * cosθ;
- add(cosθ, sinθ, cosdθ, sindθ);
- }
- add(cosdθ, sindθ, cosω, sinω);
- }
- }
演算法（ Cooley-Tukey Algorithm ）
時間複雜度優於 O(N^2) 的傅立葉轉換演算法，老人家就直接稱做「快速傅立葉轉換（ Fast Fourier Transform, FFT
）」。這裡介紹最經典的快速傅立葉轉換，把公式的偶數項與奇數項分開整理，用 Divide and Conquer 解決，時間複雜度是 O(NlogN) 。由於必須剛好對半分，所以 N 必須剛好是 2 的次方。
y[n]
N-1
= Σ { x[k] ÷ ei*2π*(n/N)*k } 傅立葉轉換公式
k=0
N-1 N-1 拆成偶數項與奇數項
= Σ { x[k] ÷ ei*2π*(n/N)*k } + Σ { x[k] ÷ ei*2π*(n/N)*k }
k=0,2,4,... k=1,3,5,...
N-1 N-1 奇數項全部左移一項，然後提出常數
= Σ { x[k] ÷ ei*2π*(n/N)*k } + Σ { x«1[k] ÷ ei*2π*(n/N)*k } ÷ ei*2π*(n/N)*1
k=0,2,4,... k=0,2,4,...
N/2-1 N/2-1 索引值調成連續，點數減半了
= Σ { x偶[k] ÷ ei*2π*(n/N)*2k } + Σ { x奇[k] ÷ ei*2π*(n/N)*2k } ÷ ei*2π*(n/N)*1
k=0,1,2,... k=0,1,2,...
N/2-1 N/2-1 可視作頻率減半、波長變兩倍
= Σ { x偶[k] ÷ ei*2π*(n/(N/2))*k } + Σ { x奇[k] ÷ ei*2π*(n/(N/2))*k } ÷ ei*2π*(n/N)*1
k=0,1,2,... k=0,1,2,...
= y偶[n] + y奇[n] ÷ ei*2π*(n/N)*1 剛好是偶數項FFT與奇數項FFT
e.g. N = 4 y[0] = y偶[0] + y奇[0] ÷ ei*2π*(0/N) y[1] = y偶[1] + y奇[1] ÷ ei*2π*(1/N) y[2] = y偶[0] + y奇[0] ÷ ei*2π*(2/N) y[3] = y偶[1] + y奇[1] ÷ ei*2π*(3/N)
計算原理並不難解釋，反而是實作比較複雜。資料結構採用陣列，偶數點和奇數點並不是連續的元素──該如何整理成連續的元素，然後遞迴地解決問題呢？最常見的方式，是觀察分治法的分割順序與合併順序，然後一開始就重新排列每個元素至定位，配合分治法的順序。
- const double π = 2.0 * acos(0);
- const int N = 8;
- complex<double> x[N];
- void FFT()
- {
- // reverse bit and replace
- for (int i=1, j=0; i<N; ++i)
- {
- for (int k=N>>1; !((j^=k)&k); k>>=1) ;
- // for (int k=N>>1; k>(j^=k); k>>=1) ;
- if (i>j) swap(x[i], x[j]);
- // if (i<j) swap(x[i], x[j]);
- }
- for (int k=2; k<=N; k<<=1)
- {
- double ω = -2.0 * π / k;
- complex<double> dθ(cos(ω), sin(ω));
- // 每k個做一次FFT
- for (int j=0; j<N; j+=k)
- {
- // 前k/2個與後k/2的三角函數值恰好對稱，
- // 因此兩兩對稱的一起做。
- complex<double> θ(1, 0);
- for (int i=j; i<j+k/2; i++)
- {
- complex<double> a = x[i];
- complex<double> b = x[i + k/2] * θ;
- x[i] = a + b;
- x[i + k/2] = a - b;
- θ *= dθ;
- }
- }
- }
- }
- void FFT()
- {
- // reverse bit and replace
- ......
- for (int k=2; k<=N; k<<=1)
- {
- double ω = -2.0 * π / k;
- complex<double> dθ(cos(ω), sin(ω));
- // 對調內外迴圈，讓θ少乘幾次。
- // 缺點則是索引值更容易跳動，更容易產生cache miss。
- complex<double> θ(1, 0);
- for (int i=0; i<k/2; i++)
- {
- for (int j=i; j<N; j+=k)
- {
- complex<double> a = x[j];
- complex<double> b = x[j + k/2] * θ;
- x[j] = a + b;
- x[j + k/2] = a - b;
- }
- θ *= dθ;
- }
- }
- }
- void FFT()
- {
- // 整個反過來算
- for (int k=N; k>=2; k>>=1)
- {
- double ω = -2.0 * π / k;
- complex<double> dθ(cos(ω), sin(ω));
- complex<double> θ(1, 0);
- for (int i=0; i<k/2; i++)
- {
- for (int j=i; j<N; j+=k)
- {
- // θ變成相減後才乘上。
- complex<double> a = x[j];
- complex<double> b = x[j + k/2];
- x[j] = a + b;
- x[j + k/2] = (a - b) * θ;
- }
- θ *= dθ;
- }
- }
- // reverse bit and replace
- ......
- }
Hartley Transform
程度★★ 難度★★
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
= Σ { x[k] * e-i*2π*(n/N)*k } ÷ sqrt(N)
k=0
N-1
= Σ { x[k] * ( cos(2π*(n/N)*k) - i * sin(2π*(n/N)*k) ) } ÷ sqrt(N)
k=0
哈特利轉換：
N-1
y[n] = Σ { x[k] * ( cos(2π*(n/N)*k) + sin(2π*(n/N)*k) ) } ÷ sqrt(N)
k=0
N-1
= Σ { x[k] * cas(2π*(n/N)*k) } ÷ sqrt(N)
k=0
哈特利轉換的輸出，可以調整成傅立葉轉換的輸出， O(N) ：http://mathworld.wolfram.com/HartleyTransform.html
實數運算比複數運算還要簡單，所以哈特利轉換比傅立葉轉換還要快速。聲音處理、影像處理時，訊號都是實數、甚至是整數，剛好也符合哈特利轉換的輸入格式。因此一般都是套用哈特利轉換進行計算，再把結果調整成傅立葉轉換。
演算法（ Bracewell's Algorithm ）
由於哈特利轉換與傅立葉轉換的公式幾乎相同，所以兩者的演算法也是一一對應。這裡介紹的也是運用 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/(N/2))*k) * cos(2π*(n/N)*1)
k=0,1,2,... + cas(-2π*(n/(N/2))*k) * sin(2π*(n/N)*1) ) }
N/2-1
= Σ { x[2k+1] * cas( 2π*(n/(N/2))*k) } * cos(2π*(n/N)*1)
k=0,1,2,...
N/2-1
+ Σ { x[2k+1] * cas(-2π*(n/(N/2))*k) } * sin(2π*(n/N)*1)
k=0,1,2,...
= y奇[n] * cos(2π*(n/N)*1) + y奇[-n] * sin(2π*(n/N)*1)
e.g. N = 8, θ = 2π / N y[0] = y偶[0] + y奇[0] * cos0θ + y奇[0] * sin0θ y[1] = y偶[1] + y奇[1] * cos1θ + y奇[3] * sin1θ y[2] = y偶[2] + y奇[2] * cos2θ + y奇[2] * sin2θ y[3] = y偶[3] + y奇[3] * cos3θ + y奇[1] * sin3θ y[4] = y偶[0] + y奇[0] * cos4θ + y奇[0] * sin4θ y[5] = y偶[1] + y奇[1] * cos5θ + y奇[3] * sin5θ y[6] = y偶[2] + y奇[2] * cos6θ + y奇[2] * sin6θ y[7] = y偶[3] + y奇[3] * cos7θ + y奇[1] * sin7θ下面是筆者所能找到的效率最好的實作程式碼：
http://home.iae.nl/users/mhx/fft.c
http://reocities.com/ResearchTriangle/8869/fft_summary.html
Number Theoretic Transform
程度★★ 難度★★
Number Theoretic Transform
數論轉換的輸出入都是一串餘數，是一對一函數。e 的純虛數次方會不斷繞圈，餘數的次方也會不斷繞圈。於是有人把傅立葉轉換的基底，從 e i*2π*(n/N) 改成原根的次方。
http://www.apfloat.org/prim.html
輸出入每個數值都是餘數，皆大於等於零、小於模數。當輸入數值很大，那麼模數也必須足夠大。
- // 預先算好的，支援到n = 1<<18。
- const long long m = 50000000001507329LL; // 質數 k*n+1
- const long long primitive_root = 3;
- // 快速數論轉換
- void FNT(long long* z, long long n, long long ω)
- {
- for (int k=n; k>=2; k>>=1)
- {
- long long θ = 1;
- for (int i=0; i<k/2; i++)
- {
- for (int j=i; j<n; j+=k)
- {
- long long a = z[j];
- long long b = z[j + k/2];
- z[j] = (a + b) % m;
- z[j + k/2] = mul((a - b + m) % m, θ, m);
- }
- θ = mul(θ, ω, m);
- }
- ω = mul(ω, ω, m);
- }
- for (int i=1, j=0; i<n; ++i)
- {
- for (int k=n>>1; k<(j^=k); k>>=1) ;
- if (i > j) swap(z[i], z[j]);
- }
- }
- void demo()
- {
- long long z[1024];
- long long n = 1024;
- long long r = pow(primitive_root, m/n, m);
- long long invr = pow(r, m-2, m);
- long long invn = pow(n, m-2, m);
- // 正向轉換
- FNT(z,n,r);
- // 逆向轉換
- FNT(z,n,invr);
- for (int i=0; i<n; ++i)
- z[i] = mul(z[i], invn, m);
- }
Convolution
程度★ 難度★★★
Multiplication
兩串數列做乘法運算，得到一個數列：對應項相乘。 a0 a1 aN-1
(a0,a1,...,aN-1) × (b0,b1,...,bN-1) = ( × , × ,..., × )
b0 b1 bN-1
a-1 a0 a1
(...,a-1,a0,a1,...) × (...,b-1,b0,b1,...) = (..., × , × , × ,...)
b-1 b0 b1
Dot Product
兩串數列做內積運算，得到一個值：對應項相乘後再求和。橫著寫叫數列，直著寫叫向量。數列內積就是向量內積。
a0 a1 aN-1
(a0,a1,...,aN-1) dot (b0,b1,...,bN-1) = × + × + ... + ×
b0 b1 bN-1
a-1 a0 a1
(...,a-1,a0,a1,...) dot (...,b-1,b0,b1,...) = ... + × + × + × + ...
b-1 b0 b1
內積是線性變換！有限長數列的內積，對應的矩陣，是橫的一條；無限長數列的內積，對應的矩陣，是無限大的矩陣。 [a0 ]
[a1 ] a0 a1 aN-1
[b0 b1 b2 ⋯ bN-1] * [a2 ] = × + × + ... + ×
[⋮ ] b0 b1 bN-1
[aN-1]
[⋮ ]
[a-1] a-1 a0 a1
[⋯ b-1 b0 b1 ⋯] * [a0 ] = ... + × + × + × + ...
[a1 ] b-1 b0 b1
[⋮ ]
Cross Product
外積，本篇用不到，略過。
Convolution
兩串數列做摺積運算，得到一串數列：一次內積算得一項，第一個數列保持不變，第二個數列頭尾顛倒，並且由第 0 項到第 N-1 項輪流作為首項。有限長數列的摺積，第二個數列會頭尾循環，因而也有人特別稱做循環摺積（ circular convolution ）。
(a0,a1,...,aN-1) convolution (b0,b1,...,bN-1) = (c0,c1,...,cN-1) c0 = (a0,a1,a2,...,aN-1) dot (b0 ,bN-1,bN-2,...,b1) c1 = (a0,a1,a2,...,aN-1) dot (b1 ,b0 ,bN-1,...,b2) c2 = (a0,a1,a2,...,aN-1) dot (b2 ,b1 ,b0 ,...,b3) ⋮ cN-1 = (a0,a1,a2,...,aN-1) dot (bN-1,bN-2,bN-3,...,b0)
(...,a-1,a0,a1,...) convolution (...,b-1,b0,b1,...) = (...,c-1,c0,c1,...) ⋮ c-1 = (...,a-1,a0,a1,...) dot (...,b0 ,b-1,b-2,...) c0 = (...,a-1,a0,a1,...) dot (...,b1 ,b0 ,b-1,...) c1 = (...,a-1,a0,a1,...) dot (...,b2 ,b1 ,b0 ,...) ⋮摺積是線性變換！有限長數列的（循環）摺積，對應的矩陣，正好是循環的 Toeplitz Matrix ；無限長數列的摺積，對應的矩陣，是無限大的矩陣。
(a0,a1,...,aN-1) convolution (b0,b1,...,bN-1) = (c0,c1,...,cN-1) [b0 b-1 b-2 ⋯ bN-1] [a0 ] [c0 ] [b1 b0 b-1 ⋯ bN-2] [a1 ] [c1 ] [b2 b1 b0 ⋯ bN-3] * [a2 ] = [c2 ] [⋮ ⋮ ⋮ ⋮ ] [⋮ ] [⋮ ] [bN-1 bN-2 bN-3 ⋯ b0 ] [aN-1] [cN-1]
(...,a-1,a0,a1,...) convolution (...,b-1,b0,b1,...) = (...,c-1,c0,c1,...) [ ⋮ ⋮ ⋮ ] [⋮ ] [⋮ ] [⋯ b0 b-1 b-2 ⋯] [a-1] [c-1] [⋯ b1 b0 b-1 ⋯] * [a0 ] = [c0 ] [⋯ b2 b1 b0 ⋯] [a1 ] [c1 ] [ ⋮ ⋮ ⋮ ] [⋮ ] [⋮ ]摺積運算的單位元素是脈衝函數。摺積運算擁有交換律、結合律、加法分配律、線性。 Linear Time-Invariant System 皆可寫成摺積。不懂沒關係，本篇用不到。
Z-transform
一串數列做 Z-transform ，得到一個複數多項式：數列每一項乘上複數 z
的次方，次方值是負的索引值，最後全部加起來。可以看作是一串數列與一串複數數列進行內積。z 是一個複數，可以代入數值；不過一般都是維持原本外貌，以代數 z 示人。
Z-transform of (a0,a1,...,aN-1) = a0z0 + a1z-1 + ... + aN-1z-(N-1)
Z-transform of (...,a-1,a0,a1,...) = ... + a-1z1 + a0z0 + a1z-1 + ...
數列摺積與多項式乘法（普通的演算法）
因為電腦無法處理無限長數列，所以此處只討論有限長數列。兩串數列做 Z-transform ，變成多項式，然後相乘一下。因為是有限長數列，我們特意讓次方循環。
a0 * z0 + a1 * z-1 + a2 * z-2
×) b0 * z0 + b1 * z-1 + b2 * z-2
——————————————————————————————————————————————————————————————————
a0b2 * z-2 + a1b2 * z-3 + a2b2 * z-4
a0b1 * z-1 + a1b1 * z-2 + a2b1 * z-3
a0b0 * z0 + a1b0 * z-1 + a2b0 * z-2
——————————————————————————————————————————————————————————————————
a1b2 * z0 + a2b2 * z-1 + a0b2 * z-2
a2b1 * z0 + a0b1 * z-1 + a1b1 * z-2
a0b0 * z0 + a1b0 * z-1 + a2b0 * z-2
——————————————————————————————————————————————————————————————————
(a0,a1,a2) (a0,a1,a2) (a0,a1,a2)
dot * z0 + dot * z-1 + dot * z-2
(b0,b2,b1) (b1,b0,b2) (b2,b1,b0)
把相乘結果做逆向 Z-transform ，竟是原本兩串數列的摺積！數列摺積就是多項式乘法、多項式乘法就是數列摺積，藉由 Z-transform 改變計算順序罷了（遮住 z 不看，就很明顯了）。
數列摺積與多項式乘法的時間複雜度都是 O(N^2) 。
【註：計算學家重視計算，數列摺積與多項式乘法是主角， Z-transform 只是輔助；數學家重視性質， Z-transform
才是主角，數列摺積與多項式乘法只是應用。坊間書籍行文風格傾向數學家；書讀不通的時候，不妨揣摩一下數學家的觀點吧。】
數列摺積與多項式乘法（ Fourier Transform ）
傅立葉轉換其實就是一串數列進行 Z-transform ，變成複數多項式， z 分別代入 e i*2π*(0/N) 、 e
i*2π*(1/N) 、……、 e i*2π*((N-1)/N) ，求得 N 個實際數值。兩個多項式相乘，可以簡化成這 2N 個數值點對點相乘！
正向與逆向傅立葉轉換需時 O(NlogN) ， N 個數值兩兩對應相乘需時 O(N) ，總時間複雜度為 O(NlogN) 。
證明過程猜測如下：線性變換的本質，是在 eigenvector 上面進行倍率縮放，摺積也不例外。傅立葉轉換的矩陣，恰好是 orthogonal
eigenbasis ，所以每個維度就可以各自縮放了。數論轉換也有著相同的性質，哈特利轉換則有著類似的性質。
UVa 12298 ICPC 4671 5705
數列摺積與多項式乘法──循環的 Toeplitz Matrix 的乘法
循環的 Toeplitz Matrix 乘以向量，就是循環摺積。得運用傅立葉轉換、數論轉換。一般的 Toeplitz Matrix 乘以向量，只要填充一下元素，也能變成循環的 Toeplitz Matrix 乘以向量。
數列摺積與多項式乘法──大數乘法
大數乘法是直式乘法，效果同多項式乘法、同數列摺積。得運用傅立葉轉換、數論轉換。大數乘法的演算法有 Schönhage-Strassen Algorithm 與 Fürer's Algorithm ，事實上只是套用不同的數論轉換演算法而已。
注意到，此處談論的多項式乘法，次方是循環的；然而大數乘法，位數並不會循環。解決方式是：被乘數陣列與乘數陣列，預先增加長度，務必大於等於乘積陣列，如此就不受循環影響。
時域摺積等於頻域乘法，頻域摺積等於時域乘法。
傅立葉轉換的輸入稱做「時域」、輸出稱作「頻域」，呼應傅立葉轉換的功能──把波（時間軸）表示成頻譜（頻率軸）。前面段落所講的是時域摺積等於頻域乘法，事實上頻域摺積也等於時域乘法。
Wavelet
（Under Construction!）
（Under Construction!）
程度★★ 難度★
Wavelet Transform
（Under Construction!）
（Under Construction!）
程度★★ 難度★★
No comments:
Post a Comment