Wave - 演算法筆記

振動、振盪是物理學名詞，振動（ Vibration ）是來回運動，振盪（ Oscillation ）是 ... 均勻分布的粒子之中，某個粒子振動所產生的波，剛好也呈現sin 函數，英文稱 ... Wave(ℝ) 日往則月來，月往則日來，日月相推而明生焉。

《易傳》 振動、振盪 這個世界天天都在振動。

地面、空氣、海水、機械、人體等等，都是不斷振動。

振動、振盪是物理學名詞，振動（Vibration）是來回運動，振盪（Oscillation）是來回變化。

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

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

平穩的振動 最平穩的振動，就是高中物理教的簡諧運動：等速圓周運動投影到座標軸，呈sin函數。

sin函數和cos函數長相一樣，只是起點不同而已。

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

振動的快慢：頻率 單位時間振動的次數，稱作「頻率Frequency」。

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

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

振動的高低：振幅 振動的最高（低）距離，稱作「振幅Amplitude」。

人類能感知振幅，受頻率大小影響。

就聽覺而言，振幅高則大聲、振幅低則小聲；就視覺而言，振幅高則亮、振幅低則暗。

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

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

注意到，相位是圓周運動cos函數的角度，而不是振動高低位置。

物理學家喜歡用相位。

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

振動有疊加效果 現實世界當中，多個振動時常融合成一個振動，等於各個振動高低位置相加。

相同方向則增益、相反方向則抵銷。

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

振動有傳遞效果：波 一個粒子振動，就會牽引隔壁粒子振動，一傳十、十傳百。

宏觀之下，形成「波Wave」。

觀察任意一個粒子，都是在振動。

傳遞速度取決於粒子之間的作用力、粒子的質量。

作用力強、質量小，則傳遞速度快。

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

水的高低起伏，就是水波。

空氣的疏密，就是聲波。

地的高低起伏與左右晃動，就是地震波。

電場與磁場的交互作用，就是電磁波。

光波經實驗證明是電磁波。

原子的振動，也許是熱。

有人覺得氣功也許是波，就叫做氣功波。

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

混淆視聽罷了。

輸入　輸出　名稱 離散　離散　DiscreteCosineTransform 離散　連續　似乎沒有名稱 連續　離散　FourierCosineSeries 連續　連續　FourierCosineTransform 離散到離散的餘弦轉換，輸入和輸出都是一串數列。

電腦做運算，數值皆離散。

本文介紹離散版本。

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

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

DiscreteCosineTransform物理意義 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倍，分別乘上振幅，疊加成一個複雜的波。

DiscreteCosineTransform有許多版本 N個平穩的波，微調頻率振幅相位。

維基百科列出了許多版本，大家習慣以第2型作為正向轉換、以第3型作為逆向轉換。

餘弦轉換必須指定相位。

傅立葉轉換可以自動得到相位。

小波轉換可以自由設計波形，不必是平穩的波。

餘弦轉換與其他轉換相比顯得大費周折。

不過有一種說法是：大量資料，分段處理，而微調有助於銜接段落。

圖片壓縮JPEG、影片壓縮H.265和VP9、聲音壓縮MP3和AAC都使用餘弦轉換。

DiscreteCosineTransform數學公式 正向餘弦轉換 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以置中對齊。

DiscreteCosineTransform是線性函數 餘弦轉換是線性函數！可以寫成矩陣形式！ 正向餘弦轉換，視角是橫條點積投影，可以畫成這樣子： 逆向餘弦轉換，視角是直條加權總和，可以畫成這樣子： 演算法（公式解） 依照公式實作，時間複雜度O(N²)。

constintN=10; floatx[N],y[N]; constfloatπ=2.0*acos(0); constfloatθ=2.0*π/N; constfloats=sqrt(2.0/N); //cos()不是O(1)。

//內部迴圈不斷呼叫cos()，因此整體不是O(N²)。

