fft.cpp

非常不错的fft的算法

/***********************************************************************
FFT_FREQ.C - SOURCE CODE FOR DISCRETE FOURIER TRANSFORM FUNCTIONS
fft        In-place radix 2 decimation in time FFT
ifft       In-place radix 2 decimation in time inverse FFT
draw-image Draw the inmage of input data
***********************************************************************/
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <conio.h>
#include <graphics.h>
#include <dos.h>
/* function prototypes and structures for fft functions */
/* COMPLEX STRUCTUREN */
typeddf struct
{  float real,imag;
}COMPLEX;
#define PI 4.0 *atan(1.0)
void fft(COMPLEX *,int);
void ifft(COMPLEX *,int);
void draw_image(double *,int m,char *title,char *xdis);
/***********************************************************************/
void main(void)
{
int    i,length,m;
char   title[80],tmp[20];
struct time,timeps,timepe;
double *amp;
double a,tempflt;
COMPLEX * samp;
m=10;
length=1<<m;
amp=(double *)calloc(length+1,sizeof(double));
samp=(COMPLEX *)calloc(length+1,sizeof(COMPLEX));
    if(!samp)
    {printf("\nUnable to allocate complex array for fft\n");
    exit(1);
    }
/* Input sampling data for processing */
printf ("Waitting for sampling data...");
for(i=0;i<length;i++)
{
  samp[i].real=sin(8 * PI * i / length);
  amp[i]=samp[i].real;
}
for(i=length/2-10;i<=length/2+10;i++)
{
	samp[i].real=1;
	amp[i]=samp[i].real;
}
strcpy(title,"The Sampling Signal Data");
draw_image(amp,length,title,itoa(length,tmp,10));
/* Find the spectrum of the data */
printf("Waitting for the FFT calculation...\n");
gettime(&timeps);
fft(samp,m);
gettime(&tmepe);
printf("start time:%d:%d:%d.%d\n",timeps.ti_hour,timeps.ti_min,timeps.ti_sec,timeps.ti_hund);
printf("end time:%d:%d:%d.%d\n",timepe.ti_hour,timepe.ti_min,timepe.ti_sec,timepe.ti_hund);
printf("Press any key to continue....\n");
getch();
/* calculate magnitude */
for(i=0;i<length;i++)
{
	tempflt = samp[i].real * samp[i].real;
	tempflt + = samp[i].imag * samp[i].imag;
	amp[i] = sqrt(tempflt)/length;
}
/* Copy a section of the spectrum to the magnitude image. */
strcpy(title,"The Signal Spectrum Result");
draw_image(amp,length,title,itoa(length,tmp,10));
printf("Waitting for the IFFT calculation...\n");
ifft(samp,m);
strcpy(title,"The IFFT Result");
for(i=0;i<length;i++)amp[i] = samp[i].real;
draw_image(amp,length,title,itoa(length,tmp,10));
free(samp);
free(amp);
}
/***********************************************************************
fft- In -place radix 2 decimation in time FFT
Requires pointer to complex array x ang power of 2 size of FFT m
(size of FFT =2 ^ m). Places FFT output on top of input COMPLEX array.
void fft(COMPLEX *x,int m)
**********************************************************************/
void fft(COMPLEXS *x,int m)
{
	static COMPLEX *w;       /* used to store the w complex array */
	static int mstore = 0    /* stores m for future reference */
	static int n = 1 ;       /* length of fft stored for future */
	COMPLEX u,temp,tm;
	COMPLEX *xi,*xip,*xj,*wptr;
	int i,j,k,l,le,windex;
	double arg,w_real,w_image,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)
		{
		 printf("\n Unable to allocate complex W array\n");
		 exit(1);
		}
	/* calculate the w values recursively */
	arg = 4.0*atan(1.0)/le; /* PI/le calculation */
	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 nultiplies */
		 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 stored 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 - tm.imag * u.imag;
			 xip->imag = tm.real * u.imag + tm.imag * u.real;
			 *xi=temp;
			}
		     wptr = wptr + windex;
		    }
		 windex = 2 * windex;
		}
	/* rearrange data by bit reversing */
	j = 0;
	for(i = 1;i < (n-1);i++)
	   {
	    k = n/2;
	    while(k <= j)
		 {
		  j = j - k;
		  k = k/2;
		 }
	    j = j + k;
	    if(i < j)
	      {
	       xi = x + i;
	       xj = x + j;
	       temp = *xj;
	       *xj = *xi;
	       *xi = temp;
	      }
	   }
}
/**************************************************************************
ifft - In - place radix 2 decimation in time inverse FFT
Requires pointer to complex array,x and power of 2 size of FFT,m
(size of FFT =2^m).Places inverse FFT output on top of input
frequency domain COMPLEX array.
void ifft(COMPLEX **,int m)
**************************************************************************/
void ifft(x,m)
   COMPLEX *x;
   int m;
   {
    static COMPLEX *w;            /*used to store the w complex array */
    static int mstore = 0;        /* stores m for future reference */
    static int n = 1;             /* length of ifft stored for future */
    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;
    float scale;
    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=inverse fft length */
	n = 1 << m;
	le = n/2;
    /* allocate the storage for w */
	w=(COMPLEX *)calloc(le - 1,sizeof(COMPLEX));
	if(!w)
		{
		printf("\nUnable to allocate complex W array\n");
		exit(1);
		}
	/* calculate the w values recursively */
	
	arg = 4.0 * atan(1.0)/le;   /* PI/le calculation */
	wrecur_real = w_real = cos(arg);
	wrecur_imag = w_imag = sin(arg);   /* opposite sign from fft */
	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;
		wtemp_imag = wrecur_real * w_imag - wrecur_imag * w_real;
		wrecur_real = wtemp_real;
		}
	}
	/* start inverse fft */

	le = n;
	windex = 1;
	for(l = 0;i < 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 stored */
		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 - tm.imag * u.imag;
				  xip->imag = tm.real * u.imag + tm.imag * u.real;
				  *xi = temp;
				}
			 wptr =wptr + windex;
			}
		windex = 2 * windex;
		}
	/* rearrange data by bit reversing */
		j = 0;
		for(i = 1;i < (n-1);i++)
			{
			 k = n/2;
			 while(k <= j)
				{
			     j = j - k;
				 k = k/2;
				}
			 j = j + k;
			 if(i < j)
				{
				 xi = x + i;
				 xj = x + j;
				 temp = *xj;
				 *xj = *xi;
				 *xi = temp;
				}
			}
	/* scale all results bu 1/n */
		scale = (float)(1.0/n);
		for(i = 0;i < n;i++)
			{
			 x->real = scale * x->real;
			 x->imag = scale * x->imag;
			 x++;
			}

   }