DCTI: Difference between revisions
imported>Dmitrii Kouznetsov m (add "references" section) |
imported>Dmitrii Kouznetsov (spacing..) |
||
| Line 11: | Line 11: | ||
DCTI, or Discrete Cos Transform of kind I, is orthogonal transform, that repalces an array <math>F</math> of length <math>N\!+\!1</math> with elements <math>F_n~</math>, <math>~n=0 .. N</math> to the array <math>G</math> with elements | DCTI, or Discrete Cos Transform of kind I, is orthogonal transform, that repalces an array <math>F</math> of length <math>N\!+\!1</math> with elements <math>F_n~</math>, <math>~n=0 .. N</math> to the array <math>G</math> with elements | ||
: <math> \!\!\!\!\!\!\!\!\!\!(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} n k \right)~</math> for <math>~k=0, .. N</math> | : <math> \!\!\!\!\!\!\!\!\!\!(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} n k \right) ~ ~ </math> for <math>~ ~ k=0, .. N</math> | ||
==Normalized form== | ==Normalized form== | ||
Revision as of 07:19, 2 September 2012
DCTI is one of realizations of the Discrete Cos transform operator. The name is created in analogy with DCT by 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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle F} of length Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N\!+\!1} with elements Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle F_n~} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle ~n=0 .. N} to the array Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle G} with elements
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \!\!\!\!\!\!\!\!\!\!(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} n k \right) ~ ~ } for Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle ~ ~ k=0, .. N}
Normalized form
The orthonormaized transform can be defined with operator Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Phi_{\mathrm C1,N}} , that acts on array Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle F} in the following way:
Operator Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \Phi_{\mathrm C1 , N}} is its own inverse; Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle (\Phi_{\mathrm C1,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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle F} of length Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N\!+\!1} has form zcosft1(F-1,N); after such a call, values of the elements of array Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle F} are replaced with values calculated with the expression (1) above. For Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N=2^q} , the evaluation requires of order of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle Nq} operations
Approximation of the CosFourier
The CosFourier operator transforms a function Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle F} of non–negative argument to function Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle G} in the following way:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \!\!\!\!\!\!\!\!\!\!\! (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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N} . Let Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x_n=\sqrt{\pi/N}~ n} . Let function Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle F} be smooth and quickly decay at infinity. Then, the transform of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle F} can be approximated as follows:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \!\!\!\!\!\!\!\!\!\!\! (4) ~ ~ ~ \displaystyle G(x) = \mathrm{CosFourier} F(x) \approx \sqrt{\frac{2}{\pi}} \left ( \frac{1}{2} F(0) + \sum_{n=1}^{N} ~\cos(x x_n)~ F(x_n) \right)~ \sqrt{{\pi}/N} }
For Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x=x_m} , this can be written as follows:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \!\!\!\!\!\!\!\!\!\!\! (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) } Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle F(x)} can be neglected as Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x\ge x_N} . In such a way,
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \!\!\!\!\!\!\!\!\!\!\! (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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle F} be even periodic function of real argument with period Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 2\pi} :
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle F(x)=F(-x)~} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle ~F(x+2\pi)=F(x)}
Then, the expansion into the Fourier series can be written as
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \displaystyle F(x)=\sum_{n=0}^\infty ~a_n~ \cos(n x)}
The Fourier–coefficients
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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}
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \displaystyle a_n=\frac{2}{\pi} ~ \int_0^\pi ~ f(x) ~ \cos(n x)~ \mathrm d x= \frac{1}{\pi} ~ \int_{-\pi}^\pi ~ f(x) ~ \cos(n x)~ \mathrm d x~} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle ~n\!>\!0}
Assume some large natural number Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N} . Let Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \displaystyle x_n=\frac{\pi}{N} n} . For approximation of coefficeins Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a} , replace the integral with the finite sum:
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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)~}
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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(n x_n) \right)~} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle ~ ~ n\!>\!0}
Comparison to equation (1) gives
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \displaystyle a_0 \approx \frac{1}{N} ~(\mathrm {DCFI}~F)_0 ~}
- Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \displaystyle a_n \approx \frac{2}{N} ~(\mathrm {DCFI}~F)_n ~} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle ~ ~ n\!>\!0}
Given expansion coefficients Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle a} , 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
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
- ↑ http://en.wikipedia.org/wiki/Discrete_cosine_transform
- ↑ 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.
The content of this article is adopted from http://tori.ils.uec.ac.jp/TORI/index.php/DCTI