fft.cpp

计算FFT 的程序

#include <cmath.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <conio.h>


#include <graphics.h>

/* COMPLEX STRUCTURE */
typedef struct {
float real, imag;
} COMPLEX;

#define PI (4.0*atan(1.0))

void fft(COMPLEX *,int);
void ifft(COMPLEX *,int);
void draw_image(double *x,int m,char *title1,char *title2,
char *xdis1,char *xdis2,int dis_type);

/********************************************************/
void main(void)
{
int i,length,m,j;
char title[80],tmp[20];
double *amp;
double a,tempflt;
COMPLEX *samp;

m=8;
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++) {
amp[i] = cos(2.0*PI*(double)i*2.0/length);
printf("*");
}
for (i = 0; i < length; i++) {
samp[i].real = amp[i]*cos(16.0*PI*(double)i*2.0/length)*amp[i];
amp[i]= samp[i].real;
printf("*");
}

strcpy(title,"The Sampling Signal Data");
draw_image(amp,length,title,"The Magnitude","0",itoa(length,tmp,10),0);

/* Find the spectrum of the data */
printf("Waitting for the FFT calculation...\n");
fft(samp,m);

/* calculate magnitude */
tempflt=0;
for (i = 0 ; i < length ; i++) {
tempflt = samp[i].real * samp[i].real + samp[i].imag * samp[i].imag;
amp[i] = tempflt;
}

strcpy(title,"The Signal Freq Magnitude Result");
draw_image(amp,length,title,"The Magnitude","0",itoa(length,tmp,10),0);

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,"The Magnitude","0",itoa(length,tmp,10),0);

free(samp);
free(amp);
}

void fft(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 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_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) {
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);
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 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;
}
}
}
//////////////////////////
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;
wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
wrecur_real = wtemp_real;
}
}

/* start inverse 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 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;
}
}

/* scale all results by 1/n */
scale = (float)(1.0/n);
for(i = 0 ; i < n ; i++) {
x->real = scale*x->real;
x->imag = scale*x->imag;
x++;
}
}