重疊-儲存之摺積法

Template:NoteTA 重疊-儲存之摺積法 ( Overlap-save method, Overlap-discard method ) 是一種區塊摺積 ( block convolution, sectioned convolution )，可以有效的計算一個很長的信號 x[n] 和一個 FIR 濾波器 h[n] 的離散摺積。

$\begin{array}{c}y[n]=x[n]\star h[n]\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}\sum _{m=-\mathrm{\infty }}^{\mathrm{\infty }}h[m]\cdot x[n-m]=\sum _{m=1}^{M}h[m]\cdot x[n-m]\end{array}$

其中 h[m] 在 [1, M] 之外為零。

與重疊-相加之摺積法不同之處在於，重疊-儲存之摺積法所算出的輸出區塊並不重疊 (因此計算上少了將輸出區塊相加所需的加法)，而是每次用的輸入區塊有所重疊。因此實作時每次讀取輸入後需將和下一個輸入重疊的部分儲存起來，作為下一輸入區塊的開頭部份，因此稱為重疊-儲存之摺積法。另外重疊-儲存之摺積法也不需補零。

算法

概念上，這個做法是選用一個較短的適當長度 L 來切割 y[n] ，則因為 h[n] 是有限長度，因此在某一區塊內的 y[n] 也只被有限長的 x[n] 區塊（會比 y[n] 分割成的區塊長一點）所決定。因此只要選擇有影響的輸入區塊和 h[n] 摺積，再選擇結果中適當的部分即可得到正確的輸出區塊。

$x_k[n]\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}\{\begin{array}{cc}x[n+kL-M]& 1\le n\le L+M-1\\ 0& \text{otherwise}.\end{array}$

$y_k[n]\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}{x}_k[n]\star h[n]$

則對於在 $[kL+1,(k+1)L]$ 內的 n ， 輸出 y[n] 可寫成

$\begin{array}{cc}y[n]=\sum _{m=1}^{M}h[m]\cdot x[n-m]=\sum _{m=1}^{M}h[m]\cdot x_k[n+M-kL-m]& =x_k[n+M-kL]\star h[n]\\ & \stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}{y}_k[n+M-kL].\end{array}$

所以只需計算 n 在 $[kL+1,(k+1)L]$ 中的 y_k[n + M - kL] ，亦即 n 在 $[M+1,L+M]$ 的 y_k[n] 部份即可。因此每一段輸出區塊 y_k[n] 的前 M-1 點可丟棄（discard）。

儘管一時看不出切割成區塊的好處為何，但將 x_k[n] 做  $N\ge L+M-1$  的週期延伸，

$x_{k,N}[n]\stackrel{\mathrm{d}\mathrm{e}\mathrm{f}}{=}\sum _{k=-\mathrm{\infty }}^{\mathrm{\infty }}{x}_k[n-kN],$

則  $(x_{k,N})\star h$  和  $x_k\star h$  這兩個摺積在 $[M+1,L+M]$ 的部份相等。所以可以將線性摺積改以 $N$ 點圓周摺積計算，結果的 $[M+1,L+M]$ 部分作為輸出 y[n] 在 $[kL+1,(k+1)L]$ 的部份。由於每段 x_k[n] 原本就有  $L+M-1$  長，所以選擇  $N=L+M-1$  的話輸入 x[n] 就不需補零。 改以圓周摺積計算後即可藉圓周摺積定理

$y_k[n]=\text{IFFT}(\text{FFT}(x_k[n])\cdot \text{FFT}(h[n]))$

轉換成三次 $N$ 點快速傅立葉變換和 $N$ 次乘法，使原本每段 O(N²) 的運算量減少至 O(N logN)，速度大幅增加。

準程式碼

   (Overlap-save algorithm for linear convolution)
   //////// revised by fantastic ////////
   N = length(x), M = length(h)
   O = M – 1;	// overlap length must be M-1
   L = M;	// >=1 is OK
   P = O + L;
   H = FFT(h, P);	// just calc once
   idx = - (O - 1);	// starting index which is offset M-1 in matlab
   while (idx <= N)
       i1 = max(1, idx);	// must be >= 1
       i2 = min(N, idx+P-1);	// must be <= N
       yt = IFFT( FFT(x(i1:i2), P).*H, P );
       y(idx:idx+P-M) = yt(M:P);	// discard first M-1 values and concatenate the remaining
       idx = idx + L;
   end
   y = y(1:M+N-1);	// the first M+N-1 values are the convolution result

區塊長度的選擇

當 x[n] 的長度 N' 和 h[n] 的長度 M 相差太大時（例如 M < log₂N' ），直接摺積（不透過圓周摺積和 FFT ）反而最快。而當 N' 和 M 差不多在同一個數量級時，不用分割，也就是只有一塊長度 L = N' 的區塊去做 FFT 即可。而當 N' 比 M 大了不少，卻沒大太多時，區塊長度 L 就需要選擇。除了與 N' 和 M 相關以外，也要考慮當兩者相除有餘數時，剩下一小段的輸入可能會造成浪費。

相關條目

• 離散摺積
• 圓周摺積
• 快速傅立葉變換
• 重疊-相加之摺積法

參考文獻

• Template:Cite book
• Template:Citation

外部連結

• DSP class Fall 2005 Lecture08Template:Dead link at The University of Texas at Arlington]