yin_fft.c.txt

我的数字信号处理 作业。FFT的C语言实现。运行时将2个c文件放在同一目录

/*         文件名为yin_fft.c,功能为计算FFT的子函数，以便被其他程序调用       */
#include "math.h"
#include "stdio.h"

void change(float xr[],float xi[],int NN){
  int i,j,k,lh;
  double t;
      j=NN/2;

for(i=1;i<NN-1;i++){
  if(i>j){
    t=xr[i];
    xr[i]=xr[j];
    xr[j]=t;
    t=xi[i];
    xi[i]=xi[j];
    xi[j]=t;
    }
  k=NN/2;

  while(j>=k){
    j=j-k;
    k=k/2;
    }
  j=j+k;
/*  printf("j=%-3d",j);
*/  }

  }



void yin_fft(float xr[],float xi[],int N){
  int L,M,B,P,i,j,k;
  double rpart,ipart;
  double rwn[64],iwn[64];
  double pi2=8.0*atan(1.0);

  M=(int)(log10(N)/log10(2));
  for(i=0;i<N;i++){
    rwn[i]=cos(i*pi2/N);
    iwn[i]=-sin(i*pi2/N);
    }

  change(xr,xi,N);

  for(L=1;L<=M;L++){
    B=(int)(pow(2,(L-1)));
    for(j=0;j<B;j++){
      P=j*((int)(pow(2,(M-L))));
      for(k=j;k<N;k+=(int)(pow(2,L)) ){
	rpart=xr[k+B]*rwn[P]-xi[k+B]*iwn[P];
	ipart=xr[k+B]*iwn[P]+xi[k+B]*rwn[P];
	xr[k+B]=xr[k]-rpart;
	xi[k+B]=xi[k]-ipart;
	xr[k]=xr[k]+rpart;
	xi[k]=xi[k]+ipart;
	}
       }
      }
}