dft_cplx()

void fourier::dft_cplx	(	double	w,
		int	n,
		double	a,
		double	b,
		cplx *	f,
		cplx &	res,
		cplx &	err
	)

Computes the Fourier integral \(I(\omega) = \int^b_a dt\, e^{i \omega t}f(t)\) using cubically corrected DFT.

Purpose

For a function \(f(t)\), represented on an equidistant grid \(t_j= j h\), \(j=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

dft_cplx computes the integral \(I(\omega)\), using the cubic correction factor and boundary correction terms computed by complex_dftcor_cubic. Furthermore, an approximate error is returned, given the difference using n and (n+1)/2 points.

Parameters

w	[double] frequency \(\omega\)
n	[int] number of time points \(t_j=j h, j=0,\dots,n\)
a	[double] starting point of interval
b	[double] end point of interval
f	[complex] pointer containing the function values
res	[complex] Fourier integral \(I(\omega)\)
err	[complex] difference using 'n' and '(N+1)/2' points.

Definition at line 283 of file fourier.cpp.

References complex_dftcor_cubic().

Referenced by fourier::adft_func::dft(), and fourier::adft_func::sample().

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

Here is the call graph for this function:

Here is the caller graph for this function: