长序列的快速卷积

一、功能

用重叠保留法和快速傅里叶变换计算一个长序列和一个短序列的快速卷积。它通常用于数字滤波。

二、方法简介

设序列\(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的基础上添加系数长度的倒数。