fft1.c

FFT傅立叶变换,C++BUILDER.

//---------------------------------------------------------------------------
void CFFT(float *RBUF,int N,int INV)
{
	float   UC,US,WC,WS,TC,TS,UC1;
	int     I,J,K,L,LE,LE1,M,NV2,NM1,IP;


	if(INV==0) { 
	  for(I=1;I<=N;I++) { XC[I]=RBUF[I]; XS[I]=0.; }
	  
	
	  }
	else { 
	  for(I=1;I<=N/2;I++) {
	    XC[I]=RBUF[I*2-1];
	    XS[I]=RBUF[I*2];
	    }
	  for(I=(N/2+1);I<=N;I++) { XC[I]=XS[I]=0.; }
	  }
	XC[0]=XS[0]=0.;

	M=10;
	if(N==2048) M=11;
	if(N==4096) M=12;
	if(N==8192) M=13;
	if(N==16384) M=14;
	if(N==32768) M=15;
	if(N==65536) M=16;

	NV2=N/2;
	NM1=N-1;
	J=1;
	for(I=1;I<=NM1;I++) {
	  if(I>=J) goto L10;
	  TC=XC[J];    TS=XS[J];
	  XC[J]=XC[I]; XS[J]=XS[I];
	  XC[I]=TC;    XS[I]=TS;
L10:  K=NV2;
L20:  if(K>=J) goto L30;
	  J=J-K;
	  K=K/2;
	  goto L20;
L30:  J=J+K;
	  }

	for(L=1;L<=M;L++) {
	  LE=pow((double)2.,(double)L);
	  LE1=LE/2;
	  UC=1.;
	  US=0.;
	  WC=cos(PAI/LE1);
	  WS=-sin(PAI/LE1);
	  if(INV!=0) WS=-WS;
	  for(J=1;J<=LE1;J++) {
	    for(I=J;I<=N;I+=LE) {
	      IP=I+LE1;
	      TC=XC[IP]*UC-XS[IP]*US; TS=XS[IP]*UC+XC[IP]*US;
	      XC[IP]=XC[I]-TC;        XS[IP]=XS[I]-TS;
	      XC[I]=XC[I]+TC;         XS[I]=XS[I]+TS;
	      }
	    UC1=UC*WC-US*WS;
	    US=US*WC+UC*WS;
	    UC=UC1;
	    }
	  }

	if(INV==0) {
	  for(I=1;I<=N/2;I++) {
	    RBUF[I*2-1]=XC[I];
	    RBUF[I*2]=XS[I];
	    }
	  }
	else for(I=1;I<=N;I++) RBUF[I]=XC[I]*2./(float)N;


}