Multidimensional Discrete Cosine Transform (DCT)

Fast forward and inverse Multidimensional Discrete Cosine Transforms (DCT, IDCT).

Description

The function is much faster than Matlab's native (dct, idct, dct2, idct2). It also allows N-D (multidimensional) input. The function takes advantage of a) FFT b) fast permutation through indicies c) persistent precomputation.

Example:
x=randn(100,200,300);

y=mirt_dctn(x); % forward DCT
x=mirt_idctn(y); % inverse DCT

Find more at my home page:
http://www.bme.ogi.edu/~myron

Enjoy and comment below if you like it!

• MATLAB release: MATLAB 7.7 (R2008b)
• Other Files: license.txt, ; mirt_dctn.m, ; mirt_idctn.m

[Mike Berg]

Great job. Thanks.

[Paul Matthews]

 

[Paul Matthews]

It doesnt seem faster than dct to me.
In fact it is slower for row vectors-

tic; for i=1:3000;x=rand(1,512);dct(x);end; toc % 1.5
tic; for i=1:3000;x=rand(1,512);idct(x);end; toc % 1.5

tic; for i=1:3000;x=rand(1,512);mirt_dctn(x);end; toc % 3.0
tic; for i=1:3000;x=rand(1,512);mirt_idctn(x);end; toc % 2.7

But for 2d it seems about twice as fast as dct.

It is a mystery to me why Matlab have not written a fast dct

[Andriy Myronenko]

to Paul Matthews,
 thanks for the feedback,

I was surprised by your results at first, then I find out that your results in 1D are correct only for ROW vectors. Try the same with COLUMN vectors, please:

tic; for i=1:3000;x=rand(512,1);dct(x);end; toc % 0.81s
tic; for i=1:3000;x=rand(512,1);idct(x);end; toc % 0.91s

tic; for i=1:3000;x=rand(512,1);mirt_dctn(x);end; toc % 0.66s
tic; for i=1:3000;x=rand(512,1);mirt_idctn(x);end; toc % 0.70s

For the larger 1D vectors, mirt_dctn (and mirt_idctn) performance improvement is even bigger
 tic; for i=1:3000;x=rand(4096,1);mirt_dctn(x);end; toc % 2.1s
 tic; for i=1:3000;x=rand(4096,1);dct(x);end; toc % 4.55s.

PS: I'll take a look why my 'row' vector processing takes longer, and correct it shortly.

[Andriy Myronenko]

A 'row' vector bug solved and corrected.

[candelabro candelabro]

These functions can be changed to allows for zero padding?
It is an important concern in FFT...
I'm not sure if it can be done by simply changing in your codes the term "fft(x)" by "fft(x,n)", where n is the size of the resulting padded output...

[candelabro candelabro]

In tree dimensions how I can obtain the part corresponding to the negative frequencies, that is, the output can be centered with the fftshift fucntion?

[Andriy Myronenko]

To candelabro
1) As for the zero-padding, modifying fft(x,n) is incorrect.
Just pad the input 3D vector x before taking the dct.
For instance, use :

newsize=[M N K]; % Define padded size of the 3D array;
tmp=zeros([M N K]); % padded data buffer
cutsize=min([size(x); newsize]);
tmp(1:cutsize(1),1:cutsize(2),1:cutsize(3))=x(1:cutsize(1),1:cutsize(2),1:cutsize(3));
x=tmp;
y=mirt_dctn(x); take 3D DCT (padded or reduced to the M N K size)

2) The positive and negative frequencies is only appropriate for fft. If you really want to have the analogy with fft, then remember that DCT gives an output for positive frequencies. And negative frequencies can be obtained by mirror flip.