Working with the DCT

Question

I am having a very hard time to implement the DCT algorithm. I have quite a few requirements like it has to work with NxN matrix or at least power of 2, it has to be 2D, it has to produce same output as FFTW fftwf_plan_r2r_2d(FFTW_REDFT10) it has to use real data and I need DCT I, II and III. And it has to be fast!

1) I could use FFT to compute DCT and wikipedia mentions this:

"One can also compute DCTs via FFTs combined with O(N) pre- and post-processing steps. In general, O(N log N) methods to compute DCTs are known as fast cosine transform (FCT) algorithms."

What are those 2 steps?

2) If I use Apples vDSP library is that even a good idea? since it needs an array twice as big as the DCT array (2N with radix2). And also vDSP is 1D only so I would end up doing it for each row and column to get a 2D result.

3) Is it possible to use an algorithm that works on 8x8 blocks and adapt it for my needs?

I am very new to learning about DSP so any help is appreciated! Thanks!

ps: Does anyone have any sample code for what I need? would help me a lot

Answer 1

You want to implement a 2D NxN Fast DCT-i (Discrete Cosine Transform) where -i refers to the type of DCT with slight modifications and type-II (DCT-II) being the most widely used one in such as old JPEG image codecs.

Note: since both 2D-FTT and 2D-DCT has a separable kernel, a 2D DCT-II is implemented from a 1D DCT-II in practice (as you described)

In the literature you may find many fast implementations of DCT. Here I will put an example how to use FFTs to implement a DCT-II in 1D (from K.R.Rao). It's not the most efficient one though.

step-1: given the N-point (length of N) signal xN[n] for n=0 to n=N-1
step-2: padd it with N-point zeros at the end to make it a 2N-point signal x2N[n]
step-3: compute 2N-point FFT of this 2N-point x2N[n] call it X2N[k]
step-4: multiply X2N[k] by c(k) * exp(-jkPI/2*N), where c(0)=1/sqrt(2) and c(k)=1 for k=1,...,N-1
step-5: naming the output of step4 as Y[k] then 1D DCT[k]= Re{Y[k]}, k=0,N-1
step 6: multiply DCT[k] with sqrt(2/N) if you need to get a Normalized DCT[k]

This is for 1D DCT-II computation. After getting it you should adopt it into row-column decomposition to get a 2D result.

here is a sample matlab code to verify the steps:

------------------------------------------------
N=8;               % Length of DCT-II                   
x=randn(1,N);      % create a test signal

c = ones(1,N);     % prepare the coefficients c[k]       
c(1)= 1/sqrt(2);   %      

x2 = [x zeros(1,N)];    % padd xN[n] with zeros  (step-2)   
X2 = fft(x2 , 2*N);     % get FFT of padded signal x2N[n]  (step-3)     

Y = c.* exp(-j* pi* [0:N-1]/(2*N)).* X2(1:N);   % (step:4)  
XdctNormalized = sqrt(2/N).* real(Y);    % (steps 5-6)   
-----------------------------------------------

You can actually find these in the following books:
1- 2D Image and Signal Processing (Jae S. Lim)
2- Techniques & Standards for Image-Video & Audio Coding (K.R.Rao)

Answer 2

While this is not much of an answer, here's a few comments while you wait for an expert to reply.

If this is one of your first DSP projects, it might be a bit ambitious just to dive straight into the fast implementations - you might want to start out with a naive implementation and then work your way up to the faster algorithms.

For the FFT-based approaches, you could check out how scipy does it: http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.fftpack.dct.html

There's a lot of papers on the subject, but they might be a bit tricky to chew through. A couple of examples through some brief searches:

• http://www.cin.ufpe.br/~vak/tg/papers/Chen-dct.pdf (1977 paper on a fast DCT)
• http://arxiv.org/pdf/1203.3442.pdf (fast DCT-II)

It would probably be easier to understand from a book - I recall reading about the DCT and some fast algorithms with pretty elaborate explanations in "JPEG Still image data compression standard" by W. B. Pennebaker and J. L. Mitchell (I'm sure you can find loads of books on the DCT, specifically, though).

Also, this seems to be how they do it over at FFmpeg: http://www.ffmpeg.org/doxygen/2.1/dct_8c_source.html