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
The orthonormaized transform can be defined with operator , that acts on array in the following way:
Operator is its own inverse;
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
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,
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;
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 .
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.
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).
The content of this article is adopted from http://tori.ils.uec.ac.jp/TORI/index.php/DCTI