complex_dftcor_cubic()

void fourier::complex_dftcor_cubic	(	double	w,
		double	delta,
		double	a,
		double	b,
		cplx *	endpts,
		cplx *	endcor,
		double *	corfac
	)

Returns correction factor \(W(\theta)\) and evaluate boundary correction terms needed for cubically corrected DFT.

Purpose

For a function \(f(t)\), represented on an equidistant grid \(t_n= n h\), \(n=0,\dots,N\), we would like to compute the Fourier integral \(I(\omega) = \int^b_a dt\, e^{i \omega t}f(t)\). Using piecewie cubic interpolation of \(f(t)\), the \(I(\omega)\) can be computed by by cubically corrected discrete Fourier transformation (DFT). The algorithm is explained in

W. H. Press, S. A. Teukolosky, W. T. Vetterling, B. P. Flannery, Numerical Recipes 3rd Edition: The Art of Scientific Computing, Cambridge University Press, 2007, chapter 13.9

complex_dftcor_cubic computes the required correction factor \(W(\theta)\) and the boundary terms \(\alpha_j(\theta)f_j\), \(j=0,\dots,3\), and \( e^{i\omega(b-a)}\alpha^*_j(\theta)f_{N-j} \) where \(\theta = \omega h\).

Parameters

w	[double] frequency \(\omega\)
delta	[double] grid spacing \(h\)
a	[double] starting point of interval
b	[double] end point of interval
endpts	[complex] 'endpts[j]', 'j=0,1,2,3', corresponds to \(f_j\), while 'endpts[7-j]', 'j=0,1,2,3', corresponds to \(f_{N-j}\)
endcor	[complex] on return, boundary correction terms, stored in the order as 'endpts'
corfac	[double] on return, correction factor \(W(\theta)\)

Definition at line 182 of file fourier.cpp.

Referenced by dft_cplx().

{
    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;
}

Here is the caller graph for this function: