傅里叶和互相关算法及c程序

当x（n）=y（n）时，得到x（n）的自相关函数为：

上面的推导表明，互相关和自相关函数的计算可利用FFT实现。由于离散付里叶变换隐含着周期性，所以用FFT计算离散相关函数也是对周期序列而言的。直接做N点FFT相当于对两个N点序列x(n)、y(n)作周期延拓，作相关后再取主值（类似圆周卷积）。而实际一般要求的是两个有限长序列的线性相关，为避免混淆，需采用与圆周卷积求线性卷积相类似的方法，先将序列延长补0后再用上述方法。

利用FFT求两个有限长序列线性相关的步骤： 设x（n）长为N1，y（n）长为N2，求线性相关。

（1）为了使两个有限长序列的线性相关可用其圆周相关代替而不产生混淆，选择周期N≥N1+N2-1，，且N=2m，以便使用FFT，将x（n），y（n）补零至长为N。即：

（2）用FFT计算X(k)，Y(k)（k=0，1…，N-1）

（3）R（k）=X（k）Y（k）

（4）对R（k）作IFFT，得到r（n）（n=0，1，…，N-1）

*

/**************fft*****************/

//先创建两个txt格式的文件，分别为x.txt和y.txt，用来存放互相关变换的两组数据；运行完产生result.txt 和result2.txt和result.txt，用来保存结果。 #include #include #include #include #include using namespace std;

#define PI 3.1415926535897932384626433832795028841971 #define N 128

typedef struct {

double real;/*实部*/ double img;/*虚部*/ }complex;

complex x[N], y[N],R[N], *W;/*输出序列的值*/ int size_x=0,size_y=0,length=0;//序列长度

void add(complex a,complex b,complex *c) {

c->real=a.real+b.real; c->img=a.img+b.img; }

void sub(complex a,complex b,complex *c) {

c->real=a.real-b.real; c->img=a.img-b.img; }

void mul(complex a,complex b,complex *c) {

c->real=a.real*b.real - a.img*b.img; c->img=a.real*b.img + a.img*b.real; }

void divi(complex a,complex b,complex *c) {

c->real=( a.real*b.real+a.img*b.img )/( b.real*b.real+b.img*b.img); c->img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);

}

void initW(int size_x)//核运算 {

int i;

W=(complex *)malloc(sizeof(complex) * size_x); for(i=0;i

W[i].real=cos(2*PI/size_x*i); W[i].img=-1*sin(2*PI/size_x*i); } }

void change()//变址运算 {

complex temp;

unsigned short i=0,j=0,k=0; double t;

for(i=0;i

k=i;j=0;

t=(log(size_x)/log(2)); while( (t--)>0 ) {

j=j<<1;

j|=(k & 1); k=k>>1; } if(j>i) {

temp=x[i]; x[i]=x[j]; x[j]=temp; } } }

void fftx()//快速傅里叶函数 {

int i=0,j=0,k=0,l=0;

complex up,down,product; change();

for(i=0;i

l=1<

for(j=0;j

mul(x[j+k+l],W[size_x*k/2/l],&product); add(x[j+k],product,&up); sub(x[j+k],product,&down); x[j+k]=up;

x[j+k+l]=down; } } } }

void ffty()//快速傅里叶函数 {

int i=0,j=0,k=0,l=0;

complex up,down,product; change();

for(i=0;i

for(j=0;j

mul(y[j+k+l],W[size_x*k/2/l],&product); add(y[j+k],product,&up); sub(y[j+k],product,&down); y[j+k]=up;

y[j+k+l]=down; } } } }

void ifft() //傅里叶逆变换 { length=size_x+size_y-1; int i=0,j=0,k=0,l=length; complex up,down;

for(i=0;i<(int)(log(length)/log(2));i++) /*一级蝶形运算*/