特别长序列的快速卷积

一、功能

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

二、方法简介

设序列\(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}, \ n=0,1,...,N-1 \]

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)] \]

由于数据量很大，因此序列\(h(n)\)和特别长序列\(x(n)\)要存储到磁盘中，计算完毕后，卷积序列\(y(n)\)也由程序自动存储到磁盘中。

三、使用说明

快速卷积的C语言实现方式如下

/************************************
	hfname		----字符指针。它指向短序列h(i)的文件名。
	xfname		----字符指针。它指向长序列x(i)的文件名。
	yfname		----字符指针。它指向卷积y(i)的文件名。
	n			----整型变量。对长序列x(i)进行分段的长度。一般选取n大于短序列h(i)长度m的两倍以上，
				且必须是2的整数次幂，即n=2^gamma。	
************************************/
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "rfft.c"
#include "irfft.c"
void convold(char *hfname, char *xfname, char *yfname, int n)
{
	int i, j, m, i1, i2, ls0, len, num, nblks;
	double *h, *r, *s;
	FILE *fp, *fp1, *fp2;
	r = malloc(n * sizeof(double));
	h = malloc(n * sizeof(double));
	n2 = n / 2;
	fopen(hfname, "r");
	i = 0;
	while(!fopen(fp)){
		fscanf(fp, "%lf", &h[i++]);
	}
	fclose(fp);
	m = i - 1;
	for(i = m; i < n; i++){
		h[i] = 0.0;
	}	
	rfft(h, n);
	s = malloc((n - m + 1) * sizeof(double));
	num = n - m + 1;
	fp1 = fopen(xfname, "r");
	fp2 = fopen(yfname, "r");
	len = 5000;
	ls0 = 0;
	do{
		for(i = ls0; i < 5000; i++){
			if(feof(fp1)){
				len = i;
				break;
			}
			fscanf(fp1, "%lf", &x[i]);
		}
		nblks = floor((len - n + m) / (double)num) + 1;
		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 = 0; i < i1; i++)
			fprintf(fp2, "%lf\n", x[i]);

		for(i = i1; i < len; i++){
			i2 = i - i1;
			x[i2] = x[i];
		}
		ls0 = len - i1;
	}while(!feof(fp1))
	fclose(fp);
	fclose(fp1);
	fclose(fp2);
	free(r);
	free(h);
	free(s);
}

其中rfft.c文件请参考实序列快速傅里叶变换（一）

irfft.c在rfft.c的基础上添加系数长度的倒数。