fft2.txt

来自「快速付里叶变换子程序所需 RAM 空间以输入的首地址为基址,向增加的方向扩展」· 文本 代码 · 共 70 行

/* fftdi.c: fast Fourier transform, direct and inverse              */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int inverseq( int i1, int imax );
extern void f(double tt, double fri[2] );
void creata( double a, double b, double A[2][N] );
void creata( double a, double b, double A[2][N] )
{
int i;
double tt,fri[2],deltat;
deltat=(b-a)/N;
for(i=0;i<N;i++)
  {
  tt=a+i*deltat;
  f( tt, fri );
  A[0][i]=fri[0]*deltat;
  A[1][i]=fri[1]*deltat;
  }
}
/********************************************************************/
void fftdi( int sign, double A[2][N], double x[2][N] )
{
int i,j,i0,ii,iloop,step[M],loop[M],least,jstep,reseq;
double pi,unit,alpha,w[2],ar,ai;
if( sign == 1 )
    for(i=0;i<N;i++)  { A[0][i]=x[0][i]; A[1][i]=x[1][i]; }
step[0]=N/2;
loop[0]=1;
pi=4.0*atan(1.0);
unit=pi*2.0/((double)N);
for(i=1;i<M;i++)  { step[i]=step[i-1]/2; loop[i]=loop[i-1]*2; }
for(i=0;i<M;i++)
  {
  for(iloop=0;iloop<loop[i];iloop++)
    {
    reseq=inverseq(iloop,i);
    alpha=unit*step[i]*reseq;
    least=2*step[i]*iloop;
    w[0]=cos(alpha);
    w[1]=sign*sin(alpha);
    for(j=0;j<step[i];j++)
      {
      jstep=least+step[i]+j;
      ar=A[0][jstep]*w[0]-A[1][jstep]*w[1];
      ai=A[0][jstep]*w[1]+A[1][jstep]*w[0];
      A[0][jstep]=A[0][least+j]-ar;
      A[1][jstep]=A[1][least+j]-ai;
      A[0][least+j]=A[0][least+j]+ar;
      A[1][least+j]=A[1][least+j]+ai;
      }
    }
  }
for(i=0;i<N;i++)
  { j=inverseq(i,M);  x[0][i]=A[0][j];  x[1][i]=A[1][j]; }
if( sign == 1 )
  for(i=0;i<N;i++)
    { A[0][i]=x[0][i]/((double)N); A[1][i]=x[1][i]/((double)N); }
}
/********************************************************************/
int inverseq( int i1, int imax )
{
int i,j,ii,i0,inv;
i0=i1;  inv=0;
for(j=0;j<imax;j++)
  { ii=i0%2; inv=2*inv+ii; i0=(i0-ii)/2; }
return(inv);
}