卷積與線性卷積：從周期延拓到FFT補(bǔ)零實(shí)現(xiàn))
搞數(shù)字信號處理的沒被循環(huán)卷積和線性卷積折磨過那基本不太正常。尤其是學(xué)到DFT那一章書上突然冒出來一句“循環(huán)卷積可以通過DFT計算”緊跟著又說“線性卷積可以通過補(bǔ)零然后用循環(huán)卷積實(shí)現(xiàn)”很多人直接就暈了兩個東西明明長得不一樣怎么繞一圈又等價了這不矛盾嗎我當(dāng)時學(xué)這塊也卡了很久。后來自己親手把兩個序列的卷積結(jié)果列出來補(bǔ)零、延拓、對位、相加一步步走一遍才發(fā)現(xiàn)這件事其實(shí)特別直觀就是“周期延拓之后發(fā)生混疊”和“補(bǔ)零補(bǔ)夠長度讓混疊自然消失”的區(qū)別而已。搞明白這個DFT做快速卷積、分段卷積這些實(shí)用技巧就都通了。這篇文章就圍繞“循環(huán)卷積和線性卷積的關(guān)系”來寫適合剛學(xué)完DFT、正在做課程設(shè)計或者是被“為什么FFT卷積要先補(bǔ)零”這個問題困擾的同學(xué)。我會從兩者定義講起配合具體的數(shù)值演算把等價條件為什么是 L ≥ N1 N2 - 1 講透再補(bǔ)充快速卷積的實(shí)現(xiàn)步驟和工程里常見的坑。1. 先搞懂兩個卷積的底層邏輯1.1 線性卷積我們最熟悉的“反折-移位-相乘-累加”線性卷積的定義我不多念書上的公式了直接給個直觀理解。假設(shè)有兩個有限長序列一個長度 N一個長度 M那么它們的線性卷積結(jié)果長度是 N M - 1。怎么算出來的就是拿一個序列逐點(diǎn)滑動和另一個序列相乘再累加滑一次出一個點(diǎn)。舉例來說取 x[n] [1, 2, 3]h[n] [4, 5]它倆的線性卷積用手寫一遍y[0] 1×4 4y[1] 1×5 2×4 13y[2] 2×5 3×4 22y[3] 3×5 15所以結(jié)果 y [4, 13, 22, 15]長度是 3 2 - 1 4。線性卷積的物理意義很明確它對應(yīng)的是“輸入信號通過一個線性時不變系統(tǒng)的完整響應(yīng)”每個輸出點(diǎn)都是系統(tǒng)對所有歷史輸入的記憶疊加。濾波、系統(tǒng)響應(yīng)、信號通過信道這些場景下我們通常要的都是線性卷積。1.2 循環(huán)卷積周期延拓以后的對位相乘累加循環(huán)卷積的定義看起來也很簡單只是把“反折-移位”里的移位換成了“循環(huán)移位”。所謂循環(huán)移位就是平移之后超出序列范圍的元素不從另一邊補(bǔ)回來而是繞回去像在一個圓環(huán)上轉(zhuǎn)圈。設(shè)兩個序列長度都為 L循環(huán)卷積要求兩個序列按同一個長度 L 來算那么循環(huán)卷積定義為y[k] Σ_{n0}^{L-1} x[n] · h[(k-n) mod L]重點(diǎn)在 (k-n) mod L這就是循環(huán)的意思。算出來的結(jié)果長度也是 L不會變長。所有下標(biāo)都是對 L 取模之后的值所以序列被看成是長度為 L 的周期序列。還是用 x[n] [1, 2, 3]h[n] [4, 5] 來舉例但做循環(huán)卷積時兩個序列都得先擴(kuò)展成同樣長度。假設(shè)取 L 4那就要給 x 補(bǔ)一個零、給 h 補(bǔ)兩個零x [1, 2, 3, 0]h [4, 5, 0, 0]按循環(huán)卷積公式一項(xiàng)項(xiàng)算y[0] x[0]h[0] x[1]h[3] x[2]h[2] x[3]h[1] 1×4 2×5 3×4 0×5 26y[1] x[0]h[1] x[1]h[0] x[2]h[3] x[3]h[2] 1×5 2×4 3×5 0×4 28y[2] x[0]h[2] x[1]h[1] x[2]h[0] x[3]h[3] 1×4 2×5 3×4 0×5 26y[3] x[0]h[3] x[1]h[2] x[2]h[1] x[3]h[0] 1×5 2×4 3×5 0×4 28結(jié)果是 [26, 28, 26, 28]跟線性卷積的 [4, 13, 22, 15] 完全對不上。問題就出在循環(huán)卷積里 h 的下標(biāo)是取模得到的比如 h[(0-1) mod 4] 取到的其實(shí)是 h[3]而在原來的線性卷積里 h[3] 根本不存在它應(yīng)該是補(bǔ)零后的點(diǎn)不是繞回來的值。這就是循環(huán)卷積和線性卷積最本質(zhì)的區(qū)別線性卷積干干凈凈地算完結(jié)果變長循環(huán)卷積把超出邊界的部分“繞回來”疊加到了前面的輸出點(diǎn)上導(dǎo)致結(jié)果被“污染”。1.3 為什么DFT天然對應(yīng)循環(huán)卷積這里要給剛學(xué)到這兒的同學(xué)提一個重要結(jié)論時域循環(huán)卷積對應(yīng)頻域乘積即 DFT(x ? h) X(k) · H(k)。同樣的頻域循環(huán)卷積對應(yīng)時域乘積。但很多人容易忽略一個點(diǎn)DFT處理的是有限長序列可它在數(shù)學(xué)上有一個隱含的“周期延拓”假設(shè)。教科書上會說“DFT是DFS取主值區(qū)間”意思就是你先想象這個序列在無限長周期延拓然后取其中一個周期做變換。既然信號是周期的那么時域卷積自然也是周期的也就是循環(huán)卷積。所以當(dāng)你直接用 FFT 做兩個序列的頻域相乘再逆變換回來得到的結(jié)果本質(zhì)上是循環(huán)卷積不是線性卷積。工程當(dāng)中大多數(shù)場景要的是線性卷積這就直接導(dǎo)致我們在用 FFT 之前必須做一些處理把循環(huán)卷積的結(jié)果“修正”成線性卷積。2. 核心關(guān)系兩者什么時候相等為什么是這個條件2.1 循環(huán)卷積是怎么“混疊”線性卷積結(jié)果的上面那個例子已經(jīng)看到現(xiàn)象了用 L 4 做循環(huán)卷積結(jié)果跟線性卷積完全不是一回事。那如果把 L 逐步變大會發(fā)生什么如果取 L 5那么兩個序列補(bǔ)零成 x [1, 2, 3, 0, 0]h [4, 5, 0, 0, 0]。此時再按循環(huán)卷積公式算y[0] 1×4 4y[1] 1×5 2×4 13y[2] 2×5 3×4 22y[3] 3×5 15y[4] 0結(jié)果是 [4, 13, 22, 15, 0]。你看前面 4 個點(diǎn)跟線性卷積一模一樣就多了一個 0。為什么因?yàn)楫?dāng) L 足夠大之后循環(huán)移位 (k-n) mod L 取到的那些“繞回”位置都正好落在補(bǔ)零的區(qū)域內(nèi)。也就是說原本會從尾部繞到頭部造成混疊的那些項(xiàng)現(xiàn)在全部乘的是 0混疊也就隨之消失了。換個角度想線性卷積結(jié)果長度為 NM-1只要循環(huán)長度 L 不小于這個長度那么線性卷積結(jié)果的 NM-1 個點(diǎn)可以被完整放進(jìn)一個周期里不會出現(xiàn)前后周期相互重疊的情況。反之如果 L NM-1周期延拓后的各個周期就會“疊”在一起重疊部分的數(shù)值加到一起導(dǎo)致輸出跟線性卷積不同。2.2 條件推導(dǎo)L ≥ NM-1 是怎么來的這里可以很直觀地推一下。線性卷積的非零范圍是x 的非零范圍0 ≤ n ≤ N-1h 的非零范圍0 ≤ n ≤ M-1卷積結(jié)果 y[n] 的非零范圍0 ≤ n ≤ NM-2總長度 NM-1。如果用循環(huán)卷積長度 L 來做兩個序列都補(bǔ)零到 L 點(diǎn)。循環(huán)卷積在取模移位時如果它訪問到的 h[(k-n) mod L] 對應(yīng)的是補(bǔ)零區(qū)間就不會產(chǎn)生混疊如果它訪問到的是 h 原本的非零區(qū)間才是正常貢獻(xiàn)。判斷混疊是否發(fā)生關(guān)鍵看線性卷積結(jié)果的“尾巴”會不會在周期延拓時再退回到頭部。線性卷積結(jié)果是一個長度為 NM-1 的有限長序列現(xiàn)在要以 L 為周期延拓。做周期延拓就是把它不斷復(fù)制平移每個周期的位置相隔 L。如果 NM-1 ≤ L那么相鄰兩個周期的結(jié)果剛好首尾相鄰不會重疊取主值區(qū)間后就是完整的線性卷積結(jié)果。如果 NM-1 L那相鄰兩個周期就會有長度為 NM-1-L 的重疊段重疊部分的數(shù)值必須相加這在時域上等價于“尾部折疊疊加到頭部”。所以條件就是循環(huán)卷積長度 L 必須滿足 L ≥ NM-1此時循環(huán)卷積的結(jié)果在 0 到 NM-2 這些點(diǎn)上等于線性卷積其余點(diǎn)為零。這也是為什么用 DFT 做卷積時通常會取 L NM-1 或者取一個不小于它的 2 的冪方便 FFT這個 2 的冪還要大于等于 NM-1。2.3 改變循環(huán)卷積長度的本質(zhì)是改變周期大小這個點(diǎn)可以再深挖一下對真正理解 DFT 卷積特別有幫助。循環(huán)卷積長度 L 一拍腦袋隨便取本質(zhì)上是人為設(shè)定周期延拓的周期。同一個序列周期越短相鄰周期之間的重疊越嚴(yán)重混疊加的項(xiàng)越多。比如同樣是 x[n] [1, 2, 3]h[n] [4, 5]你取 L 3兩個序列各只有 3 個點(diǎn)循環(huán)卷積算出來是y[0] x[0]h[0] x[1]h[2] x[2]h[1] 1×4 2×5 3×4 26這里 h[2] 取的是 h[0] 繞回h[1] 是正常的y[1] x[0]h[1] x[1]h[0] x[2]h[2] 1×5 2×4 3×5 28y[2] x[0]h[2] x[1]h[1] x[2]h[0] 1×4 2×5 3×4 26結(jié)果是 [26, 28, 26]和 L4 時的前三個值一樣因?yàn)檎嬲斐刹町惖木褪俏膊坷@回的那一項(xiàng)。取 L3 時尾部繞回的項(xiàng)更多、疊加得更嚴(yán)重整段結(jié)果都亂了。我自己的理解是循環(huán)卷積的長度 L 決定了周期有多長周期越長線性卷積結(jié)果延拓后互相疊的概率越小L 只要超過結(jié)果總長度就“一點(diǎn)重疊都不?！贝藭r循環(huán)卷積和線性卷積完全一致。所以你不需要記住復(fù)雜的數(shù)學(xué)證明記住“周期延拓重疊則混疊”這八個字就夠了。3. 實(shí)戰(zhàn)用循環(huán)卷積實(shí)現(xiàn)線性卷積的完整流程3.1 核心步驟補(bǔ)零、DFT、相乘、IDFT既然 DFT 做的是循環(huán)卷積那要實(shí)現(xiàn)線性卷積思路就一句話把循環(huán)卷積的周期拉長長到不重疊為止。假設(shè) x[n] 長度 Nh[n] 長度 M目標(biāo)是算 y[n] x[n] * h[n]線性卷積。標(biāo)準(zhǔn)做法如下計算最小循環(huán)長度L_min N M - 1。選實(shí)際 FFT 長度 L為了用 FFT 加速一般取 L 2^ceil(log2(NM-1))即不小于 L_min 的最接近的 2 的冪。當(dāng)然如果直接用 DFTL 取任意不小于 L_min 的數(shù)都行。給 x[n] 后面補(bǔ) L-N 個零給 h[n] 后面補(bǔ) L-M 個零讓兩個序列長度都變成 L。分別做 L 點(diǎn) DFT得到 X[k] 和 H[k]。逐點(diǎn)相乘 Y[k] X[k] * H[k]。對 Y[k] 做 L 點(diǎn) IDFT得到 y[n]。取 y[0] 到 y[NM-2] 這 NM-1 個點(diǎn)作為最終線性卷積結(jié)果后面從 y[NM-1] 開始應(yīng)該是 0數(shù)值上有浮點(diǎn)誤差時會是非常小的數(shù)。第 3 步是這個方法里最“靈魂”的一步。很多新手犯的錯就是“我已經(jīng)補(bǔ)零到兩個序列一樣長了”但實(shí)際上補(bǔ)的零根本不夠。比如 N1000M50有人直接把兩個序列都補(bǔ)零到 1000 點(diǎn)就去做 FFT結(jié)果出來前一段是對的、后一段全亂了就是因?yàn)?L1000 小于 NM-11049尾部混疊已經(jīng)發(fā)生。3.2 一個可以照著跑的代碼示例為了把流程固定下來這里給一段 Python 代碼核心邏輯用 numpy 實(shí)現(xiàn)并把每一步都注釋清楚。這段代碼可以直接復(fù)制到一個 .py 文件里運(yùn)行。import numpy as np # 兩個測試序列 x np.array([1.0, 2.0, 3.0]) h np.array([4.0, 5.0]) N len(x) M len(h) # 1. 最少需要的循環(huán)卷積長度 L_min N M - 1 # 2. 實(shí)際FFT長度為了加速取不小于L_min的2的冪 L 1 while L L_min: L 1 # 3. 補(bǔ)零到長度L x_pad np.zeros(L) x_pad[:N] x h_pad np.zeros(L) h_pad[:M] h # 4. 各自FFT X np.fft.fft(x_pad) H np.fft.fft(h_pad) # 5. 頻域逐點(diǎn)相乘 Y X * H # 6. 逆變換回時域 y_ifft np.fft.ifft(Y).real # 理論上應(yīng)該是實(shí)數(shù)浮點(diǎn)誤差會產(chǎn)生極小的虛部 # 7. 取前NM-1個點(diǎn)作為線性卷積結(jié)果 y_lin y_ifft[:L_min] # 對比用numpy直接算的線性卷積 y_conv np.convolve(x, h) print(直接線性卷積: , y_conv) print(FFT循環(huán)卷積: , y_lin) print(誤差: , np.max(np.abs(y_lin - y_conv)))這段代碼跑出來誤差基本在 1e-14 這個量級可以認(rèn)為是數(shù)值誤差不是算法錯誤。很多人會問為什么 ifft 之后要取 .real因?yàn)楦↑c(diǎn)運(yùn)算里 FFT/IFFT 產(chǎn)生的虛部誤差雖然極小但你如果不取實(shí)部后面看數(shù)據(jù)時總覺得哪里不對勁甚至有人直接把虛部當(dāng)成錯誤信號。實(shí)際上對于實(shí)數(shù)信號IDFT 的理論結(jié)果是實(shí)數(shù)所以取實(shí)部是合理且常規(guī)的操作。3.3 FFT卷積到底快在哪什么時候該用直接線性卷積的復(fù)雜度是 O(NM)因?yàn)槊總€輸出點(diǎn)要做 M 次乘法累加總共 NM-1 個輸出點(diǎn)。FFT 方法復(fù)雜度是 O(L log L)其中 L 取 2 的冪。補(bǔ)零后 L 約等于 NM所以復(fù)雜度近似 O((NM) log(NM))。當(dāng) N 和 M 都大時FFT 方法的優(yōu)勢非常明顯。比如 NM1024直接卷積大約要做 100 萬次乘加而 L2048 的 FFT 大約是 2048×11 ≈ 22528 次蝶形運(yùn)算量差距接近 50 倍。信號越長差距越夸張。但要注意如果其中一個序列特別短比如濾波器長度只有 8 個點(diǎn)輸入也只有 64 個點(diǎn)那這種短序列直接用 conv 反而更快。因?yàn)?FFT 有固定開銷包含補(bǔ)零、正變換兩次、反變換一次、頻域復(fù)數(shù)乘法光這些調(diào)用本身的成本就超過直接卷積了。我實(shí)測過閾值大概在 N、M 都超過 30~50 以后FFT 才開始體現(xiàn)優(yōu)勢。工程上到底用哪種可以先算一下量級再決定別迷信“FFT 一定快”。4. 工程里的延伸長序列濾波和分段卷積4.1 一個輸入無限長濾波器有限長怎么搞實(shí)際工程里經(jīng)常遇到這種情況濾波器 h[n] 長度固定為 M但輸入信號 x[n] 非常長甚至是一個流式信號不能等全部收到再做 FFT。比如實(shí)時音頻處理、雷達(dá)回波處理信號一直在來內(nèi)存也存不下整個序列。這時候就需要分段卷積。分段卷積的基本思路是把長輸入切成長度 L_block 的小塊每塊分別和 h[n] 做線性卷積然后把各塊的卷積結(jié)果按重疊部分相加。這就是重疊相加法overlap-add。做法如下取塊長 L_block要求 L_block M - 1 不超過 FFT 長度通常直接讓 FFT 長度 L L_block M - 1。對每一塊 x_i[n] 補(bǔ)零后做 FFT乘以 H[k]H 是 h 補(bǔ)零后的 FFT只需要算一次再 IFFT。每塊的輸出長度是 L_block M - 1相鄰塊的輸出會有 M-1 個點(diǎn)的重疊。將當(dāng)前塊輸出與上一次保留下來的 M-1 個點(diǎn)相加再輸出有效的前 L_block 個點(diǎn)后 M-1 個點(diǎn)留給下一次疊加。這個過程用代碼實(shí)現(xiàn)并不復(fù)雜但有點(diǎn)繞。核心點(diǎn)就一個每塊獨(dú)立卷積的結(jié)果在重疊區(qū)要相加而不是直接覆蓋。與之對偶的方法是重疊保留法overlap-save它不重疊相加而是重疊保留輸入丟棄輸出中混疊的那部分。兩種方法效果一樣只是實(shí)現(xiàn)細(xì)節(jié)和邊界處理習(xí)慣不同。分段卷積是 FFT 卷積真正發(fā)揮價值的地方因?yàn)槿绻w做 FFT需要等信號全部到齊這在實(shí)時系統(tǒng)里根本不可能。4.2 頻域?yàn)V波其實(shí)也是循環(huán)卷積思想自適應(yīng)濾波、OFDM 調(diào)制解調(diào)、聲學(xué)回聲消除這類系統(tǒng)里塊處理方式非常常見。它們的共同點(diǎn)是先把時域信號分塊轉(zhuǎn)頻域處理再轉(zhuǎn)回時域。如果塊長度和濾波器長度考慮不當(dāng)就很容易出現(xiàn)“看起來處理了但數(shù)據(jù)對不上”的詭異現(xiàn)象。舉個例子OFDM 里為了保護(hù)多徑時延擴(kuò)展會在符號之間插入循環(huán)前綴。循環(huán)前綴的本質(zhì)是一種循環(huán)延拓它讓線性卷積的信道響應(yīng)變成循環(huán)卷積的形式這樣接收端就能直接用頻域均衡做單抽頭補(bǔ)償。這里面“線性卷積變循環(huán)卷積”的思路跟本文講的補(bǔ)零恰好是相反方向的操作但都是同一個思想在不同場景下的應(yīng)用。另一個例子是很多音頻處理教材里提到的“頻域卷積”實(shí)驗(yàn)把一段音頻和某個房間沖激響應(yīng)做卷積模擬混響效果。如果你不補(bǔ)零直接對兩個信號做 FFT 相乘再 IFFT出來的音頻會有一個很明顯的“回環(huán)”雜音這就是循環(huán)卷積的尾部混疊在音頻上聽起來像回聲短路一樣其實(shí)也就是周期延拓導(dǎo)致的時域 aliasing。4.3 圓周卷積譜與快速卷積理論的統(tǒng)一視角如果仔細(xì)琢磨可以感覺到這里面的理論閉環(huán)DFT 是傅里葉變換的離散化周期化版本它處理的所有序列都隱含周期延拓因此在 DFT 域做卷積天然是循環(huán)卷積而線性卷積是物理世界對“有限長輸入經(jīng)過有限長沖激響應(yīng)系統(tǒng)”的描述。要在 DFT 域?qū)崿F(xiàn)物理世界的線性卷積就必須通過補(bǔ)零擴(kuò)大周期消除周期重疊。補(bǔ)零這個操作看似簡單實(shí)質(zhì)是把“循環(huán)卷積”的周期拉大到足以承載“線性卷積”的結(jié)果長度讓兩個數(shù)學(xué)對象在同一個計算平臺上統(tǒng)一起來。這也是很多專業(yè)課程里把 DFT、循環(huán)卷積、快速卷積、重疊相加法串在一起講的原因。理解了這條主線遇到相關(guān)問題時你腦子里會有一個統(tǒng)一的判斷標(biāo)準(zhǔn)這個系統(tǒng)里線性卷積結(jié)果會不會發(fā)生周期重疊重疊了怎么處理是補(bǔ)零防混疊還是故意利用循環(huán)卷積的特性5. 學(xué)習(xí)與實(shí)操中常見的坑5.1 補(bǔ)零長度只補(bǔ)到兩個輸入一樣長這個錯誤簡直太典型了。比如 x 長度 500h 長度 200很多人覺得“我把兩個都補(bǔ)零到 512然后用 512 點(diǎn) FFT”看起來挺對稱。但注意線性卷積長度是 699512 根本不夠。結(jié)果就是卷出來的前 500 個點(diǎn)可能還能看從第 500 個點(diǎn)開始后續(xù)響應(yīng)全部消失尾部混疊值直接疊在頭部附近整體結(jié)果和預(yù)期差一大截。正確做法是先算 NM-1再決定 FFT 長度。哪怕不用 FFT、直接調(diào)庫函數(shù)也要對輸出長度的預(yù)期心里有數(shù)。5.2 忘記取 IFFT 的前 NM-1 個點(diǎn)補(bǔ)零補(bǔ)夠了FFT 也做了IFFT 也做了結(jié)果發(fā)現(xiàn)輸出數(shù)組長度為 L而 L 往往大于預(yù)期的 NM-1。如果直接把整個數(shù)組拿去跟線性卷積結(jié)果比后面多出來的一串近似為 0 的數(shù)會影響誤差判斷有時候甚至有人把多出來的零當(dāng)成有效信號導(dǎo)致后面處理數(shù)組長度不匹配。習(xí)慣做法是IFFT 后只取前 NM-1 個點(diǎn)。多出來的點(diǎn)要么是非常小的浮點(diǎn)噪聲要么是嚴(yán)格為 0。處理信號數(shù)組時長度不匹配的 bug 很難排查最好一開始就指定輸出區(qū)間。5.3 實(shí)數(shù)序列做完 FFT 后埋著復(fù)數(shù)虛部不管很多信號是實(shí)數(shù)FFT 之后變成復(fù)數(shù)IFFT 回來理論上是實(shí)數(shù)但由于浮點(diǎn)誤差虛部會有 1e-15 量級的殘留。如果不取實(shí)部后面繪制波形圖時會看到一些莫名其妙的極小虛部雖然不影響幅值但有些工具里會警告“數(shù)據(jù)是復(fù)數(shù)”甚至導(dǎo)致畫圖報錯。建議在 IFFT 后用 .real 取實(shí)部。但也別用全盤取絕對值的方式處理因?yàn)樾盘柪锟赡苡胸?fù)值取絕對值會把負(fù)半軸整體翻上去破壞波形。5.4 下標(biāo)從 1 還是從 0最容易把人繞暈MATLAB 里索引從 1 開始Python 里從 0 開始。線性卷積定義里 y[n] 的 n 也是從 0 開始標(biāo)。寫代碼的時候補(bǔ)零位置、取主值區(qū)間、分段卷積的塊索引這些地方但凡少一個 1 或者多一個 -1結(jié)果就可能整體偏移一位。而且這種偏移在波形圖上肉眼幾乎看不出來只有和參考結(jié)果做數(shù)值對比時才能發(fā)現(xiàn)。我的經(jīng)驗(yàn)是寫卷積相關(guān)代碼時不要直接看公式抄先把一個特別小的測試序列跑起來比如 x[1,2]h[3,4]手算出線性卷積結(jié)果然后程序跑一遍逐步打印中間過程來核對索引。5.5 用頻域相乘實(shí)現(xiàn)卷積時忘了 H 也要補(bǔ)零有時候 x 補(bǔ)零了H 的 FFT 卻用的是原始長度或者反過來。這會造成 X[k] 和 H[k] 的頻譜點(diǎn)數(shù)不一致FFT 長度不匹配時直接相乘會報錯或者更隱蔽的是長度碰巧一致但補(bǔ)零規(guī)則不一樣比如 X 補(bǔ)到了 1024H 只補(bǔ)到了 512然后廣播相乘結(jié)果完全沒有意義。所有參與 DFT 的序列必須補(bǔ)到同一個長度 L。這里沒有例外。5.6 把 FFT 卷積用在小規(guī)模序列上反而更慢這一點(diǎn)是性能上的坑。FFT 卷積雖然在大規(guī)模下有巨大優(yōu)勢但小規(guī)模下并不劃算。比如兩個長度都只有 8 的序列直接用雙重循環(huán) 64 次乘法就夠了FFT 反而要做補(bǔ)零、64 點(diǎn) FFT、復(fù)數(shù)乘法、IFFT耗時可能是直接卷積的 5 到 10 倍。我在實(shí)際項(xiàng)目里一般會設(shè)置一個閾值比如 N*M 2000 時才切換 FFT 方案低于閾值直接暴力卷積整體性能最優(yōu)。5.7 驗(yàn)證結(jié)果時只看“像不像”不看數(shù)值做實(shí)驗(yàn)最忌諱的就是輸出畫出來“看起來差不多”就交差。循環(huán)卷積與線性卷積的差異不是噪聲而是結(jié)構(gòu)性的。如果你補(bǔ)零長度差得不多混疊可能只出現(xiàn)在特定區(qū)間圖像整體趨勢接近但個別點(diǎn)偏差很大。這種問題肉眼不一定看得出來但數(shù)值一對比就暴露了。我自己的習(xí)慣是每次跑通 FFT 卷積后用 np.convolve 或者 MATLAB 的 conv 函數(shù)做一次對照計算 max(abs(y_fft - y_conv))只要這個值在 1e-10 以下說明實(shí)現(xiàn)沒問題。之后再去做長信號的工程實(shí)現(xiàn)心里就有底了。最后說點(diǎn)實(shí)際的循環(huán)卷積和線性卷積的關(guān)系說到底是數(shù)學(xué)在“有限長”這件事上的一次精妙折中DFT 只能處理有限長序列但它骨子里又是周期性的所以卷積必須按循環(huán)來做而物理世界大量場景需要的線性卷積恰好可以通過把周期拉長來“假裝”實(shí)現(xiàn)。補(bǔ)零不是為了好看也不是為了湊 2 的冪而是為了給線性卷積結(jié)果騰出一個不重疊的周期空間。我在實(shí)驗(yàn)室里第一次自己寫 FFT 卷積濾波器的時候就是沒把補(bǔ)零長度當(dāng)回事用的 512 點(diǎn) FFT 去處理 300 點(diǎn)信號和 200 點(diǎn)沖激響應(yīng)結(jié)果輸出從第 113 個點(diǎn)左右開始完全變形。后來畫圖才發(fā)現(xiàn)是尾部周期混疊一下就想通了“L ≥ NM-1”這個條件到底在防什么東西。從那以后再遇到頻域卷積、分段濾波我第一件事就是算清楚長度再開始寫代碼。希望這篇文章也能讓你少走這個彎路。 循環(huán)卷積與線性卷積：從周期延拓到FFT補(bǔ)零實(shí)現(xiàn)
搞數(shù)字信號處理的沒被循環(huán)卷積和線性卷積折磨過那基本不太正常。尤其是學(xué)到DFT那一章書上突然冒出來一句“循環(huán)卷積可以通過DFT計算”緊跟著又說“線性卷積可以通過補(bǔ)零然后用循環(huán)卷積實(shí)現(xiàn)”很多人直接就暈了兩個東西明明長得不一樣怎么繞一圈又等價了這不矛盾嗎我當(dāng)時學(xué)這塊也卡了很久。后來自己親手把兩個序列的卷積結(jié)果列出來補(bǔ)零、延拓、對位、相加一步步走一遍才發(fā)現(xiàn)這件事其實(shí)特別直觀就是“周期延拓之后發(fā)生混疊”和“補(bǔ)零補(bǔ)夠長度讓混疊自然消失”的區(qū)別而已。搞明白這個DFT做快速卷積、分段卷積這些實(shí)用技巧就都通了。這篇文章就圍繞“循環(huán)卷積和線性卷積的關(guān)系”來寫適合剛學(xué)完DFT、正在做課程設(shè)計或者是被“為什么FFT卷積要先補(bǔ)零”這個問題困擾的同學(xué)。我會從兩者定義講起配合具體的數(shù)值演算把等價條件為什么是 L ≥ N1 N2 - 1 講透再補(bǔ)充快速卷積的實(shí)現(xiàn)步驟和工程里常見的坑。1. 先搞懂兩個卷積的底層邏輯1.1 線性卷積我們最熟悉的“反折-移位-相乘-累加”線性卷積的定義我不多念書上的公式了直接給個直觀理解。假設(shè)有兩個有限長序列一個長度 N一個長度 M那么它們的線性卷積結(jié)果長度是 N M - 1。怎么算出來的就是拿一個序列逐點(diǎn)滑動和另一個序列相乘再累加滑一次出一個點(diǎn)。舉例來說取 x[n] [1, 2, 3]h[n] [4, 5]它倆的線性卷積用手寫一遍y[0] 1×4 4y[1] 1×5 2×4 13y[2] 2×5 3×4 22y[3] 3×5 15所以結(jié)果 y [4, 13, 22, 15]長度是 3 2 - 1 4。線性卷積的物理意義很明確它對應(yīng)的是“輸入信號通過一個線性時不變系統(tǒng)的完整響應(yīng)”每個輸出點(diǎn)都是系統(tǒng)對所有歷史輸入的記憶疊加。濾波、系統(tǒng)響應(yīng)、信號通過信道這些場景下我們通常要的都是線性卷積。1.2 循環(huán)卷積周期延拓以后的對位相乘累加循環(huán)卷積的定義看起來也很簡單只是把“反折-移位”里的移位換成了“循環(huán)移位”。所謂循環(huán)移位就是平移之后超出序列范圍的元素不從另一邊補(bǔ)回來而是繞回去像在一個圓環(huán)上轉(zhuǎn)圈。設(shè)兩個序列長度都為 L循環(huán)卷積要求兩個序列按同一個長度 L 來算那么循環(huán)卷積定義為y[k] Σ_{n0}^{L-1} x[n] · h[(k-n) mod L]重點(diǎn)在 (k-n) mod L這就是循環(huán)的意思。算出來的結(jié)果長度也是 L不會變長。所有下標(biāo)都是對 L 取模之后的值所以序列被看成是長度為 L 的周期序列。還是用 x[n] [1, 2, 3]h[n] [4, 5] 來舉例但做循環(huán)卷積時兩個序列都得先擴(kuò)展成同樣長度。假設(shè)取 L 4那就要給 x 補(bǔ)一個零、給 h 補(bǔ)兩個零x [1, 2, 3, 0]h [4, 5, 0, 0]按循環(huán)卷積公式一項(xiàng)項(xiàng)算y[0] x[0]h[0] x[1]h[3] x[2]h[2] x[3]h[1] 1×4 2×5 3×4 0×5 26y[1] x[0]h[1] x[1]h[0] x[2]h[3] x[3]h[2] 1×5 2×4 3×5 0×4 28y[2] x[0]h[2] x[1]h[1] x[2]h[0] x[3]h[3] 1×4 2×5 3×4 0×5 26y[3] x[0]h[3] x[1]h[2] x[2]h[1] x[3]h[0] 1×5 2×4 3×5 0×4 28結(jié)果是 [26, 28, 26, 28]跟線性卷積的 [4, 13, 22, 15] 完全對不上。問題就出在循環(huán)卷積里 h 的下標(biāo)是取模得到的比如 h[(0-1) mod 4] 取到的其實(shí)是 h[3]而在原來的線性卷積里 h[3] 根本不存在它應(yīng)該是補(bǔ)零后的點(diǎn)不是繞回來的值。這就是循環(huán)卷積和線性卷積最本質(zhì)的區(qū)別線性卷積干干凈凈地算完結(jié)果變長循環(huán)卷積把超出邊界的部分“繞回來”疊加到了前面的輸出點(diǎn)上導(dǎo)致結(jié)果被“污染”。1.3 為什么DFT天然對應(yīng)循環(huán)卷積這里要給剛學(xué)到這兒的同學(xué)提一個重要結(jié)論時域循環(huán)卷積對應(yīng)頻域乘積即 DFT(x ? h) X(k) · H(k)。同樣的頻域循環(huán)卷積對應(yīng)時域乘積。但很多人容易忽略一個點(diǎn)DFT處理的是有限長序列可它在數(shù)學(xué)上有一個隱含的“周期延拓”假設(shè)。教科書上會說“DFT是DFS取主值區(qū)間”意思就是你先想象這個序列在無限長周期延拓然后取其中一個周期做變換。既然信號是周期的那么時域卷積自然也是周期的也就是循環(huán)卷積。所以當(dāng)你直接用 FFT 做兩個序列的頻域相乘再逆變換回來得到的結(jié)果本質(zhì)上是循環(huán)卷積不是線性卷積。工程當(dāng)中大多數(shù)場景要的是線性卷積這就直接導(dǎo)致我們在用 FFT 之前必須做一些處理把循環(huán)卷積的結(jié)果“修正”成線性卷積。2. 核心關(guān)系兩者什么時候相等為什么是這個條件2.1 循環(huán)卷積是怎么“混疊”線性卷積結(jié)果的上面那個例子已經(jīng)看到現(xiàn)象了用 L 4 做循環(huán)卷積結(jié)果跟線性卷積完全不是一回事。那如果把 L 逐步變大會發(fā)生什么如果取 L 5那么兩個序列補(bǔ)零成 x [1, 2, 3, 0, 0]h [4, 5, 0, 0, 0]。此時再按循環(huán)卷積公式算y[0] 1×4 4y[1] 1×5 2×4 13y[2] 2×5 3×4 22y[3] 3×5 15y[4] 0結(jié)果是 [4, 13, 22, 15, 0]。你看前面 4 個點(diǎn)跟線性卷積一模一樣就多了一個 0。為什么因?yàn)楫?dāng) L 足夠大之后循環(huán)移位 (k-n) mod L 取到的那些“繞回”位置都正好落在補(bǔ)零的區(qū)域內(nèi)。也就是說原本會從尾部繞到頭部造成混疊的那些項(xiàng)現(xiàn)在全部乘的是 0混疊也就隨之消失了。換個角度想線性卷積結(jié)果長度為 NM-1只要循環(huán)長度 L 不小于這個長度那么線性卷積結(jié)果的 NM-1 個點(diǎn)可以被完整放進(jìn)一個周期里不會出現(xiàn)前后周期相互重疊的情況。反之如果 L NM-1周期延拓后的各個周期就會“疊”在一起重疊部分的數(shù)值加到一起導(dǎo)致輸出跟線性卷積不同。2.2 條件推導(dǎo)L ≥ NM-1 是怎么來的這里可以很直觀地推一下。線性卷積的非零范圍是x 的非零范圍0 ≤ n ≤ N-1h 的非零范圍0 ≤ n ≤ M-1卷積結(jié)果 y[n] 的非零范圍0 ≤ n ≤ NM-2總長度 NM-1。如果用循環(huán)卷積長度 L 來做兩個序列都補(bǔ)零到 L 點(diǎn)。循環(huán)卷積在取模移位時如果它訪問到的 h[(k-n) mod L] 對應(yīng)的是補(bǔ)零區(qū)間就不會產(chǎn)生混疊如果它訪問到的是 h 原本的非零區(qū)間才是正常貢獻(xiàn)。判斷混疊是否發(fā)生關(guān)鍵看線性卷積結(jié)果的“尾巴”會不會在周期延拓時再退回到頭部。線性卷積結(jié)果是一個長度為 NM-1 的有限長序列現(xiàn)在要以 L 為周期延拓。做周期延拓就是把它不斷復(fù)制平移每個周期的位置相隔 L。如果 NM-1 ≤ L那么相鄰兩個周期的結(jié)果剛好首尾相鄰不會重疊取主值區(qū)間后就是完整的線性卷積結(jié)果。如果 NM-1 L那相鄰兩個周期就會有長度為 NM-1-L 的重疊段重疊部分的數(shù)值必須相加這在時域上等價于“尾部折疊疊加到頭部”。所以條件就是循環(huán)卷積長度 L 必須滿足 L ≥ NM-1此時循環(huán)卷積的結(jié)果在 0 到 NM-2 這些點(diǎn)上等于線性卷積其余點(diǎn)為零。這也是為什么用 DFT 做卷積時通常會取 L NM-1 或者取一個不小于它的 2 的冪方便 FFT這個 2 的冪還要大于等于 NM-1。2.3 改變循環(huán)卷積長度的本質(zhì)是改變周期大小這個點(diǎn)可以再深挖一下對真正理解 DFT 卷積特別有幫助。循環(huán)卷積長度 L 一拍腦袋隨便取本質(zhì)上是人為設(shè)定周期延拓的周期。同一個序列周期越短相鄰周期之間的重疊越嚴(yán)重混疊加的項(xiàng)越多。比如同樣是 x[n] [1, 2, 3]h[n] [4, 5]你取 L 3兩個序列各只有 3 個點(diǎn)循環(huán)卷積算出來是y[0] x[0]h[0] x[1]h[2] x[2]h[1] 1×4 2×5 3×4 26這里 h[2] 取的是 h[0] 繞回h[1] 是正常的y[1] x[0]h[1] x[1]h[0] x[2]h[2] 1×5 2×4 3×5 28y[2] x[0]h[2] x[1]h[1] x[2]h[0] 1×4 2×5 3×4 26結(jié)果是 [26, 28, 26]和 L4 時的前三個值一樣因?yàn)檎嬲斐刹町惖木褪俏膊坷@回的那一項(xiàng)。取 L3 時尾部繞回的項(xiàng)更多、疊加得更嚴(yán)重整段結(jié)果都亂了。我自己的理解是循環(huán)卷積的長度 L 決定了周期有多長周期越長線性卷積結(jié)果延拓后互相疊的概率越小L 只要超過結(jié)果總長度就“一點(diǎn)重疊都不?！贝藭r循環(huán)卷積和線性卷積完全一致。所以你不需要記住復(fù)雜的數(shù)學(xué)證明記住“周期延拓重疊則混疊”這八個字就夠了。3. 實(shí)戰(zhàn)用循環(huán)卷積實(shí)現(xiàn)線性卷積的完整流程3.1 核心步驟補(bǔ)零、DFT、相乘、IDFT既然 DFT 做的是循環(huán)卷積那要實(shí)現(xiàn)線性卷積思路就一句話把循環(huán)卷積的周期拉長長到不重疊為止。假設(shè) x[n] 長度 Nh[n] 長度 M目標(biāo)是算 y[n] x[n] * h[n]線性卷積。標(biāo)準(zhǔn)做法如下計算最小循環(huán)長度L_min N M - 1。選實(shí)際 FFT 長度 L為了用 FFT 加速一般取 L 2^ceil(log2(NM-1))即不小于 L_min 的最接近的 2 的冪。當(dāng)然如果直接用 DFTL 取任意不小于 L_min 的數(shù)都行。給 x[n] 后面補(bǔ) L-N 個零給 h[n] 后面補(bǔ) L-M 個零讓兩個序列長度都變成 L。分別做 L 點(diǎn) DFT得到 X[k] 和 H[k]。逐點(diǎn)相乘 Y[k] X[k] * H[k]。對 Y[k] 做 L 點(diǎn) IDFT得到 y[n]。取 y[0] 到 y[NM-2] 這 NM-1 個點(diǎn)作為最終線性卷積結(jié)果后面從 y[NM-1] 開始應(yīng)該是 0數(shù)值上有浮點(diǎn)誤差時會是非常小的數(shù)。第 3 步是這個方法里最“靈魂”的一步。很多新手犯的錯就是“我已經(jīng)補(bǔ)零到兩個序列一樣長了”但實(shí)際上補(bǔ)的零根本不夠。比如 N1000M50有人直接把兩個序列都補(bǔ)零到 1000 點(diǎn)就去做 FFT結(jié)果出來前一段是對的、后一段全亂了就是因?yàn)?L1000 小于 NM-11049尾部混疊已經(jīng)發(fā)生。3.2 一個可以照著跑的代碼示例為了把流程固定下來這里給一段 Python 代碼核心邏輯用 numpy 實(shí)現(xiàn)并把每一步都注釋清楚。這段代碼可以直接復(fù)制到一個 .py 文件里運(yùn)行。import numpy as np # 兩個測試序列 x np.array([1.0, 2.0, 3.0]) h np.array([4.0, 5.0]) N len(x) M len(h) # 1. 最少需要的循環(huán)卷積長度 L_min N M - 1 # 2. 實(shí)際FFT長度為了加速取不小于L_min的2的冪 L 1 while L L_min: L 1 # 3. 補(bǔ)零到長度L x_pad np.zeros(L) x_pad[:N] x h_pad np.zeros(L) h_pad[:M] h # 4. 各自FFT X np.fft.fft(x_pad) H np.fft.fft(h_pad) # 5. 頻域逐點(diǎn)相乘 Y X * H # 6. 逆變換回時域 y_ifft np.fft.ifft(Y).real # 理論上應(yīng)該是實(shí)數(shù)浮點(diǎn)誤差會產(chǎn)生極小的虛部 # 7. 取前NM-1個點(diǎn)作為線性卷積結(jié)果 y_lin y_ifft[:L_min] # 對比用numpy直接算的線性卷積 y_conv np.convolve(x, h) print(直接線性卷積: , y_conv) print(FFT循環(huán)卷積: , y_lin) print(誤差: , np.max(np.abs(y_lin - y_conv)))這段代碼跑出來誤差基本在 1e-14 這個量級可以認(rèn)為是數(shù)值誤差不是算法錯誤。很多人會問為什么 ifft 之后要取 .real因?yàn)楦↑c(diǎn)運(yùn)算里 FFT/IFFT 產(chǎn)生的虛部誤差雖然極小但你如果不取實(shí)部后面看數(shù)據(jù)時總覺得哪里不對勁甚至有人直接把虛部當(dāng)成錯誤信號。實(shí)際上對于實(shí)數(shù)信號IDFT 的理論結(jié)果是實(shí)數(shù)所以取實(shí)部是合理且常規(guī)的操作。3.3 FFT卷積到底快在哪什么時候該用直接線性卷積的復(fù)雜度是 O(NM)因?yàn)槊總€輸出點(diǎn)要做 M 次乘法累加總共 NM-1 個輸出點(diǎn)。FFT 方法復(fù)雜度是 O(L log L)其中 L 取 2 的冪。補(bǔ)零后 L 約等于 NM所以復(fù)雜度近似 O((NM) log(NM))。當(dāng) N 和 M 都大時FFT 方法的優(yōu)勢非常明顯。比如 NM1024直接卷積大約要做 100 萬次乘加而 L2048 的 FFT 大約是 2048×11 ≈ 22528 次蝶形運(yùn)算量差距接近 50 倍。信號越長差距越夸張。但要注意如果其中一個序列特別短比如濾波器長度只有 8 個點(diǎn)輸入也只有 64 個點(diǎn)那這種短序列直接用 conv 反而更快。因?yàn)?FFT 有固定開銷包含補(bǔ)零、正變換兩次、反變換一次、頻域復(fù)數(shù)乘法光這些調(diào)用本身的成本就超過直接卷積了。我實(shí)測過閾值大概在 N、M 都超過 30~50 以后FFT 才開始體現(xiàn)優(yōu)勢。工程上到底用哪種可以先算一下量級再決定別迷信“FFT 一定快”。4. 工程里的延伸長序列濾波和分段卷積4.1 一個輸入無限長濾波器有限長怎么搞實(shí)際工程里經(jīng)常遇到這種情況濾波器 h[n] 長度固定為 M但輸入信號 x[n] 非常長甚至是一個流式信號不能等全部收到再做 FFT。比如實(shí)時音頻處理、雷達(dá)回波處理信號一直在來內(nèi)存也存不下整個序列。這時候就需要分段卷積。分段卷積的基本思路是把長輸入切成長度 L_block 的小塊每塊分別和 h[n] 做線性卷積然后把各塊的卷積結(jié)果按重疊部分相加。這就是重疊相加法overlap-add。做法如下取塊長 L_block要求 L_block M - 1 不超過 FFT 長度通常直接讓 FFT 長度 L L_block M - 1。對每一塊 x_i[n] 補(bǔ)零后做 FFT乘以 H[k]H 是 h 補(bǔ)零后的 FFT只需要算一次再 IFFT。每塊的輸出長度是 L_block M - 1相鄰塊的輸出會有 M-1 個點(diǎn)的重疊。將當(dāng)前塊輸出與上一次保留下來的 M-1 個點(diǎn)相加再輸出有效的前 L_block 個點(diǎn)后 M-1 個點(diǎn)留給下一次疊加。這個過程用代碼實(shí)現(xiàn)并不復(fù)雜但有點(diǎn)繞。核心點(diǎn)就一個每塊獨(dú)立卷積的結(jié)果在重疊區(qū)要相加而不是直接覆蓋。與之對偶的方法是重疊保留法overlap-save它不重疊相加而是重疊保留輸入丟棄輸出中混疊的那部分。兩種方法效果一樣只是實(shí)現(xiàn)細(xì)節(jié)和邊界處理習(xí)慣不同。分段卷積是 FFT 卷積真正發(fā)揮價值的地方因?yàn)槿绻w做 FFT需要等信號全部到齊這在實(shí)時系統(tǒng)里根本不可能。4.2 頻域?yàn)V波其實(shí)也是循環(huán)卷積思想自適應(yīng)濾波、OFDM 調(diào)制解調(diào)、聲學(xué)回聲消除這類系統(tǒng)里塊處理方式非常常見。它們的共同點(diǎn)是先把時域信號分塊轉(zhuǎn)頻域處理再轉(zhuǎn)回時域。如果塊長度和濾波器長度考慮不當(dāng)就很容易出現(xiàn)“看起來處理了但數(shù)據(jù)對不上”的詭異現(xiàn)象。舉個例子OFDM 里為了保護(hù)多徑時延擴(kuò)展會在符號之間插入循環(huán)前綴。循環(huán)前綴的本質(zhì)是一種循環(huán)延拓它讓線性卷積的信道響應(yīng)變成循環(huán)卷積的形式這樣接收端就能直接用頻域均衡做單抽頭補(bǔ)償。這里面“線性卷積變循環(huán)卷積”的思路跟本文講的補(bǔ)零恰好是相反方向的操作但都是同一個思想在不同場景下的應(yīng)用。另一個例子是很多音頻處理教材里提到的“頻域卷積”實(shí)驗(yàn)把一段音頻和某個房間沖激響應(yīng)做卷積模擬混響效果。如果你不補(bǔ)零直接對兩個信號做 FFT 相乘再 IFFT出來的音頻會有一個很明顯的“回環(huán)”雜音這就是循環(huán)卷積的尾部混疊在音頻上聽起來像回聲短路一樣其實(shí)也就是周期延拓導(dǎo)致的時域 aliasing。4.3 圓周卷積譜與快速卷積理論的統(tǒng)一視角如果仔細(xì)琢磨可以感覺到這里面的理論閉環(huán)DFT 是傅里葉變換的離散化周期化版本它處理的所有序列都隱含周期延拓因此在 DFT 域做卷積天然是循環(huán)卷積而線性卷積是物理世界對“有限長輸入經(jīng)過有限長沖激響應(yīng)系統(tǒng)”的描述。要在 DFT 域?qū)崿F(xiàn)物理世界的線性卷積就必須通過補(bǔ)零擴(kuò)大周期消除周期重疊。補(bǔ)零這個操作看似簡單實(shí)質(zhì)是把“循環(huán)卷積”的周期拉大到足以承載“線性卷積”的結(jié)果長度讓兩個數(shù)學(xué)對象在同一個計算平臺上統(tǒng)一起來。這也是很多專業(yè)課程里把 DFT、循環(huán)卷積、快速卷積、重疊相加法串在一起講的原因。理解了這條主線遇到相關(guān)問題時你腦子里會有一個統(tǒng)一的判斷標(biāo)準(zhǔn)這個系統(tǒng)里線性卷積結(jié)果會不會發(fā)生周期重疊重疊了怎么處理是補(bǔ)零防混疊還是故意利用循環(huán)卷積的特性5. 學(xué)習(xí)與實(shí)操中常見的坑5.1 補(bǔ)零長度只補(bǔ)到兩個輸入一樣長這個錯誤簡直太典型了。比如 x 長度 500h 長度 200很多人覺得“我把兩個都補(bǔ)零到 512然后用 512 點(diǎn) FFT”看起來挺對稱。但注意線性卷積長度是 699512 根本不夠。結(jié)果就是卷出來的前 500 個點(diǎn)可能還能看從第 500 個點(diǎn)開始后續(xù)響應(yīng)全部消失尾部混疊值直接疊在頭部附近整體結(jié)果和預(yù)期差一大截。正確做法是先算 NM-1再決定 FFT 長度。哪怕不用 FFT、直接調(diào)庫函數(shù)也要對輸出長度的預(yù)期心里有數(shù)。5.2 忘記取 IFFT 的前 NM-1 個點(diǎn)補(bǔ)零補(bǔ)夠了FFT 也做了IFFT 也做了結(jié)果發(fā)現(xiàn)輸出數(shù)組長度為 L而 L 往往大于預(yù)期的 NM-1。如果直接把整個數(shù)組拿去跟線性卷積結(jié)果比后面多出來的一串近似為 0 的數(shù)會影響誤差判斷有時候甚至有人把多出來的零當(dāng)成有效信號導(dǎo)致后面處理數(shù)組長度不匹配。習(xí)慣做法是IFFT 后只取前 NM-1 個點(diǎn)。多出來的點(diǎn)要么是非常小的浮點(diǎn)噪聲要么是嚴(yán)格為 0。處理信號數(shù)組時長度不匹配的 bug 很難排查最好一開始就指定輸出區(qū)間。5.3 實(shí)數(shù)序列做完 FFT 后埋著復(fù)數(shù)虛部不管很多信號是實(shí)數(shù)FFT 之后變成復(fù)數(shù)IFFT 回來理論上是實(shí)數(shù)但由于浮點(diǎn)誤差虛部會有 1e-15 量級的殘留。如果不取實(shí)部后面繪制波形圖時會看到一些莫名其妙的極小虛部雖然不影響幅值但有些工具里會警告“數(shù)據(jù)是復(fù)數(shù)”甚至導(dǎo)致畫圖報錯。建議在 IFFT 后用 .real 取實(shí)部。但也別用全盤取絕對值的方式處理因?yàn)樾盘柪锟赡苡胸?fù)值取絕對值會把負(fù)半軸整體翻上去破壞波形。5.4 下標(biāo)從 1 還是從 0最容易把人繞暈MATLAB 里索引從 1 開始Python 里從 0 開始。線性卷積定義里 y[n] 的 n 也是從 0 開始標(biāo)。寫代碼的時候補(bǔ)零位置、取主值區(qū)間、分段卷積的塊索引這些地方但凡少一個 1 或者多一個 -1結(jié)果就可能整體偏移一位。而且這種偏移在波形圖上肉眼幾乎看不出來只有和參考結(jié)果做數(shù)值對比時才能發(fā)現(xiàn)。我的經(jīng)驗(yàn)是寫卷積相關(guān)代碼時不要直接看公式抄先把一個特別小的測試序列跑起來比如 x[1,2]h[3,4]手算出線性卷積結(jié)果然后程序跑一遍逐步打印中間過程來核對索引。5.5 用頻域相乘實(shí)現(xiàn)卷積時忘了 H 也要補(bǔ)零有時候 x 補(bǔ)零了H 的 FFT 卻用的是原始長度或者反過來。這會造成 X[k] 和 H[k] 的頻譜點(diǎn)數(shù)不一致FFT 長度不匹配時直接相乘會報錯或者更隱蔽的是長度碰巧一致但補(bǔ)零規(guī)則不一樣比如 X 補(bǔ)到了 1024H 只補(bǔ)到了 512然后廣播相乘結(jié)果完全沒有意義。所有參與 DFT 的序列必須補(bǔ)到同一個長度 L。這里沒有例外。5.6 把 FFT 卷積用在小規(guī)模序列上反而更慢這一點(diǎn)是性能上的坑。FFT 卷積雖然在大規(guī)模下有巨大優(yōu)勢但小規(guī)模下并不劃算。比如兩個長度都只有 8 的序列直接用雙重循環(huán) 64 次乘法就夠了FFT 反而要做補(bǔ)零、64 點(diǎn) FFT、復(fù)數(shù)乘法、IFFT耗時可能是直接卷積的 5 到 10 倍。我在實(shí)際項(xiàng)目里一般會設(shè)置一個閾值比如 N*M 2000 時才切換 FFT 方案低于閾值直接暴力卷積整體性能最優(yōu)。5.7 驗(yàn)證結(jié)果時只看“像不像”不看數(shù)值做實(shí)驗(yàn)最忌諱的就是輸出畫出來“看起來差不多”就交差。循環(huán)卷積與線性卷積的差異不是噪聲而是結(jié)構(gòu)性的。如果你補(bǔ)零長度差得不多混疊可能只出現(xiàn)在特定區(qū)間圖像整體趨勢接近但個別點(diǎn)偏差很大。這種問題肉眼不一定看得出來但數(shù)值一對比就暴露了。我自己的習(xí)慣是每次跑通 FFT 卷積后用 np.convolve 或者 MATLAB 的 conv 函數(shù)做一次對照計算 max(abs(y_fft - y_conv))只要這個值在 1e-10 以下說明實(shí)現(xiàn)沒問題。之后再去做長信號的工程實(shí)現(xiàn)心里就有底了。最后說點(diǎn)實(shí)際的循環(huán)卷積和線性卷積的關(guān)系說到底是數(shù)學(xué)在“有限長”這件事上的一次精妙折中DFT 只能處理有限長序列但它骨子里又是周期性的所以卷積必須按循環(huán)來做而物理世界大量場景需要的線性卷積恰好可以通過把周期拉長來“假裝”實(shí)現(xiàn)。補(bǔ)零不是為了好看也不是為了湊 2 的冪而是為了給線性卷積結(jié)果騰出一個不重疊的周期空間。我在實(shí)驗(yàn)室里第一次自己寫 FFT 卷積濾波器的時候就是沒把補(bǔ)零長度當(dāng)回事用的 512 點(diǎn) FFT 去處理 300 點(diǎn)信號和 200 點(diǎn)沖激響應(yīng)結(jié)果輸出從第 113 個點(diǎn)左右開始完全變形。后來畫圖才發(fā)現(xiàn)是尾部周期混疊一下就想通了“L ≥ NM-1”這個條件到底在防什么東西。從那以后再遇到頻域卷積、分段濾波我第一件事就是算清楚長度再開始寫代碼。希望這篇文章也能讓你少走這個彎路。