DCTI

DCTI is one of realizations of the discrete analogies of the CosFourier operator $\displaystyle (\mathrm {CosFourier} ~F)(x)$ $={\sqrt {\frac {2}{\pi }}}\int _{0}^{\infty }\cos(xy)~F(y)~\mathrm {d} y$

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

DCTI, or Discrete Cos Transform of kind I, is orthogonal transform, that repalces an array $F$ of length $N\!+\!1$ with elements $F_{n}~$, $~n=0..N$ to the array $G$ with elements

$\!\!\!\!\!\!\!\!\!\!(1)~~~\displaystyle G_{k}=(\mathrm {DCTI} _{N}F)_{k}={\frac {1}{2}}F_{0}+{\frac {(-1)^{k}}{2}}F_{N}+\sum _{n=1}^{N-1}F_{n}\cos \!\left({\frac {\pi }{N}}nk\right)~~$ for $~~k=0,..N$

Normalized form

The orthonormaized transform can be defined with operator $\Phi _{\mathrm {C} 1,N}$, that acts on array $F$ in the following way:

$\!\!\!\!\!\!\!\!\!\!(2)~~~~(\Phi _{\mathrm {C} 1,N}~F)_{k}={\sqrt {\frac {2}{N}}}~(\mathrm {DCT1} _{N}~F)_{k}={\sqrt {\frac {2}{N}}}\left({\frac {1}{2}}F_{0}+{\frac {(-1)^{k}}{2}}F_{N}+\sum _{n=1}^{N-1}F_{n}\cos \left({\frac {\pi }{N}}nk\right)\right)$

Operator $\Phi _{\mathrm {C} 1,N}$ is its own inverse; $(\Phi _{\mathrm {C} 1,N})^{2}~F=F$

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 $F$ of length $N\!+\!1$ has form zcosft1(F-1,N); after such a call, values of the elements of array $F$ are replaced with values calculated with the expression (1) above. For $N=2^{q}$, the evaluation requires of order of $Nq$ operations

Approximation of the CosFourier

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

$\!\!\!\!\!\!\!\!\!\!\!(3)~~~\displaystyle G(x)=\mathrm {CosFourier} F(x)={\sqrt {\frac {2}{\pi }}}\int _{0}^{\infty }~\cos(xy)~F(y)~\mathrm {d} y$

For the discrete approximation of this operator, assume some large natural number $N$. Let $x_{n}={\sqrt {\pi /N}}~n$. Let function $F$ be smooth and quickly decay at infinity. Then, the transform of $F$ can be approximated as follows:

$\!\!\!\!\!\!\!\!\!\!\!(4)~~~\displaystyle G(x)=\mathrm {CosFourier} F(x)\approx {\sqrt {\frac {2}{\pi }}}\left({\frac {1}{2}}F(0)+\sum _{n=1}^{N}~\cos(xx_{n})~F(x_{n})\right)~{\sqrt {{\pi }/N}}$

For $x=x_{m}$, this can be written as follows:

$\!\!\!\!\!\!\!\!\!\!\!(5)~~~G_{m}=\displaystyle G(x_{m})=\mathrm {CosFourier} F(x_{m})\approx {\sqrt {\frac {2}{N}}}\left({\frac {1}{2}}F(0)+{\frac {(-1)^{m}}{2}}F(x_{N})+\sum _{n=1}^{N}~\cos(x_{m}x_{n})~F(x_{n})~\right)$ $\displaystyle \approx {\sqrt {\frac {2}{N}}}\left({\frac {1}{2}}F_{0}+{\frac {(-1)^{m}}{2}}F_{N}+\sum _{n=1}^{N-1}~\cos \left({\frac {\pi }{N}}mn\right)~F_{n}\right)$

At the transformation, it is assumed, that $F(x)$ can be neglected as $x\geq x_{N}$. In such a way,

$\!\!\!\!\!\!\!\!\!\!\!(6)~~~\displaystyle G_{m}\approx {\sqrt {\frac {2}{N}}}~(\mathrm {DCTI} _{N}~F)_{m}=(\Phi _{C1,N}~F)_{m}$

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 $F(x)=\exp(-x^{2}/2)$ 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 $x_{n}$, 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 $F$ be even periodic function of real argument with period $2\pi$:

$F(x)=F(-x)~$, $~F(x+2\pi )=F(x)$

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

$\displaystyle F(x)=\sum _{n=0}^{\infty }~a_{n}~\cos(nx)$

The Fourier–coefficients

$\displaystyle a_{0}={\frac {1}{\pi }}~\int _{0}^{\pi }~f(x)~\mathrm {d} x={\frac {1}{2\pi }}~\int _{-\pi }^{\pi }~f(x)~\mathrm {d} x$

$\displaystyle a_{n}={\frac {2}{\pi }}~\int _{0}^{\pi }~f(x)~\cos(nx)~\mathrm {d} x={\frac {1}{\pi }}~\int _{-\pi }^{\pi }~f(x)~\cos(nx)~\mathrm {d} x~$, $~n\!>\!0$

Assume some large natural number $N$. Let $\displaystyle x_{n}={\frac {\pi }{N}}n$. For approximation of coefficeins $a$, replace the integral with the finite sum:

$\displaystyle a_{m}\approx {\frac {1}{N}}\left({\frac {1}{2}}F(x_{0})+{\frac {1}{2}}F(x_{N})+\sum _{n=1}^{N-1}F(x_{n})\right)~$

$\displaystyle a_{m}\approx {\frac {2}{N}}\left({\frac {1}{2}}F(x_{0})+{\frac {1}{2}}F(x_{N})\cos(\pi m)+\sum _{n=1}^{N-1}F(x_{n})\cos(nx_{n})\right)~$, $~~n\!>\!0$

Comparison to equation (1) gives

$\displaystyle a_{0}\approx {\frac {1}{N}}~(\mathrm {DCFI} ~F)_{0}~$

$\displaystyle a_{n}\approx {\frac {2}{N}}~(\mathrm {DCFI} ~F)_{n}~$, $~~n\!>\!0$

Given expansion coefficients $a$, the function at the mesh $0,{\frac {\pi }{N}},{\frac {2\pi }{N}},..\pi$ can be evaluated with $(\mathrm {DCFI} ~G)$, with $G_{0}=2a_{0}$ and $G_{n}=a_{n}$ for $n\!>\!0$.

Numerical test of the expansion to the Fourier series

Let $N=8$. Let $F(x)=1+.1\cos(x)+.01\cos(2x)$. The expansion coefficients are expected to be $a_{0}=1$, $a_{1}=0.1$, $a_{2}=0.01$; 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 $N/2=4$. 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

$\mathrm {DCTI} _{N}$ 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 $N\!+\!1$ 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 $N=2^{q}$ for some natural number $q$. 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

The content of this article is adopted from http://tori.ils.uec.ac.jp/TORI/index.php/DCTI