一、功能
用重疊保留法和快速傅里葉變換計算一個長序列和一個短序列的快速卷積。它通常用於數字濾波。
二、方法簡介
設序列\(x(n)\)的長度為\(L\)，序列\(h(n)\)的長度為\(M\)，序列\(x(n)\)與\(h(n)\)的線性卷積定義為
\[y(n)=\sum_{i=0}^{M-1}x(i)h(n-i) \]
用重疊保留法和快速傅里葉變換計算線性卷積的算法如下：
1、將序列\(h(n)\)按如下方式補零，形成長度為\(N=2^{\gamma }\)的序列
\[\begin{matrix}h(n)=\left\{\begin{matrix}\begin{align*}h(n) &, n=0,1,...,M-1 \\ 0 &, n=M,M+1,...,N-1\end{align*}\end{matrix}\right.\\ \end{matrix} \]
2、用FFT算法計算序列\(h(n)\)的離散傅里葉變換\(H(k)\)
\[H(k)=\sum_{n=0}^{N-1}h(n)e^{-j2\pi nk/N} \]
3、將長序列\(x(n)\)分為長度為\(N\)的小段\(x_i(n)\)，相鄰段間重疊\(M-1\)點，即
\[\begin{matrix}x_i(n)=\left\{\begin{matrix}\begin{align*}x[n+i(N-M+1)-(M-1)] &, 0 \leqslant n \leqslant N-1 \\ 0 &, others\end{align*}\end{matrix}\right.\\ i=0,1,...\end{matrix} \]
注意：對於第一段的\(x_0(n)\)，由於沒有前一段的保留信號，因此要在其前面填補\(M-1\)個零。
4、用FFT算法計算序列\(x_i(n)\)的離散傅里葉變換\(X_i(k)\)
\[X_i(k)=\sum_{n=0}^{N-1}x_i(n)e^{-j2\pi nk/N} \]
5、計算\(X_i(k)\)與\(H(k)\)的乘積
\[Y_i(k)=X_i(k)H(K) \]
6、用FFT算法計算\(Y_i(k)\)的離散傅里葉反變換，得到序列\(x_i(n)\)與\(h(n)\)的卷積\(y_i(n)\)
\[y_i(n)=\frac{1}{N}\sum_{k=0}^{N-1}Y_i(k)e^{j2\pi nk/N} \]
7、將\(y_i(n)\)的前\(M-1\)點舍去，僅保留后面的\(N-M+1\)個樣本。
\[\begin{matrix}\bar{y}_i(n)=\left\{\begin{matrix}\begin{align*}y(n) &, M-1 \leqslant n \leqslant N-1 \\ 0 &, others\end{align*}\end{matrix}\right.\\ \end{matrix} \]
8、重復步驟3~7，直到所有分段算完為止。
9、考慮到\(x(n)\)分段時，\(x_i(n)\)有\(M-1\)點與前一段重疊，新添的樣本只有\(N-M+1\)個，所以\(y(n)\)由\(\bar{y}_i(n)\)首尾銜接而成，即
\[y(n)=\sum_{i=0}^{\infty}\bar{y}_i[n-i(N-M+1)+(M-1)] \]
三、使用說明
快速卷積的C語言實現方式如下
/************************************
x ----雙精度一維數組，長度為len。開始時存放實序列x(i)，最后存放線性卷積的值。
h ----雙精度一維數組，長度為m。開始時存放實序列h(i)。
len ----整型變量。長序列x(i)的長度。
m ----整型變量。短序列h(i)的長度。
n ----整型變量。對長序列x(i)進行分段的長度。一般選取n大於短序列h(i)長度m的兩倍以上，
且必須是2的整數次冪，即n=2^gamma。
************************************/
#include "math.h"
#include "stdlib.h"
#include "rfft.c"
#include "irfft.c"
void convols(double *x, double *h, int len, int m, int n)
{
int i, j, i1, n2, num, nblks;
double t, *r, *s;
r = malloc(n * sizeof(double));
s = malloc((n - m + 1) * sizeof(double));
n2 = n / 2;
num = n - m + 1;
nblks = floor((len - n + m) / (double)num) + 1;
for(i = m; i < n; i++){
h[i] = 0.0;
}
rfft(h, n);
for(j = 0; j < nblks; j++){
if(j == 0){
for(i = 0; i < (m - 1); i++)
r[i] = 0.0;
for(i = m - 1; i < n; i++){
i1 = i - m + 1;
r[i] = x[i1];
}
} else {
for(i = 0; i < n; i++){
i1 = i + j * num - m + 1;
r[i] = x[i1];
}
for(i = 0; i < num; i++){
i1 = i + (j - 1) * num;
x[i1] = s[i];
}
}
rfft(r, n);
r[0] = r[0] * h[0];
r[n2] = r[n2] * h[n2];
for(i = 1; i < n2; i++){
t = h[i] * r[i] - h[n -i] * r[n - i];
r[n - i] = h[i] * r[n - i] + h[n - i] * r[i];
r[i] = t;
}
irfft(r, n);
for(i = (m - 1); i < m; i++){
i1 = i - m + 1;
s[i1] = r[i];
}
}
for(i = 0; i < num; i++){
i1 = i + (j - 1) * num;
x[i1] = s[i];
}
i1 = j * num;
for(i = i1; i < len; i++)
x[i] = 0.0;
free(r);
free(s);
}
其中rfft.c文件請參考實序列快速傅里葉變換（一）
irfft.c在rfft.c的基礎上添加系數長度的倒數。