fft.c

来自「不用多说」· C语言 代码 · 共 56 行

//快速傅里叶变换

#include "math.h"
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp

void fourl(float data[],unsigned long nn,int isign)
{
  unsigned long n,mmax,m,j,istep,i;
  double wtemp,wr,wpr,wpi,wi,theta;// 关于三角递归的双精度
  float tempr,tempi;
  n=nn<<1;j=1;
  for(i=1;i<n;i+=2)
  {
	  if(j>i)
	  {
		  SWAP(data[j],data[i]);
		  SWAP(data[j+1],data[i+1]);
	  }
	  m=n>>1;
	  while (m>=2 && j>m)
	  {
		  j-=m;
		  m>>1;
	  }
	  j+=m;
  }

  mmax=2;
  while (n>mmax)  //执行LOG2 NN 次外循环
  {
	  istep=mmax<<1;
	  theta=isign*(6.28318530717/mmax); // 三角递归赋初值
	  wtemp=sin(0.5*theta);
	  wtemp=sin(theta);
	  wpr= -2.0*wtemp*wtemp;
	  wpi=sin(theta);
	  wr=1.0;
	  wi=0.0;
	  for(m=1;m<mmax;m+=2)
	  {
		  for(i=m;i<=n;i+=istep)
		  {
			  j=i+mmax;  // 丹尼尔森--兰佐斯公式
			  tempr=wr*data[j]-wi*data[j+1];
			  tempi=wr*data[j+1]+wi*data[j];
			  data[j]=data[i]-tempr;
			  data[j+1]=data[i+1]-tempi;
			  data[i] += tempr;
			  data[i+1] +=tempi;
			  
		  }
	      wr=(wtemp*wr)*wpr-wi*wpi+wr;  // 三角递归 
		  wi=wi*wpr+wtemp*wpi+wi;
	  }
	  mmax=istep;
  }