fft.c

FFT,单片机求FFT程序

/**
 *	FFT(IFFT) - Fast Fourier transform. The length of X must be a 
 *	power of two, for a fast radix-2 fast-Fourier transform algorithm 
 *	is used.	spadger@bmy <[e-mail]> 2007.9.2
 */

#include <math.h>
#include <stdio.h>

//#include "fft.h"
// fft.h
#if 1	
typedef double real;
//#define	real float
typedef struct { real r,i; } cplx_t;
void cplx_exp(cplx_t *x, cplx_t *r);
void cplx_mul(cplx_t *x, cplx_t *y, cplx_t *r);
void bit_reverse(cplx_t *x, int N);
void fft(cplx_t *x, int N);
void ifft(cplx_t *x, int N);
#endif

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

void cplx_exp(cplx_t *x, cplx_t *r)
{
	real expx=exp(x->r);	
	r->r=expx*cos(x->i);
	r->i=expx*sin(x->i);
}

void cplx_mul(cplx_t *x, cplx_t *y, cplx_t *r)
{
	r->r=x->r*y->r-x->i*y->i;
	r->i=x->r*y->i+x->i*y->r;
}

void bit_reverse(cplx_t *x, int N)
{
	real t;
	cplx_t tmp;
	unsigned int i=0,j=0,k=0;
	for(i=0; i<N; i++) {
		k=i;
		j=0;
		t=log(0.0+N)/log(2.0);
		while((t--)>0) {
			j<<=1;
			j|=k&1;
			k>>=1;
		}
		if(j>i) {
			tmp=x[i];
			x[i]=x[j];
			x[j]=tmp;
		}
	}
}

void fft(cplx_t *x, int N)
{
	cplx_t u,d,p,W,tmp;
	int i=0,j=0,k=0,l=0,M=floor(log(0.0+N)/log(2.0));
	if(log(0.0+N)/log(2.0)-M > 0){
		printf("The length of x (N) must be a power of two!!!\n");
		return;
	}
	bit_reverse(x,N);
	for(i=0; i<M; i++) {
		l=1<<i;
		for(j=0; j<N; j+=2*l ) {
			for(k=0; k<l; k++) {
				tmp.r=0.0;
				tmp.i=-2*M_PI*k/2/l;
				cplx_exp(&tmp,&W);
				cplx_mul(&x[j+k+l],&W,&p);
				u.r=x[j+k].r+p.r;
				u.i=x[j+k].i+p.i;
				d.r=x[j+k].r-p.r;
				d.i=x[j+k].i-p.i;
				x[j+k]=u;
				x[j+k+l]=d;
			}
		}
	}
}

void ifft(cplx_t *x, int N)
{
	unsigned int i=0;
	for(i=0;i<N;i++)
		x[i].i=-x[i].i;
	fft(x,N);
	for(i=0;i<N;i++) {
		x[i].r=x[i].r/(N+0.0);
		x[i].i=-x[i].i/(N+0.0);
	}
}

/**
 *	for test and demonstation
 */  
#if 1
#define DATA_LEN 64
int main()
{
	int i;
	cplx_t x[DATA_LEN];
	for(i=0;i<DATA_LEN;i++){
		x[i].r=i;
		x[i].i=i;
	}
	
	printf("Before...\nReal\t\tImag\n");
	for(i=0;i<DATA_LEN;i++)
		printf("%f\t%f\n",x[i].r,x[i].i);

	//cplx_exp(&x[1],&x[0]);
	//bit_reverse(x,8);
	fft(x,DATA_LEN);

	printf("After...\nReal\t\tImag\n");
	for(i=0;i<DATA_LEN;i++)
		printf("%f\t%f\n",x[i].r,x[i].i);

	ifft(x,DATA_LEN);

	printf("After...\nReal\t\tImag\n");
	for(i=0;i<DATA_LEN;i++)
		printf("%f\t%f\n",x[i].r,x[i].i);

	return 0;
}
#endif