fourier.cpp

Go to the documentation of this file.

/*--------------------------------------------------
       fourier.c
       Martin Eckstein, Dec. 7, 2006
-----------------------------------------------------
PURPOSE:   fourier transforms of complex functions         */
#include "./fourier.hpp"
#include <cmath>

namespace fourier{

typedef std::complex<double> cplx;

void get_dftcorr_cubic(double th,double *corfac,cplx *endcor)
{
    double ai[4],ar[4],t;
    double t2,t4,t6;
    double cth,ctth,spth2,sth,sth4i,stth,th2,th4,tmth2,tth4i;
        int j;
    
    if (fabs(th) < 5.0e-2) {
        t=th;
        t2=t*t;
        t4=t2*t2;
        t6=t4*t2;
        *corfac=1.0-(11.0/720.0)*t4+(23.0/15120.0)*t6;
        ar[0]=(-2.0/3.0)+t2/45.0+(103.0/15120.0)*t4-(169.0/226800.0)*t6;
        ar[1]=(7.0/24.0)-(7.0/180.0)*t2+(5.0/3456.0)*t4-(7.0/259200.0)*t6;
        ar[2]=(-1.0/6.0)+t2/45.0-(5.0/6048.0)*t4+t6/64800.0;
        ar[3]=(1.0/24.0)-t2/180.0+(5.0/24192.0)*t4-t6/259200.0;
        ai[0]=t*(2.0/45.0+(2.0/105.0)*t2-(8.0/2835.0)*t4+(86.0/467775.0)*t6);
        ai[1]=t*(7.0/72.0-t2/168.0+(11.0/72576.0)*t4-(13.0/5987520.0)*t6);
        ai[2]=t*(-7.0/90.0+t2/210.0-(11.0/90720.0)*t4+(13.0/7484400.0)*t6);
        ai[3]=t*(7.0/360.0-t2/840.0+(11.0/362880.0)*t4-(13.0/29937600.0)*t6);
    } else {
        cth=cos(th);
        sth=sin(th);
        ctth=cth*cth-sth*sth;
        stth=2.0e0*sth*cth;
        th2=th*th;
        th4=th2*th2;
        tmth2=3.0e0-th2;
        spth2=6.0e0+th2;
        sth4i=1.0/(6.0e0*th4);
        tth4i=2.0e0*sth4i;
        *corfac=tth4i*spth2*(3.0e0-4.0e0*cth+ctth);
        ar[0]=sth4i*(-42.0e0+5.0e0*th2+spth2*(8.0e0*cth-ctth));
        ai[0]=sth4i*(th*(-12.0e0+6.0e0*th2)+spth2*stth);
        ar[1]=sth4i*(14.0e0*tmth2-7.0e0*spth2*cth);
        ai[1]=sth4i*(30.0e0*th-5.0e0*spth2*sth);
        ar[2]=tth4i*(-4.0e0*tmth2+2.0e0*spth2*cth);
        ai[2]=tth4i*(-12.0e0*th+2.0e0*spth2*sth);
        ar[3]=sth4i*(2.0e0*tmth2-spth2*cth);
        ai[3]=sth4i*(6.0e0*th-spth2*sth);
    }
    for(j=0;j<=3;j++) endcor[j]=cplx(ar[j],ai[j]);
}
void get_dftcorr_linear(double th,double *corfac,cplx *endcor)
{
    double ai,ar;
    double th2,th4,th6,cth,sth;
    
    if (fabs(th) < 5.0e-2) {
        th2=th*th;
        th4=th2*th2;
        th6=th4*th2;
        *corfac=1.0-(1.0/12.0)*th2+(1.0/360.0)*th4-(1.0/20160.0)*th6;
        ar=-0.5+th2/24.0-(1.0/720.0)*th4+(1.0/40320.0)*th6;
        ai=th*(1.0/6.0-(1.0/120.0)*th2+(1.0/5040.0)*th4-(1.0/362880.0)*th6);
    } else {
        cth=cos(th);
        sth=sin(th);
        th2=th*th;
        *corfac=2.0*(1.0-cth)/th2;
        ar=(cth-1.0)/th2;
        ai=(th-sth)/th2;
    }
    endcor[0]=cplx(ar,ai);
    endcor[1]=0;
}


void complex_dftcor_cubic(double w, double delta,double a,double b, cplx *endpts,cplx *endcor,double *corfac)
{
    double ai[4],ar[4],arg,t;
    double t2,t4,t6;
    double cth,ctth,spth2,sth,sth4i,stth,th,th2,th4,tmth2,tth4i;
        cplx expfac,cr,cl,cor,al;
    int j;
    
    th=w*delta;
    if (fabs(th) < 5.0e-2) {
        t=th;
        t2=t*t;
        t4=t2*t2;
        t6=t4*t2;
        *corfac=1.0-(11.0/720.0)*t4+(23.0/15120.0)*t6;
        ar[0]=(-2.0/3.0)+t2/45.0+(103.0/15120.0)*t4-(169.0/226800.0)*t6;
        ar[1]=(7.0/24.0)-(7.0/180.0)*t2+(5.0/3456.0)*t4-(7.0/259200.0)*t6;
        ar[2]=(-1.0/6.0)+t2/45.0-(5.0/6048.0)*t4+t6/64800.0;
        ar[3]=(1.0/24.0)-t2/180.0+(5.0/24192.0)*t4-t6/259200.0;
        ai[0]=t*(2.0/45.0+(2.0/105.0)*t2-(8.0/2835.0)*t4+(86.0/467775.0)*t6);
        ai[1]=t*(7.0/72.0-t2/168.0+(11.0/72576.0)*t4-(13.0/5987520.0)*t6);
        ai[2]=t*(-7.0/90.0+t2/210.0-(11.0/90720.0)*t4+(13.0/7484400.0)*t6);
        ai[3]=t*(7.0/360.0-t2/840.0+(11.0/362880.0)*t4-(13.0/29937600.0)*t6);
    } else {
        cth=cos(th);
        sth=sin(th);
        ctth=cth*cth-sth*sth;
        stth=2.0e0*sth*cth;
        th2=th*th;
        th4=th2*th2;
        tmth2=3.0e0-th2;
        spth2=6.0e0+th2;
        sth4i=1.0/(6.0e0*th4);
        tth4i=2.0e0*sth4i;
        *corfac=tth4i*spth2*(3.0e0-4.0e0*cth+ctth);
        ar[0]=sth4i*(-42.0e0+5.0e0*th2+spth2*(8.0e0*cth-ctth));
        ai[0]=sth4i*(th*(-12.0e0+6.0e0*th2)+spth2*stth);
        ar[1]=sth4i*(14.0e0*tmth2-7.0e0*spth2*cth);
        ai[1]=sth4i*(30.0e0*th-5.0e0*spth2*sth);
        ar[2]=tth4i*(-4.0e0*tmth2+2.0e0*spth2*cth);
        ai[2]=tth4i*(-12.0e0*th+2.0e0*spth2*sth);
        ar[3]=sth4i*(2.0e0*tmth2-spth2*cth);
        ai[3]=sth4i*(6.0e0*th-spth2*sth);
    }
    cl=0;
    for(j=0;j<=3;j++){
       al=cplx(ar[j],ai[j]);
       cor=endpts[j]*al;
       cl+=cor;
    }
    cr=0;
    for(j=0;j<=3;j++){
       al=cplx(ar[j],-ai[j]);
       cor=endpts[7-j]*al;
       cr += cor;
    }
    arg=w*(b-a);
    expfac=cplx(cos(arg),sin(arg));
    cr *= expfac;
    *endcor = cl+cr;
}

void dft_cplx(double w,int n,double a,double b,cplx *f,cplx &res,cplx &err)
{
  /*same as fourier transform complex_ft_cubic, but
  an approximate error is returned, given the difference between 
  complex_ft_cubic using m and (m+1)/2 points*/
  double theta,delta,len,corfac,corfac2,arg;
  cplx endpoints[8],endcor,endcor2,res2,expfac,z;
  int j;
  
  len=b-a;
  if(n<=16  || (n/2)*2!=n  || len<=0.0) std::cerr << "n=" << n << " len= " << len << std::endl;
  if(n<=16  || (n/2)*2!=n  || len<=0.0) {std::cerr << "complex_ft_cubic_err: wrong input"  << std::endl;abort();}
  delta=len/((double) n);
  theta=w*delta;
  /*endcorrections using all points*/
  /*store endpoints*/
  for(j=0;j<=3;j++){
    endpoints[j]=f[j];
    endpoints[7-j]=f[n-j];
  }
  /*get corrections*/
  complex_dftcor_cubic(w,delta,a,b,endpoints,&endcor,&corfac);
  /*endcorrections using every second points*/
  /*store endpoints*/
  for(j=0;j<=3;j++){
    endpoints[j]=f[2*j];
    endpoints[7-j]=f[n-2*j];
  }
  /*get corrections*/
  complex_dftcor_cubic(w,2.0*delta,a,b,endpoints,&endcor2,&corfac2);
  /*compute discrete fourier transform 
  (naive version, without FFT)*/
  /*first using every second point*/
  res2=0;
  for(j=0;j<=n;j+=2){
    arg=((double) j)*theta;
    expfac=cplx(cos(arg),sin(arg));
    res2 += expfac*f[j];
  }
  /*using remaining points*/
  res=res2;
  for(j=1;j<n;j+=2){
    arg=((double) j)*theta;
    expfac=cplx(cos(arg),sin(arg));
    res += expfac*f[j];
  }
  /*apply corrections*/
  res *= corfac;
  res += endcor;
  arg=w*a;
  expfac=cplx(delta*cos(arg),delta*sin(arg));
  res *= expfac;
  
  res2 *= corfac2;
  res2 += endcor2;
  expfac *=2;
  res2 *=  expfac;
  err=res-res2;
}

} //namespace