fft16.c

DSP培训 作业【源码】

#include "fft16.h"
#include "math.h"

struct compx EE(struct compx b1,struct compx b2)
{
	struct compx b3;
	b3.real=b1.real*b2.real-b1.imag*b2.imag;
	b3.imag=b1.real*b2.imag+b1.imag*b2.real;
	return(b3);
}

void FFT(struct compx *xin,int N)
{
	int f,m,nv2,nm1,i,k,j=1,l;
	/*int f,m,nv2,nm1,i,k,j=N/2,l;*/
	int le,lei,ip;
	float pi;
	struct compx v,w,t;
	nv2=N/2;
	f=N;
	for(m=1;(f=f/2)!=1;m++){;}
	nm1=N-1;    

	/*变址运算*/
	for(i=1;i<=nm1;i++)
	{
		if(i<j){t=xin[j];xin[j]=xin[i];xin[i]=t;}
		k=nv2;
		while(k<j){j=j-k;k=k/2;}
		j=j+k;
	}



	for(l=1;l<=m;l++)
	{   
		le=pow(2,l);// 这里用的是L而不是1  !!!!
		lei =le/2;
		pi=3.14159;
		v.real=1.0;
		v.imag=0.0;  
		w.real=cos(pi/lei);
		w.imag=-sin(pi/lei);
		for(j=1;j<=lei;j++)
		{
  			for(i=j;i<=N;i=i+le)
			{ 
  				ip=i+lei;
     			t=EE(xin[ip],v);
  				xin[ip].real=xin[i].real-t.real;
  				xin[ip].imag=xin[i].imag-t.imag;
  				xin[i].real=xin[i].real+t.real;
  				xin[i].imag=xin[i].imag+t.imag;
  			}
			v=EE(v,w);
		}
	}
	return;
}