dct

Arguments

A

a real or complex vector or real or complex array (vector, matrix or N-D array.

X

a real or complex array with same shape as A.

sign

an integer. with possible values 1 or -1. Select direct or inverse transform. The default value is -1 (direct transform).

selection

a vector containing index on A array dimensions. See the Description part for details.

dims

a vector of positive numbers with integer values, or a vector of positive integers. See the Description part for details.

Each element must be a divisor of the total number of elements of A.

The product of the elements must be less than the total number of elements of A.

incr

a vector of positive numbers with integer values, or a vector of positive integers. See the Description part for details.

incr must have the same number of elements than dims.

Each element must be a divisor of the total number of elements of A.

The incr elements must be in strictly increasing order.

option

a character string. with possible values "dct1", "dct2", "dct4" or "dct" for direct transform and "dct1", "dct3", "dct4" or "idct" for inverse transform. The default value is "dct" for direct transform and "idct" for inverse transform. See the Description part for details.

Description

Transform description

This function realizes direct or inverse 1-D or N-D Discrete Cosine Transforms with shift depending on the option parameter value. For a 1-D array A of length n:

For "dct1" the function computes the unnormalized DCT-I transform:

$X(k) = A(1) + (-1)^{k-1}A(n) + 2\sum_{i=2}^{n-1} {A(i) \cos\frac{\pi (i -1)(k-1)}{n-1}}, k=1\ldots n$
For "dct2" the function computes the unnormalized DCT-II transform:

$X(k) = 2 \sum_{i=1}^{n} {A(i) \cos\frac{\pi (i -1/2)(k-1)}{n}}, k=1\ldots n$
For "dct3" the function computes the unnormalized DCT-III transform:

$X(k) = A(1) + 2 \sum_{i=2}^{n} {A(i) \cos\frac{\pi (i -1)(k-1/2)}{n}}, k=1\ldots n$
For "dct4" the function computes the unnormalized DCT-IV transform:

$X(k) = 2 \sum_{i=1}^{n} {A(i) \cos\frac{\pi (i -1/2)(k-1/2)}{n}}, k=1\ldots n$
For "dct" the function computes the normalized DCT-II transform:

$X(k) = \omega(k) \sum_{i=1}^n {A(i) \cos\frac{\pi (i -1/2)(k-1)}{n}}, k=1\ldots n \quad\text{with } \omega(1)=\frac{1}{\sqrt{n}} \quad\text{and } \omega(k)=\sqrt{\frac{2}{n}} , k>1$
For "idct" the function computes the normalized DCT-III transform:

$X(i) = \sum_{k=1}^n {\omega(k) A(k) \cos\frac{\pi (i -1/2)(k-1)}{n}}, k=1\ldots n \quad\text{with } \omega(1)=\frac{1}{\sqrt{n}} \quad\text{and } \omega(k)=\sqrt{\frac{2}{n}} , k>1$

The multi-dimensional DCT transforms , in general, are the separable product of the given 1d transform along each dimension of the array. For unnormalized transforms , computing the forward followed by the backward/inverse multi-dimensional transform will result in the original array scaled by the product of the dimension sizes.

The normalized multi-dimensional DCT transform of an array A with dimensions n₁, n₂, …, n_p is given by

$\begin{array} \\X(k_1, \dots, k_p) = \omega_1(k_1) \ldots \omega_p(k_p) \sum_{i_1=1}^{n_1} \ldots \sum_{i_p=1}^{n_p} {A(i_1,\ldots ,i_p) \cos\frac{\pi (2 i_1 -1)(k_1-1)}{2 n} \ldots \cos\frac{\pi (2 i_p -1)(k_p-1)}{2 n}}, k_j=1\ldots n_j\\ \text{with}\\ \omega_j(1)=\frac{1}{\sqrt{n_j}}\\ \omega_j(k)=\sqrt{\frac{2}{n_j}} , k>1 \end{array}$

The normalized multi-dimensional DCT inverse transform of an array A with dimensions n₁, n₂, …, n_p is given by

$\begin{array} \\X(i_1, \dots, i_p) = \sum_{k_1=1}^{n_1} \ldots \sum_{k_p=1}^{n_p} {\omega_1(i_1) \ldots \omega_p(i_p) A(k_1,\ldots ,k_p) \cos\frac{\pi (2 k_1 -1)(i_1-1)}{2 n} \ldots \cos\frac{\pi (2 k_p -1)(i_p-1)}{2 n}}, i_j=1\ldots n_j\\ \text{with}\\ \omega_j(1)=\frac{1}{\sqrt{n_j}}\\ \omega_j(k)=\sqrt{\frac{2}{n_j}} , k>1 \end{array}$

Syntax description

Short syntax

direct

X=dct(A,-1 [,option]) or X=dct(A [,option]) gives a direct transform according to the option value. The default is normalized DCT-II direct transform.

If A is a vector (only one dimension greater than 1) a 1-d transform is performed and in the other cases a n-dimensional transform is done.

(the -1 argument refers to the sign of the exponent..., NOT to "inverse"),

inverse

X=dct(A,1 [,option]) or X=idct(A [,option])performs the inverse transform.

If A is a vector (only one dimension greater than 1) a 1-d transform is performed and in the other cases a n-dimensional transform is done.

Long syntax for DCT along specified dimensions

X=dct(A,sign,selection [,option]) allows to perform efficiently all direct or inverse dct of the "slices" of A along selected dimensions.

For example, if A is a 3-D array X=dct(A,-1,2) is equivalent to:
```
for i1=1:size(A,1),
  for i3=1:size(A,3),
    X(i1,:,i3)=dct(A(i1,:,i3),-1);
  end
end
```
and X=dct(A,-1,[1 3]) is equivalent to:
```
for i2=1:size(A,2),
  X(:,i2,:)=dct(A(:,i2,:),-1);
end
```
X=dct(A,sign,dims,incr) is an old syntax that also allows to perform all direct or inverse dct of the slices of A along selected dimensions.

For example, if A is an array with n1*n2*n3 elements X=dct(A,-1,n1,1) is equivalent to X=dct(matrix(A,[n1,n2,n3]),-1,1). and X=dct(A,-1,[n1 n3],[1 n1*n2]) is equivalent to X=dct(matrix(A,[n1,n2,n3]),-1,[1,3]).

Optimizing dct

Remark: function automatically stores his last parameters in memory to re-use it in a second time. This improves greatly the time computation when consecutives calls (with same parameters) are performed.

It is possible to go further in dct optimization using get_fftw_wisdom, set_fftw_wisdom functions.

Algorithms

This function uses the fftw3 library.

Examples

1-D dct

//Frequency components of a signal
//----------------------------------
// build a sampled at 1000hz  containing pure frequencies
// at 50 and 70 Hz
sample_rate = 1000;
t = 0:1/sample_rate:0.6;
N = size(t,'*'); //number of samples
s = sin(2*%pi*50*t) + sin(2*%pi*70*t+%pi/4) + grand(1,N,'nor',0,1);
d = dct(s);
// zero low energy components
d(abs(d)<1) = 0;
size(find(d<>0), '*') //only 30 non zero coefficients out of 600
clf
plot(s,'b')
plot(dct(d,1),'r')

2-D dct

function z=__milk_drop(x, y)
    sq = x.^2+y.^2;
    z = exp( exp(-sq).*(exp(cos(sq).^20)+8*sin(sq).^20+2*sin(2*(sq)).^8) );
endfunction
x = -2:0.1:2;
[X,Y] = ndgrid(x,x);
A = __milk_drop(X,Y);
d = dct(A);
d(abs(d)<1)=0;
size(find(d<>0),'*')
A1 = dct(d,1);
clf
gcf().color_map = graycolormap(128);
subplot(121), grayplot(x,x,A)
subplot(122), grayplot(x,x,A1)

Bibliography

Matteo Frigo and Steven G. Johnson, "FFTW Documentation" http://www.fftw.org/#documentation

Report an issue
<< transforms	transforms	dst >>

idct

Syntax