myfft.cpp

根据菲涅耳衍射原理编写的光学衍射模拟程序

/**###**************************************************************###**/
//注    意：此软件允许无偿使用、拷贝、修改、分发，并可用于任何用途，仅当您同意如下条款：
//          在使用此软件的任何文挡中包含下面的所有声明，以及本注意事项！ 
//Caution ：Permission to use, copy, modify, and distribute this software
//          in code form for any purpose and without fee is hereby granted, 
//          provided that the above notice ,below declare and this limited warranty 
//          appears in all copies and all supporting documentation.
//作    者：王寅庆    E_mail:[e-mail]
//指导教师：李祥教授
//部    门：贵州大学计算机软件和理论研究所 贵州大学
//地    址：贵州省 贵阳市花溪 (550025)
//时    间：2003-1-20 19:51:18
//注    释：
/**###**************************************************************###**/

#include <math.h>
#include <malloc.h>
#include "myfft.h"

//fft的实现文件
//struct Complex {double re,im;};
void fft2(struct Complex ft[],int m,int inv)
{ int n,i,j,k,mm,mp;
  double s,t,ang,tpi=6.283185307179586;
  struct Complex u,w,*p,*q,*pf;
  n=1; n<<=m; pf=ft+n-1;
  for(j=0,p=ft; p<pf ;++p){ q=ft+j;
    if(p<q){ t=p->re; p->re=q->re; q->re=t;
             t=p->im; p->im=q->im; q->im=t; }
    for(mm=n/2; mm<=j ;mm/=2) j-=mm; j+=mm;
   }
  if(inv!='d') for(p=ft,s=1./n; p<=pf ;){p->re*=s; (p++)->im*=s; }
  for(i=mp=1; i<=m ;++i){
    mm=mp; mp*=2; ang=tpi/mp; if(inv=='d') ang= -ang;
    w.re=cos(ang); w.im=sin(ang);
    for(j=0; j<n ;j+=mp){ p=ft+j;
      u.re=1.; u.im=0.;
      for(k=0; k<mm ;++k,++p){ q=p+mm;
        t=q->re*u.re-q->im*u.im;
        s=q->im*u.re+q->re*u.im;
        q->re=p->re-t; q->im=p->im-s;
        p->re+=t; p->im+=s;
        t=u.re*w.re-u.im*w.im;
        u.im=u.im*w.re+u.re*w.im; u.re=t;
       }
     }
   }
}
//////////////////
void fft2_2d(struct Complex a[],int m,int n,int f)
{ register int md,nd,i,j; struct Complex *p,*q;
  register struct Complex *r,*s;
  md=1<<m; nd=1<<n;
  for(i=0,p=a; i<md ;++i){
    fft2(p,n,f); p+=nd;
   }
  q=(struct Complex *)malloc(sizeof(*q)*md);
  for(i=0,p=a-nd; i<nd ;++i){
    for(r=q,s=p,j=0; j<md ;++j) *r++ = *(s+=nd);
    fft2(q,m,f);
    for(r=q,s=p++,j=0; j<md ;++j) *(s+=nd)= *r++;
   }
  free(q);
}