DCTI

DCTI is one of realizations of the discrete analogies of the CosFourier operator  

The name DCTII is chosen in analogy with Wikipedia [1] and notations by the Numerical recipes in C [2].

DCTI, or Discrete Cos Transform of kind I, is orthogonal transform, that repalces an array  of length  with elements ,  to the array  with elements

 for 

Normalized form

The orthonormaized transform can be defined with operator , that acts on array  in the following way:



Operator  is its own inverse; 

Numerical implementation

the C++ numerical implementation of the discrete cos transtorm of First kind consists of 3 files zfour1.cin, zrealft.cin, zcosft1.cin; these files should be loaded to the working directory in order to compile the examples. For the application in wave optics, z_type should be defined as double" or complex(double); however, for other applications, such a type may be defined in other ways too. The name of the functions and sense of the arguments are chosen following notations by the Numerical recipes in C, the call of the transform of array  of length  has form zcosft1(F-1,N); after such a call, values of the elements of array  are replaced with values calculated with the expression (1) above. For , the evaluation requires of order of  operations

Approximation of the CosFourier

The CosFourier operator transforms a function  of non–negative argument to function  in the following way:



For the discrete approximation of this operator, assume some large natural number . Let . Let function  be smooth and quickly decay at infinity. Then, the transform of  can be approximated as follows:



For , this can be written as follows:

 

At the transformation, it is assumed, that  can be neglected as . In such a way,



Numerical test of approximation of CosFourier

Files zfour1.cin, zrealft.cin, zcosft1.cin should be loaded to the working directory in order to compile the testfile below.

The example numerical implementation of the CosFourier transform with approximation described in previous section is suggested below. The C++ numerical CosFourier transform of the self-Fourier function  can be realized as follows:

#include<math.h>
#include<stdio.h>
#include <stdlib.h>
//#include <complex>
//using namespace std;
#define z_type double
#include"zfour1.cin"
#include"zrealft.cin"
#include"zcosft1.cin"
main(){ z_type *a, *b; int j,k, N=16; double d,x,y;
a=(z_type *) malloc((size_t)((N+1)*sizeof(z_type)));
b=(z_type *) malloc((size_t)((N+1)*sizeof(z_type)));
d=sqrt(M_PI/N); for(j=0;j<N+1;j++){ x=j*d; a[j]=b[j]=exp(-x*x/2.);}
zcosft1(a-1,N);  for(j=0;j<N+1;j++) a[j]*=sqrt(2./N);
for(j=0;j<N+1;j++) printf("%12.9f %12.9f %12.9f %11.4e\n",j*d, a[j],b[j], a[j]-b[j]);
free(a); free(b);
}


The code above generates the following output:

0.000000000  1.000000000  1.000000000 -2.3238e-12
0.443113463  0.906490462  0.906490462  2.3206e-12
0.886226925  0.675231907  0.675231907 -2.3094e-12
1.329340388  0.413303564  0.413303564  2.2924e-12
1.772453851  0.207879576  0.207879576 -2.2688e-12
2.215567314  0.085917370  0.085917370  2.2417e-12
2.658680776  0.029179416  0.029179416 -2.2100e-12
3.101794239  0.008143268  0.008143268  2.1780e-12
3.544907702  0.001867443  0.001867443 -2.1444e-12
3.988021165  0.000351903  0.000351903  2.1124e-12
4.431134627  0.000054491  0.000054491 -2.0815e-12
4.874248090  0.000006933  0.000006933  2.0543e-12
5.317361553  0.000000725  0.000000725 -2.0309e-12
5.760475015  0.000000062  0.000000062  2.0121e-12
6.203588478  0.000000004  0.000000004 -1.9832e-12
6.646701941  0.000000000  0.000000000  2.4651e-12
7.089815404  0.000000000  0.000000000  1.0175e-11


The 0th column represents the coodinate , the following two– its DTFI and the input function, and the last shows the error of the approximation of the CosFourier operator with the DTFI routine.

The example shows, that, for the simplest self-Fourier function, 17-nodes approximation gives of order of 12 correct decimal digits.

For the analysis of the code above, it should be noted, that the self-Fourier function is eigenfunction of the Fourier operator with eigenvalue 1;

Evaluation of the Fourier-coefficients

Let  be even periodic function of real argument with period :

, 

Then, the expansion into the Fourier series can be written as



The Fourier–coefficients


, 

Assume some large natural number . Let . For approximation of coefficeins , replace the integral with the finite sum:


, 

Comparison to equation (1) gives


, 

Given expansion coefficients , the function at the mesh  can be evaluated with , with  and  for .

Numerical test of the expansion to the Fourier series

Let . Let . The expansion coefficients are expected to be , , ; all other coefficients are expected to be zero. The approximation above can be implemented in C++ with the following code:

#include<math.h>
#include<stdio.h>
#include <stdlib.h>
#define z_type double
#include"zfour1.cin"
#include"zrealft.cin"
#include"zcosft1.cin"
main(){ z_type *a, *b, *c; int j,k, N=8; double d,x,y;
a=(z_type *) malloc((size_t)((N+1)*sizeof(z_type)));
b=(z_type *) malloc((size_t)((N+1)*sizeof(z_type)));
c=(z_type *) malloc((size_t)((N+1)*sizeof(z_type)));
d=M_PI/N;
for(j=0;j<N+1;j++){ x=j*d; a[j]=b[j]=1.+.1*cos(x)+.01*cos(2*x);}
zcosft1(a-1,N); for(j=0;j<N+1;j++) c[j]=a[j];
zcosft1(a-1,N);
for(j=0;j<N+1;j++) printf("%2d %12.9f %12.9f %12.9f\n",j, b[j],c[j],a[j]);
free(a); free(b);
}


The code above can be compiled with files zfour1.cin, zrealft.cin, zcosft1.cin and generates the output below:

0  1.110000000  8.000000000  4.440000000
1  1.099459021  0.400000000  4.397836084
2  1.070710678  0.040000000  4.282842712
3  1.031197275 -0.000000000  4.124789102
4  0.990000000  0.000000000  3.960000000
5  0.954660589 -0.000000000  3.818642356
6  0.929289322 -0.000000000  3.717157288
7  0.914683115 -0.000000000  3.658732458
8  0.910000000  0.000000000  3.640000000


The 0th column just numerates the nodes of the grid.

The 1st column shows the initial function.

The 2d column represents the DCTI transform of the initial function. The number in this column are the Fourier coefficients used in the initial function, multiplied with factor . The last column is result of the additional application of the operator DCTI and represents the same initial values of the function multiplied with the same factor.

Conclusion

 can be used for evaluation of the CosFourier operator at the equidistant array of values of function, assuming that the function is continuous and smoothly decays at infinity. The array should have  elements, and the numeration sould begin with zero.

The same discrete operator can be used for the evaluation of the Fourier coefficients of a symmetric periodic function, as well as for the evaluation of a function by given Fourier coefficients with truncated Fourier series.

The numerical implementation is especially efficient for  for some natural number . The C++ implementation is stored in routines zfour1.cin, zrealft.cin, zcosft1.cin; they can be copied and included to the user's code without any "installing". In order to deal with real numbers, type z_type can be defined as "double"; in order to deal with complex numbers, it can be defined as complex(double).

References

1. http://en.wikipedia.org/wiki/Discrete_cosine_transform
2. http://88.167.97.19/albums/files/TMTisFree/Documents/Physics/11%20-%20Fourier%20Transform%20Spectral%20Methods.pdf W.H.Press, B.P.Flannery, S.A.Teukolsky, W.T.Vetterling. Numerical Recipes in C. Fast Sine and Cosine transform.