fft.txt

用C++编写的FFT函数,可用于做类,用于波形显示

void CMyvirtualDlg::FFT(COMPLEX *x,int m)
{

	static COMPLEX *w;
	static int mstore = 0;
	static double n= 1;

	COMPLEX u,temp,tm;
	COMPLEX *xi,*xip,*xj,*wptr;

	int i,j,k,l,le,windex;
	double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
		
		if(m != mstore)
	{
		// free previously allocated storage and set new m
		if(mstore!=0) free(w);
		mstore = m;
		if(m==0) return; //if m=0 then done
		//n=2**m = fft length
		n=1<<m;
		le = n/2;
		//allocate the storage for w
		w = (COMPLEX *)calloc(le-1,sizeof(COMPLEX));
		if(!w)
		{
			//TRACE("\n unable to allocate complex w array\n");
		//	exit(1);
		}
		//caculate the w values recursively
		arg = 4.0*atan(1.0)/le;
		wrecur_real = w_real = cos(arg);
		wrecur_imag = w_imag = -sin(arg);
		xj = w;
		for(j=1;j<le;j++)
		{
			xj->real = (float)wrecur_real;
			xj->imag = (float)wrecur_imag;
			xj++;
			wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
			wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
			wrecur_real = wtemp_real;
		}
	}

	//start fft 
	le = n;
	windex = 1;
	for(l=0;l<m;l++)
	{
		le = le/2;
		//first iteration with no multiplies
		for(i = 0;i<n;i=i+2*le)
		{
			xi = x +i;
			xip = xi + le;
			temp.real = (xi->real + xip->real);
			temp.imag = (xi->imag + xip ->imag);
			xip->real = (xi->real - xip->real);
			xip->imag = (xi->imag - xip->imag);
			*xi = temp;
		}
		//remaining iterations use store w 
		wptr = w + windex-1;
		for(j=1;j<le;j++)
		{
			u = *wptr;
			for(i=j;i<n;i=i+2*le)
			{
				xi = x+i;
				xip = xi +le;
				temp.real = (xi->real + xip->real);
				temp.imag = (xi->imag + xip->imag);
				tm.real = xi->real - xip->real;
				tm.imag = xi->imag - xip->imag;
				xip->real = (tm.real*u.real*u.real-tm.imag*u.imag);
				xip->imag = (tm.real*u.imag+tm.imag*u.real);
				*xi = temp;
			}
			wptr = wptr + windex;
		}
		windex = 2*windex;
	}
	//rearrang data by bit reversing
	for(i=0;i<n;++i)
	{
		j=0;
		for(k=0;k<m;++k)
			j=(j<<1)|(1&(1>>k));
		if(i<j)
		{
			xi = x+i;
			xj = x+j;
			temp = *xj;
			*xj = *xi;
			*xi = temp;
		}
	}
}