fft.cpp

matlab版图像数字处理这本书的配套代码

#include "stdafx.h"
#include "FFT.h"
#include<math.h>
#include <complex>
using namespace std;

void  FFT(complex<double>* x,complex<double>* y,int M)
{	
	int N=1<<M;		//N=2^M
	complex<double>* w;
	w=new complex<double>[N/2];
	double PI2=8.0*atan(1.0);
	
	int i,j,p,k,L,J,P,B;
	//ji suan w
	for(i=0;i<N/2;i++)
		w[i]=complex<double>(cos(i*PI2/N),sin(i*PI2/N));
	
	//ma wei dao xu
	///////////////////////////////////////////////////////////
	//ChangeOrder(x,M);
	complex<double>* pBak;
	pBak=new complex<double>[N];
	memcpy(pBak,x,sizeof(complex<double>)*N);
	
	for(j=0;j<N;j++)//change order for N points
	{
		p=0;
		for(i=0;i<M;i++)
		{
			if(j&(1<<i))
				p+=1<<(M-i-1);
		}
		x[j]=pBak[p];
	}
	delete []pBak;
	///////////////////////////////////////////////////////////
	//die xing yun suan 
	for(L=1;L<=M;L++)// L level
	{	
		B=1<<(L-1);		//B=bfsize/2
		for(J=0;J<=B-1;J++)
		{
			P=J*(1<<(M-L));//w de zhishu
			complex<double> temp;
			for(k=J;k<=N-1;k+=(1<<L))
			{
				temp=x[k+B]*w[P];
				x[k+B]=x[k]-temp;
				x[k]=x[k]+temp;
			}
		}
	}
	memcpy(y,x,sizeof(complex<double>)*N);
	delete []w;
}