voidDCT() { for(inti=0;i constintN=10; floatx[N],y[N]; constfloatπ=2.0*acos(0); constfloatθ=2.0*π/N; constfloats=sqrt(2.0/N); //三角函數和角公式，求得θ+=dθ之後的三角函數值。

voidadd(float&cosθ,float&sinθ,floatcosdθ,floatsindθ) { floattemp_cosθ=cosθ; cosθ=cosθ*cosdθ-sinθ*sindθ; sinθ=temp_cosθ*sindθ+sinθ*cosdθ; } voidDCT() { constfloatcosdθ=cos(θ/2.0); constfloatsindθ=sin(θ/2.0); constfloatcosdφ=cos(θ/4.0); constfloatsindφ=sin(θ/4.0); floatcosθ=1,sinθ=0; floatcosφ=1,sinφ=0; for(inti=0;i 演算法（Arai–Agui–NakajimaAlgorithm） 針對輸入輸出只有8點的餘弦轉換，進行細部加速。

https://www.nayuki.io/page/fast-discrete-cosine-transform-algorithms IntegerDiscreteCosineTransform 實數運算既複雜又緩慢。

改弦易轍，以整數運算趨近正確答案： 整數運算比實數運算簡單，整數餘弦轉換比餘弦轉換迅速。

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

https://stackoverflow.com/questions/18621167/ 2DDiscreteCosineTransform 餘弦轉換可以推廣到高維度。

二維餘弦轉換，輸入和輸出都是一個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³)。

GibbsPhenomenon 連續到離散的餘弦轉換，有個缺點：先拆解再疊加，產生針刺。

函數曲線劇烈起伏之處（斜率很大或者很小）尤其明顯。

彷彿多項式內插的RungePhenomenon。

Wave(ℂ) 有物混成，先天地生，寂兮寥兮，獨立而不改，周行而不殆，可以為天下母。

《老子》 ComplexNumber 快速複習一下複數吧。

實數，再額外考慮𝑖=√-1，就是複數。

例如2+3𝑖、(1-√2)+(1/3)𝑖、1/(-2𝑖-4)、∛𝑖、sin(𝑖)。

凡是複數，必可重新整理成左邊實數不乘𝑖、右邊實數有乘𝑖，兩個部分相加的格式。

不乘𝑖的部分叫做實部（realpart），有乘i的部分叫做虛部（imaginarypart）。

例如1/(-2𝑖-4)可以重新整理成-0.2+0.1𝑖，其中實部是-0.2，虛部是0.1𝑖。

複數亦可以畫成圖形。

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

兩個複數相加，就是實部加實部、虛部加虛部。

在複數平面上，外觀宛如向量相加。

兩個複數相乘，就是實乘實、虛乘虛、實乘虛、虛乘實，再累加這四個乘積。

在複數平面上，外觀宛如長度相乘、角度相加。

一個複數可以重新表示成一個長度與一個角度，叫做極座標表示法。

長度可以用畢氏定理求得，角度可以用arctan函數求得。

一個長度與一個角度也可以還原成一個複數。

實部可以用cos函數求得，虛部可以用sin函數求得。

附帶一提，長度也有人稱作強度（magnitude），角度也有人稱作相位（phase）。

Euler'sFormula 強者歐拉發現這世界上有一個神奇數字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。

離題了。

Wave 因為e𝑖θ長得像波，所以用e𝑖θ將實數波推廣成複數波。

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

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

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

FourierTransform FourierTransform 「傅立葉轉換」是雙射函數，輸入輸出都是一串複數，可以是離散數列或者連續函數，各有對應名稱。

混淆視聽罷了。

輸入　輸出　名稱 離散　離散　DiscreteFourierTransform 離散　連續　Discrete-timeFourierTransform 連續　離散　FourierSeries 連續　連續　FourierTransform 離散到離散的傅立葉轉換，輸入和輸出都是一串複數數列。

電腦做運算，數值皆離散。

本文介紹離散版本。

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

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

DiscreteFourierTransform物理意義 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倍，分別乘上振幅、添上相位，疊加成一個複雜的波。

DiscreteFourierTransform數學公式 正向傅立葉轉換 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 DiscreteFourierTransform是線性函數 傅立葉轉換是線性函數，恰是正規正交矩陣。

ω=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] 複數波，變成離散數列，可以畫成這樣子： 傅立葉轉換的矩陣，可以畫成這樣子： 演算法（公式解） 依照公式實作，時間複雜度O(N²)。

constintN=10; complexx[N],y[N]; //cos()和sin()不是O(1)。

//內部迴圈不斷呼叫cos()與sin()，因此整體不是O(N²)。

voidDFT() { constfloatπ=2.0*acos(0); constfloatθ=2.0*π/N; for(inti=0;i(cos(θ*i*j),sin(θ*i*j)); } } constintN=10; floatxR[N],xI[N]; //輸入，實部與虛部。

floatyR[N],yI[N]; //輸出，實部與虛部。

//三角函數和角公式，求得θ+=dθ之後的三角函數值。

voidadd(float&cosθ,float&sinθ,floatcosdθ,floatsindθ) { floattemp_cosθ=cosθ; cosθ=cosθ*cosdθ-sinθ*sindθ; sinθ=temp_cosθ*sindθ+sinθ*cosdθ; } voidDFT() { constfloatπ=2.0*acos(0); constfloatθ=2.0*π/N; constfloatcosdθ=cos(-θ); constfloatsindθ=sin(-θ); floatcosθ=1,sinθ=0; for(inti=0;i 演算法（Horner'sRule） 依照公式實作，點積視作多項式，時間複雜度O(N²)。

constintN=10; complexx[N],y[N]; //exp()不是O(1)。

//內部迴圈不斷呼叫exp()，因此整體不是O(N²)。

voidDFT() { constfloatπ=2.0*acos(0); constfloatθ=2.0*π/N; complexim(0,1); //虛數𝑖 for(inti=0;iωn=exp(-im*θ*(float)i); y[i]=0; for(intj=N-1;j>=0;--j) y[i]=y[i]*ωn+x[j]; } } constintN=10; complexx[N],y[N]; voidDFT() { constfloatπ=2.0*acos(0); constfloatθ=2.0*π/N; complexim(0,1); //虛數𝑖 complexω=exp(-im*θ); complexωn=1; for(inti=0;i=0;--j) y[i]=y[i]*ωn+x[j]; ωn*=ω; } } 演算法（Cooley–TukeyAlgorithm） 時間複雜度優於O(N²)的傅立葉轉換演算法，老人家稱作「快速傅立葉轉換FastFourierTransform,FFT」。

這裡介紹最經典的快速傅立葉轉換。

公式的偶數項與奇數項分開整理，採用DynamicProgramming，時間複雜度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的遞推過程，偶數項與奇數項分開處理，索引值不連續，不易取值。

預先重新排列陣列元素，符合遞推過程，減少cachemiss；還可以重複使用記憶體、節省空間。

如何重新排列呢？索引值的二進位表示法，高低位數顛倒之後，恰是正確結果！ 重新排列的時間複雜度是O(N)。

假設高低位數顛倒是O(1)。

顛倒整個過程，不必每回合重算ω。

constfloatπ=2.0*acos(0); constintN=8; complexx[N]; voidFFT() { /*bit-reversalpermutation*/ for(inti=1,j=0;i>1;!((j^=k)&k);k>>=1); // for(intk=N>>1;k>(j^=k);k>>=1); if(i>j)swap(x[i],x[j]); // if(iω(cos(θ),sin(θ)); //每k個做一次FFT for(intj=0;jωn(1,0); for(inti=j;ia=x[i]; complexb=x[i+k/2]*ωn; x[i]=a+b; x[i+k/2]=a-b; ωn*=ω; } } } } //顛倒整個過程 voidFFT() { //不必每回合重算ω floatθ=-2.0*π/N; complexω(cos(θ),sin(θ)); /*dynamicprogramming*/ for(intk=N;k>=2;k>>=1) { for(intj=0;jωn(1,0); for(inti=j;ia=x[i]; complexb=x[i+k/2]; x[i]=a+b; x[i+k/2]=(a-b)*ωn; ωn*=ω; } } //ω每回合翻倍 ω*=ω; } /*bit-reversalpermutation*/ ...... } voidFFT() { floatθ=-2.0*π/N; complexω(cos(θ),sin(θ)); /*dynamicprogramming*/ for(intk=N;k>=2;k>>=1) { //對調內外迴圈，讓ωn少乘幾次。

//缺點則是索引值更容易跳動，更容易產生cachemiss。

complexωn(1,0); for(inti=0;ia=x[j]; complexb=x[j+k/2]; x[j]=a+b; x[j+k/2]=(a-b)*ωn; } ωn*=ω; } ω*=ω; } /*bit-reversalpermutation*/ ...... } HartleyTransform 哈特利轉換是雙射函數，輸入和輸出都是一串實數。

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

傅立葉轉換： 2πft2πft-𝑖2πft/N cos————-𝑖⋅sin————=e NN 哈特利轉換： 2πft2πft2πft cos————+sin————=cas———— NNN 另一個哈特利轉換，比較沒人用： 2πft2πft2πft cos————-sin————=cis———— NNN 傅立葉轉換： 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-ConquerMethod。

不一樣的是奇數項的處理方式，提出常數的步驟變複雜了。

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://home.iae.nl/users/mhx/fft.html 2DDiscreteFourierTransform 傅立葉轉換可以推廣到高維度。

二維傅立葉轉換，輸入輸出都是一個N×N複數方陣。

輸入方陣，分別除以N×N種二維複數波，得到N×N個輸出數值，形成輸出方陣。

（由於二維複數波很難畫，以下改畫二維實數波。

） 依照公式實作，時間複雜度O(N⁴)。

高速演算法是每一橫條各自傅立葉轉換，然後每一直條各自傅立葉轉換，時間複雜度O(NNlogN+NNlogN)=O(N²logN)。

FourierTransform的性質 FrequencySpectrum 傅立葉轉換，輸出數列有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倍，頻譜呈現「人」型，所有頻率皆有強度，漏得到處都是。

這個現象稱作「spectralleakage」。

這個現象有兩種解讀： 一、1倍和2倍頻率波，疊加之後，結果是1倍，不是1.5倍。

更明確來說是最大公因數。

傅立葉轉換的0倍頻率波到N-1倍頻率波，皆無法組合出1.5倍。

只好湊合各種頻率，盡量趨近1.5倍。

二、離散版本的傅立葉轉換，輸入輸出是循環數列。

1.5倍，循環之後，其實不是平穩的振動，因而產生許多高頻波。

當振動次數不是整數次、頻率不是整數倍，那麼傅立葉轉換無法精準量測！這是重大缺點！ WindowFunction 然而數學家尚未發明更好的方式。

當今主流仍是傅立葉轉換。

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

原本數列，乘上一個窗函數：中央高、兩端趨近零的數列。

如此令原本數列左右兩端連續，抑制頻譜多餘強度。

窗函數非常多種，功效略有差異。

請讀者自行研究。

voidspectrum(floatarray[],intN) { //Hannwindow constfloatπ=3.1415926f; floatwindow[N]; for(inti=0;i WindowFunction的快速演算法 傅立葉轉換，有時候輸入稱作「時域TimeDomain」、輸出稱作「頻域FrequencyDomain」，呼應傅立葉轉換的功能：把波（時間軸）表示成頻譜（頻率軸）。

時域乘法等於頻域循環卷積。

請見本站文件「Filter」。

數列與窗函數相乘，等於數列與窗函數在頻域的循環卷積。

窗函數，多由cos波組成；窗函數在頻域，只有少數幾點有值。

例如Hann窗，從時域轉頻域，只有三點有值。

因此，與其在時域套用窗函數，不如在頻域套用窗函數。

過程非常簡單：每個值減去兩側的值（相位差不多是π），附帶權重。

這揭露了窗函數的真正功效──頻譜之中，平者更平，尖者更尖。

voidspectrum(floatarray[],intN) { //discretefouriertransform complexf[N]; fft(array,f); //HannwindowwhenN>=4 //f'[i]=(1/2*f[i])-(1/4*f[i+1])-(1/4*f[i-1]) for(inti=0;ix=f[i]*0.5 -f[(i-1+N)%N]*0.25 -f[(i+1)%N]*0.25; magnitude[i]=abs(x); //sqrtlength phase[i]=arg(x); //atan } } 最後額外補充一下。

連續版本的傅立葉轉換，窗函數頻譜，外觀是一個尖峰。

取abs和log，外觀是一個大圓丘（mainlobe），附帶連綿小矮丘（sidelobe）。

很多資工系老師上課只教連續版本，但是我們根本不會用到連續版本！ 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波的加權總和。

傅立葉轉換是線性函數。

換句話說，輸入數列們的加權總和，經過傅立葉轉換，等於輸出數列們的加權總和；但是不等於頻譜們的加權總和！ 輸出數列是複數。

複數加法是向量相加，複數倍率是向量伸縮。

向量相加不等於長度相加、角度相加。

（唯一例外：所有波都是整數次的平穩振動。

因為頻譜幾乎都是零。

） 多個波形疊加，不會正確反映於頻譜！這是重大缺點！ 然而大家仍用頻譜分解頻率，無法可管。

自己保重。

SparseFourierTransform 只計算特定頻率的強度與相位。

速度較快。

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}] FourierTransform的性質 輸入輸出對應 連續到連續的傅立葉轉換，輸入輸出有著特殊對應關係。

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

運算對應 加法-加法a+b=aft+bft 倍率-倍率a⋅k=aft⋅k 乘法-卷積a×b=aft∗bft 卷積-乘法a∗b=aft×bft 微分-角加速度a′=aft⋅2πifd/dta(t)=2πif⋅aft(f) 角加速度-微分a⋅2πit=aft′ 平方和(能量)守恆‖a‖²=‖aft‖²∑a(t)²=∑aft(f)² 因為正向轉換幾乎等於逆向轉換，所以運算對調之後，對應依然成立。

LaplaceTransform 拉普拉斯轉換是傅立葉轉換的推廣版本。

有兩個地方不同： 一、e-𝑖⋅(2π/N)⋅f⋅t的次方值，改成任意複數。

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

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

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

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

計算學家不使用拉普拉斯轉換。

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

Wavelet（UnderConstruction!） Wavelet 「小波」，自訂特殊造型的波。

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

哪些小波可以用於小波轉換呢？以線性代數的觀點來看，N個波必須構成N維空間，才有反矩陣。

換句話說，線性獨立導致雙射函數。

餘弦轉換、傅立葉轉換，N個波線性獨立，是雙射函數。

哪些小波可以用於小波轉換呢？以線性代數的觀點來看，最簡潔的方式是正規正交基底，反矩陣就是轉置矩陣。

正規是指個個範數為1。

正交是指兩兩內積為0。

換句話說，N個小波能量皆是1（一個小波的N個位置的平方的總和是1），N個小波互相垂直。

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

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

https://en.wikipedia.org/wiki/Fast_wavelet_transform https://www.jjj.de/fxt/fxtbook.pdf