ffttest.cpp

分步傅利叶方法解非线性方程

/*
    This program is to test the vality of the split-step Fourier 
	method. Please note that the fiber loss is not included.
*/
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <malloc.h>
#include <math.h>
#include <complex.h>

int const n=1024;                    //n: the sampling point
double const  PI=3.1415926;
void FFT(double *, double *, int, int);  //FFT algorithm

void main()
{	
   int     i, j, flag,k,i1,i2;
   double  fr[n], fi[n],frr[n],Tw,dt,t,nn,z0,dz,T0;
   double  fr1[20][n],z,omiga,Tr,s1,s2;
   m_complex  cc1, cc2;
   FILE    *fp,*fp1;

   Tw=40;               //Tw: the width of the sampling window in time domain
   omiga=2.0*PI/Tw;     //omiga: the radian frequency space 
   z0=PI/2.0;           //z0: the soliton period
   dz=z0/200;           //dz: the distance space 
   dt=Tw/n;             //dt: time step
   Tr=0.01;             //Tr: RSFS coefficient
   nn=4.0;              //nn: the square of the soliton order

   //sample the pulse in time domain
   for(i=0;i<n/2;i++)
   {
	   fr[i]=1.0/cosh(dt*i);
	   fi[i]=0.0;
   }
   for(i=n/2;i<n;i++)
   {
	   fr[i]=1.0/cosh(dt*(i-n));
	   fi[i]=0.0;
   }


   //fr1[][] store the data for output
   for(i=0;i<n;i++)
	   fr1[0][i]=fr[i]*fr[i]+fi[i]*fi[i]; 

   //the split-step algorithm
   j=1;
   for(z=0;z<5.0*z0+0.01;z=z+dz)
   {
	   if(fabs(z-z0*j)<1.0e-4)
	   {
		   for(i=0;i<n;i++)
			   fr1[j][i]=fr[i]*fr[i]+fi[i]*fi[i];
		   j++;
	   }

	   //FFT transform, signal in time domain--> in frequency domain
	   flag=0;
	   FFT(fr,fi,n,flag);
	   for(i=0;i<n/2;i++)
	   {
		   cc1=m_complex(0.0,-omiga*omiga*i*i*dz/2.0);
		   cc1=exp(cc1);
		   cc2=m_complex(fr[i],fi[i]);
		   fr[i]=real(cc1*cc2);
		   fi[i]=imag(cc1*cc2);
	   }
	   for(i=n/2;i<n;i++)
	   {
		   cc1=m_complex(0.0,-omiga*omiga*(i-n)*(i-n)*dz/2.0);
		   cc1=exp(cc1);
		   cc2=m_complex(fr[i],fi[i]);
		   fr[i]=real(cc1*cc2);
		   fi[i]=imag(cc1*cc2);
	   }
	   
	   //inverse FFT transform, signal in frequency domain--> in time domain
	   flag=1;
	   FFT(fr,fi,n,flag);
   	
	   //consider the effect of the RSFS
	   for(i=0;i<n;i++)
		   frr[i]=fr[i]*fr[i]+fi[i]*fi[i];
	   for(i=0;i<n;i++)
	   {
		   if(i==0)
			   i1=n-1;
		   else
			   i1=(i-1)%n;
		   if(i==n-1)
			   i2=0;
		   else
			   i2=(i+1)%n;
		   s1=frr[i1];
		   s2=frr[i2];
		   s2=(s2-s1)/(2.0*dt);
		   cc1=m_complex(0.0,nn*(frr[i]-Tr*s2)*dz);
		   cc1=exp(cc1);
		   cc2=m_complex(fr[i],fi[i]);
		   cc2=cc1*cc2;
		   fr[i]=real(cc2);
		   fi[i]=imag(cc2);
	   }
   }

   //testt.dat store the signals in time domain
   fp=fopen("d:\\temp\\testt.dat","wt");
   if(fp==NULL)
   {
	   printf("File open error!");
	   exit(1);
   }
   //testf.dat store the signals in frequency doamin
   fp1=fopen("d:\\temp\\testf.dat","wt");
   if(fp1==NULL)
   {
	   printf("File open error!");
	   exit(1);
   }

   for(i=n/2;i<n;i++)
   {
	   fprintf(fp,"%lf\t",dt*(i-n));
	   for(k=0;k<j;k++)
		   fprintf(fp,"%lf\t",fr1[k][i]+1.0*k);
	   fprintf(fp,"\n");
   }
   for(i=0;i<n/2;i++)
   {
	   fprintf(fp,"%lf\t",dt*i);
	   for(k=0;k<j;k++)
		   fprintf(fp,"%lf\t",fr1[k][i]+1.0*k);
	   fprintf(fp,"\n");
   }
   FFT(fr,fi,n,0);
   for(i=n/2;i<n;i++)
	   fprintf(fp1,"%lf\t%e\n",(i-n)/Tw,sqrt(fr[i]*fr[i]+fi[i]*fi[i]));
   for(i=0;i<n/2;i++)
	   fprintf(fp1,"%lf\t%e\n",i/Tw,sqrt(fr[i]*fr[i]+fi[i]*fi[i]));
   fclose(fp);
   fclose(fp1);
}


/*************************************************************
the FFT program--Fast Fourier Transform Program
**************************************************************/
void FFT(double *fr, double *fi, int n, int flag)
//n: the number of point
//flag=0: FFT transform
//flag=1: inverse FFT transform
//fr[]:store the real part of data for both input and output
//fi[]:store the imaginary part of data for both input and output
{
   int     mp, arg, cntr;
   int      i, j, a, b, k;
   double   sign, harm,t;
   m_complex  u, w, tt;

   j=0;
   if(flag!=0)
   {
	   sign=1.0;
	   for(i=0;i<=n-1;++i)
	   {
		   fr[i]=fr[i]/n;
		   fi[i]=fi[i]/n;
	   }
   }
   else
	   sign=-1.0;
   for(i=0;i<n-2;++i)
   {
	   if(i<j)
	   {
		   t=fr[i];
		   fr[i]=fr[j];
		   fr[j]=t;
		   t=fi[i];
		   fi[i]=fi[j];
		   fi[j]=t;
	   }
	   k=n/2;
	   while(k<=j)
	   {
		   j-=k;
		   k/=2;
	   }
	   j+=k;
   }
   mp=0;
   i=n;
   while(i!=1)
   {
	   mp+=1;
	   i/=2;
   }
   harm=2*asin(1);
   a=2;
   b=1;
   for(cntr=1;cntr<=mp;++cntr)
   {
	   u=m_complex(1.0,0.0);
	   w=m_complex(cos(harm/b),-sign*sin(harm/b));
	   for(k=0;k<=b-1;++k)
	   {
		   i=k;
		   while(i<n)
		   {
			   arg=i+b;
			   tt=m_complex(fr[arg],fi[arg])*u;
			   fr[arg]=fr[i]-real(tt);
			   fi[arg]=fi[i]-imag(tt);
			   fr[i]=fr[i]+real(tt);
			   fi[i]=fi[i]+imag(tt);
			   i+=a;
		   }
		   u=u*w;
	   }
	   a*=2;
	   b*=2;
   }
   if(flag==0)
   {
	   for(i=0;i<n;i++)
	   {
		   fr[i]=fr[i]/n;
		   fi[i]=fi[i]/n;
	   }
   }
   else
   {
	   for(i=0;i<n;i++)
	   {
		   fr[i]=fr[i]*n;
		   fi[i]=fi[i]*n;
	   }
   }
}