fft.cpp

实数fft代码

#include "stdio.h"
#include "iostream.h"
#include "fstream.h"
#include "fft_real.c"
#include "Math.h"
#define change(x,y) {float tmp=(x);(x)=(y);(y)=tmp;}



void main()
{

	FILE *cfPtr;
	short array[7000],xx[513];
	
	if( ( cfPtr = fopen ("subarraydata.dat","rb")) == NULL)
		printf("Can't open data file");
	else{
	
			fread(array, sizeof(short), 6000, cfPtr);
	}
	int N =512;
	for (int i=0;i<N;i++)
		xx[i]=array[i*7];
	FFT((float *)(xx),N);

}
/*
void FFT(float xxx[], int N)
{
	float temp,tempr,tempi; 
	float wtemp,wr,wpr,wpi,wi,theta;
	float *data;
	float  x[2050];

	int n,mmax,m,q,istep,p,i; 
	float a,c,e,s,fr,fi,gr,gi;
	data=&xxx[-1];
	n=N<<1; q=1; 
	for(p=1;p<n;p+=2) //倒码序
	{
		if(q>p)
		{ 
			change(data[q],data[p]); 
			change(data[q+1],data[p+1]); 
		} 
		m=n>>1;
		while(m>=2&&q>m)
		{
			q-=m;m>>=1;
		} 
		q+=m; 
	} 
 
	mmax=2;
	while(n>mmax)
	{ 
		istep=2*mmax; 
		theta=-6.28318530717959/mmax; 
		wtemp=sin((double)(0.5*theta)); 
		wpr=-2.0*wtemp*wtemp; 
		wpi=sin((double)theta);
		
		for(wr=1.0,wi=0.0,m=1;m<mmax;m+=2)
		{ 
			for(p=m;p<=n;p+=istep)
			{ 
				q=p+mmax; 
				tempr=wr*data[q]-wi*data[q+1]; 
				tempi=wr*data[q+1]+wi*data[q]; 
				data[q]=data[p]-tempr; 
				data[q+1]=data[p+1]-tempi; 
				data[p]+=tempr; 
				data[p+1]+=tempi;
				} 
			wr=(wtemp=wr)*wpr-wi*wpi+wr; 
			wi=wi*wpr+wtemp*wpi+wi;  
		} 
		mmax=istep;
	}

	//得到xxx[0]、[2]……[2*N-2]为实部，[1]、[3]……[2*N-1]为虚部
	e=3.14159265359/N;
	x[0]=xxx[0]+xxx[1];
	x[N]=xxx[0]-xxx[1];
	for(i=1;i<N;i++) 
	{
		fr=(xxx[2*i]+xxx[n-2*i])/2;
		fi=(xxx[2*i+1]-xxx[n-2*i+1])/2;
		gi=(xxx[2*i+1]+xxx[n-2*i+1])/2;
		gr=(xxx[2*i]-xxx[n-2*i])/2;
		a=i*e;
		c=cos(a);
		s=sin(a);
		x[i]=fr+c*gi-s*gr;
		x[n-i]=fi-c*gr-s*gi;
	}
	for(i=0;i<n;i++)
		xxx[i]=x[i];
;
}
*/