一、功能
用快速傅里叶变换计算两个有限长序列的线性相关。
二、方法简介
设序列(x(n))的长度为(M)，序列(y(n))的长度为(N)，序列(x(n))与(y(n))的互相关定义为
[z(n)=sum_{i=0}^{M-1}x(i)y(n+i), n=-(M-1),...,N-1
]
序列(x(n))的自相关定义为
[z(n)=sum_{i=0}^{M-1}x(i)x(n+i), n=-(M-1),...,N-1
]
显然，如果(M=N)并且序列(x(n))与(y(n))相同，那么互相关就变成了自相关。
用快速傅里叶变换计算线性相关的算法如下：
1、选择(L)满足下述条件
[left{egin{matrix}egin{align*}L &geqslant M + N - 1\ L &= 2^{gamma }, gamma is positive integerend{align*}end{matrix}
ight.
]
2、将序列(x(n))与(y(n))按如下方式补零，形成长度为(L=2^{gamma })的序列
[ar{x}(n)=left{egin{matrix}egin{align*}x(n) & , n=0,1,...,M-1\ 0 & , n=0,1,...,M,M+1,...,L-1end{align*}end{matrix}
ight.
]
[ar{y}(n)=left{egin{matrix}egin{align*}0 & , n = 0,1,...,M-2\ y(n) & , n = M-1, M,...,M+N-2\ 0 & , n = M+N-1, M+N,...,L-1end{align*}end{matrix}
ight.
]
3、用FFT算法分别计算序列(ar{x}(n))和(ar{y}(n))的离散傅里叶变换(X(k))与(Y(k))
[X(k)=sum_{n=0}^{L-1}ar{x}(n)e^{-j2pi nk/L}
]
[Y(k)=sum_{n=0}^{L-1}ar{y}(n)e^{-j2pi nk/L}
]
4、计算(X(k))与(Y(k))的乘积
[Z(k)=X^*(k)Y(K)
]
其中(*)表示复共轭；
5、用FFT算法计算(Z(k))的离散傅里叶反变换，得到线性相关(z(n))
[z(n)=frac{1}{N}sum_{k=0}^{L-1}Z(k)e^{j2pi nk/L}, n=0,1,...,L-1
]
序列(z(n))的前(M+N-1)点的值就是序列(x(n))与(y(n))的线性相关。
三、使用说明
快速卷积的C语言实现方式如下
/************************************
x ----双精度实型一维数组，长度为m。开始时存放实序列x(i)，最后存放线性相关的值。
y ----双精度实型一维数组，长度为n。开始时存放实序列y(i)。
m ----整型变量。序列x(i)的长度。
n ----整型变量。序列y(i)的长度。
len ----整型变量。len>=m+n-1，且必须是2的整数次幂，即len=2^gamma
************************************/
#include "stdlib.h"
#include "rfft.c"
#include "irfft.c"
void correl(double *x, double *y, int m, int n, int len)
{
int i, len2;
double t, *z;
z = mlloc(len * sizeof(double));
for(i = m; i < len; i++)
x[i] = 0.0;
for(i = 0; i < (m - 1); i++)
z[i] = 0.0;
for(i = (m - 1); i <= (m + n - 2); i++)
z[i] = y[i - m + 1];
for(i = (m + n - 1); i < len; i++)
z[i] = 0.0;
rfft(x, len);
rfft(z, len);
len2 = len / 2;
x[0] = x[0] * z[0];
x[len2] = x[len2] * z[len2];
for(i = 1; i < len2; i++){
t = x[i] * z[i] + x[len - i] * z[len - i];
x[len - 1] = x[i] * z[len - i] - x[len - i] * z[i];
x[i] = t;
}
irfft(x, len);
free(z);
}
其中rfft.c文件请参考实序列快速傅里叶变换（一）
irfft.c在rfft.c的基础上添加系数长度的倒数。