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