 NESSi  v1.0.2 The NonEquilibrium Systems Simulation Library

## ◆ 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().

284 {
285  /*same as fourier transform complex_ft_cubic, but
286  an approximate error is returned, given the difference between
287  complex_ft_cubic using m and (m+1)/2 points*/
288  double theta,delta,len,corfac,corfac2,arg;
289  cplx endpoints,endcor,endcor2,res2,expfac,z;
290  int j;
291
292  len=b-a;
293  if(n<=16 || (n/2)*2!=n || len<=0.0) std::cerr << "n=" << n << " len= " << len << std::endl;
294  if(n<=16 || (n/2)*2!=n || len<=0.0) {std::cerr << "complex_ft_cubic_err: wrong input" << std::endl;abort();}
295  delta=len/((double) n);
296  theta=w*delta;
297  /*endcorrections using all points*/
298  /*store endpoints*/
299  for(j=0;j<=3;j++){
300  endpoints[j]=f[j];
301  endpoints[7-j]=f[n-j];
302  }
303  /*get corrections*/
304  complex_dftcor_cubic(w,delta,a,b,endpoints,&endcor,&corfac);
305  /*endcorrections using every second points*/
306  /*store endpoints*/
307  for(j=0;j<=3;j++){
308  endpoints[j]=f[2*j];
309  endpoints[7-j]=f[n-2*j];
310  }
311  /*get corrections*/
312  complex_dftcor_cubic(w,2.0*delta,a,b,endpoints,&endcor2,&corfac2);
313  /*compute discrete fourier transform
314  (naive version, without FFT)*/
315  /*first using every second point*/
316  res2=0;
317  for(j=0;j<=n;j+=2){
318  arg=((double) j)*theta;
319  expfac=cplx(cos(arg),sin(arg));
320  res2 += expfac*f[j];
321  }
322  /*using remaining points*/
323  res=res2;
324  for(j=1;j<n;j+=2){
325  arg=((double) j)*theta;
326  expfac=cplx(cos(arg),sin(arg));
327  res += expfac*f[j];
328  }
329  /*apply corrections*/
330  res *= corfac;
331  res += endcor;
332  arg=w*a;
333  expfac=cplx(delta*cos(arg),delta*sin(arg));
334  res *= expfac;
335
336  res2 *= corfac2;
337  res2 += endcor2;
338  expfac *=2;
339  res2 *= expfac;
340  err=res-res2;
341 }
std::complex< double > cplx
Definition: fourier.cpp:11
void complex_dftcor_cubic(double w, double delta, double a, double b, cplx *endpts, cplx *endcor, double *corfac)
Returns correction factor and evaluate boundary correction terms needed for cubically corrected DFT...
Definition: fourier.cpp:182 Here is the call graph for this function: Here is the caller graph for this function: