fft.c

快速傅立叶变换 c语言编写 要通过调试 dsp程序

#include "math.h" 
#include "stdio.h"
int		N=16;                        
int		r=4;                         
int 	i,j=1,k=0,l,m,a,b,c;
float   R1,R2,I1,I2,pi=3.1415926;
float   x_real[16];          
float   x_imag[16];
float   w_real[8];
float   w_imag[8];
float   y_real[16];
float   y_imag[16];
float   y[16];
void main()
{   
	for (i=0;i<N;i++)                
	{ 
		x_real[i]=exp(-i);
		x_imag[i]=0;
	}
	
	                                 	
	for (i=0;i<N;i++)                
	{ 
		y_real[i]=x_real[k];
		y_imag[i]=x_imag[k];
		m=N;
loop:	k=k+(m/2);
		if(k>=m)
		{
			k=k-m;
			m=m/2;
			goto    loop;
		}		
	}
	
	for (i=0;i<(N/2);i++)            
	{   
	    w_real[i]=cos(i*(2*pi)/N);
		w_imag[i]=-sin(i*(2*pi)/N);
	}
    
	for (i=0;i<r;i++)               
	{   
		if(i!=0){j=j*2;}
		m=(N/2)/j;
		for(k=0;k<m;k++)
		{
			for (l=0;l<j;l++)                
			{ 
				a=l+2*j*k;
				b=a+j;
				c=l*m;				
				R1=y_real[a];I1=y_imag[a];
				R2=y_real[b];I2=y_imag[b];
				y_real[a]=R1+R2*w_real[c]-I2*w_imag[c];
				y_imag[a]=I1+R2*w_imag[c]+I2*w_real[c];	
				y_real[b]=R1-R2*w_real[c]+I2*w_imag[c];
				y_imag[b]=I1-R2*w_imag[c]-I2*w_real[c];			
			}
		}		
	} 
	
	for (i=0;i<N;i++)              
	{ 
		y[i]=sqrt(pow(y_real[i],2)+pow(y_imag[i],2));		
	printf("y_real[%d]=%f\ty_imag[%d]=%f\n",i,y_real[i],i,y_imag[i]);  
	}
	

}