重疊-儲存之摺積法

维基百科,自由的百科全书
跳转至: 导航搜索

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


\begin{align}
y[n] = x[n] \star h[n] \ \stackrel{\mathrm{def}}{=} \ \sum_{m=-\infty}^{\infty} h[m] \cdot x[n-m]
= \sum_{m=1}^{M} h[m] \cdot x[n-m]
\end{align}

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

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

算法[编辑]

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


x_k[n]  \ \stackrel{\mathrm{def}}{=}
\begin{cases}
x[n+kL-M] & 1 \le n \le L+M-1\\
0 & \textrm{otherwise}.
\end{cases}
y_k[n] \ \stackrel{\mathrm{def}}{=} \ x_k[n] \star h[n]\,


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


\begin{align}
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{def}}{=} \ y_k[n+M-kL].
\end{align}

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


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

x_{k,N}[n] \ \stackrel{\mathrm{def}}{=} \ \sum_{k=-\infty}^{\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] 的部份。由於每段 xk[n] 原本就有  L+M-1\,  長,所以選擇  N=L+M-1\,  的話輸入 x[n] 就不需補零。 改以圓周摺積計算後即可藉圓周摺積定理

y_k[n] = \textrm{IFFT}\left(\textrm{FFT}\left(x_k[n]\right)\cdot\textrm{FFT}\left(h[n]\right)\right)

轉換成三次 N\,快速傅立葉變換N\, 次乘法,使原本每段 O(N2) 的運算量減少至 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 < log2N' ),直接摺積(不透過圓周摺積FFT )反而最快。而當 N'M 差不多在同一個數量級時,不用分割,也就是只有一塊長度 L = N' 的區塊去做 FFT 即可。而當 N'M 大了不少,卻沒大太多時,區塊長度 L 就需要選擇。除了與 N'M 相關以外,也要考慮當兩者相除有餘數時,剩下一小段的輸入可能會造成浪費。

相關條目[编辑]

參考文獻[编辑]

外部連結[编辑]

  • DSP class Fall 2005 Lecture08 at The University of Texas at Arlington]