2-一维快速傅氏变换-p47.txt

地球物理大师Claebout的几个有用的程序

  #1D fast Fourier transform
  subroutine fork(lx,cx,signi)
  complex cx(lx),carg,cexp,cw,ct
  j=1; k=1;    sc=sqrt(1./lx)
  do i=1,lx {
      if(i<=j) { ct=cx(j) *sc;    cx(j)=cx(i) *sc; cx(i)=ct }
      m=lx/2
      while(j>m) { j=j-m; m=m/2; if(m<1)break }
      j=j+m
  }
  repeat {
      istep=2*k
      do m=1,k
         carg=(0.,1.)*(3.14159265*signi*(m-1))/k; cw=cexp(carg)
         do i=m,lx,istep
            {  ct=cw*cx(i+k); cx(i+k)=cx(i)-ct; cx(i)=cx(i)+ct }
  }
  k=istep
  } until(k>lx)

relurn; end

!lx是2的整数幂，输出cx(1)是零频率，cx(lx/2+1)是所谓的Nyquist频率，而cx(lx)则是最小负频率
!signi=+1正变换;signi=-1反变换