重叠-存储之卷积法

2018-09-06 14:26:19 所属分类：数值分析

重叠-存储之卷积法 ( Overlap-save method, Overlap-discard method ) 是一种区块卷积 ( block convolution, sectioned convolution )，可以有效的计算一个很长的信号 xn 和一个 FIR 滤波器 hn 的离散卷积。

{displaystyle {begin{aligned}yn=xnstar hn {stackrel {mathrm {def} }{=}} sum _{m=-infty }^{infty }hmcdot xn-m=sum _{m=1}^{M}hmcdot xn-mend{aligned}}}

其中 hm 在 1, M 之外为零。

与重叠-相加之卷积法不同之处在于，重叠-存储之卷积法所算出的输出区块并不重叠 (因此计算上少了将输出区块相加所需的加法)，而是每次用的输入区块有所重叠。因此实现时每次读取输入后需将和下一个输入重叠的部分存储起来，作为下一输入区块的开头部分，因此称为重叠-存储之卷积法。另外重叠-存储之卷积法也不需补零。

算法

概念上，这个做法是选用一个较短的适当长度 L 来切割 yn ，则因为 hn 是有限长度，因此在某一区块内的 yn 也只被有限长的 xn 区块（会比 yn 分区成的区块长一点）所决定。因此只要选择有影响的输入区块和 hn 卷积，再选择结果中适当的部分即可得到正确的输出区块。

{displaystyle x_{k}n {stackrel {mathrm {def} }{=}}{begin{cases}xn+kL-M&1leq nleq L+M-1\0&{textrm {otherwise}}.end{cases}}}

{displaystyle y_{k}n {stackrel {mathrm {def} }{=}} x_{k}nstar hn,}

则对于在 ${displaystyle kL+1,(k+1)L}$ 内的 n ，输出 yn 可写成

{displaystyle {begin{aligned}yn=sum _{m=1}^{M}hmcdot xn-m=sum _{m=1}^{M}hmcdot x_{k}n+M-kL-m&=x_{k}n+M-kLstar hn\&{stackrel {mathrm {def} }{=}} y_{k}n+M-kL.end{aligned}}}

所以只需计算 n 在 ${displaystyle kL+1,(k+1)L}$ 中的 y_kn + M - kL ，亦即 n 在 ${displaystyle M+1,L+M}$ 的 y_kn 部分即可。因此每一段输出区块 y_kn 的前 M-1 点可丢弃（discard）。

尽管一时看不出切割成区块的好处为何，但将 x_kn 做 ${displaystyle Ngeq L+M-1,}$ 的周期延伸，

{displaystyle x_{k,N}n {stackrel {mathrm {def} }{=}} sum _{k=-infty }^{infty }x_{k}n-kN,}

则 ${displaystyle (x_{k,N})star h,}$ 和 ${displaystyle x_{k}star h,}$ 这两个卷积在 ${displaystyle M+1,L+M}$ 的部分相等。所以可以将线性卷积改以 ${displaystyle N,}$ 点循环卷积计算，结果的 ${displaystyle M+1,L+M}$ 部分作为输出 yn 在 ${displaystyle kL+1,(k+1)L}$ 的部分。由于每段 x_kn 原本就有 ${displaystyle L+M-1,}$ 长，所以选择 ${displaystyle N=L+M-1,}$ 的话输入 xn 就不需补零。改以循环卷积计算后即可藉循环卷积定理

{displaystyle y_{k}n={textrm {IFFT}}left({textrm {FFT}}left(x_{k}nright)cdot {textrm {FFT}}left(hnright)right)}

变换成三次 ${displaystyle N,}$ 点快速傅里叶变换和 ${displaystyle 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

区块长度的选择

当 xn 的长度 N' 和 hn 的长度 M 相差太大时（例如 M < log₂N' ），直接卷积（不透过循环卷积和 FFT ）反而最快。而当 N' 和 M 差不多在同一个数量级时，不用分区，也就是只有一块长度 L = N' 的区块去做 FFT 即可。而当 N' 比 M 大了不少，却没大太多时，区块长度 L 就需要选择。除了与 N' 和 M 相关以外，也要考虑当两者相除有余数时，剩下一小段的输入可能会造成浪费。

参考文献

Rabiner, Lawrence R.; Gold, Bernard. Theory and application of digital signal processing. Englewood Cliffs, N.J.: Prentice-Hall. 1975: pp 65–67. ISBN 0-13-914101-4.
Helms, H., Fast Fourier transform method of computing difference equations and simulating filters, IEEE Transactions on Audio and Electroacoustics, 1967, 15(2): 85–90