How does this 1D IDCT work?

Question

I have an implementation of the inverse discrete cosine transform, and I'm trying to figure out how they got to this code. So far I have found out that this is probably the optimized implementation of the Cooley-Tukey radix-2 Decimation-in-time for DCT instead of the DFT (Discrete Fourier Transform).

However, I still do not understand what is happening at each stage. I figured Wx constants were probably turning factors.

Can someone give a link to the explanation or give some explanation to this code?

//Twiddle factors #define W1 2841 /* 2048*sqrt(2)*cos(1*pi/16) */ #define W2 2676 /* 2048*sqrt(2)*cos(2*pi/16) */ #define W3 2408 /* 2048*sqrt(2)*cos(3*pi/16) */ #define W5 1609 /* 2048*sqrt(2)*cos(5*pi/16) */ #define W6 1108 /* 2048*sqrt(2)*cos(6*pi/16) */ #define W7 565 /* 2048*sqrt(2)*cos(7*pi/16) */ //Discrete Cosine Transform on a row of 8 DCT coefficients. NJ_INLINE void njRowIDCT(int* blk) { int x0, x1, x2, x3, x4, x5, x6, x7, x8; int t; if (!((x1 = blk[4] << 11) | (x2 = blk[6]) | (x3 = blk[2]) | (x4 = blk[1]) | (x5 = blk[7]) | (x6 = blk[5]) | (x7 = blk[3]))) { blk[0] = blk[1] = blk[2] = blk[3] = blk[4] = blk[5] = blk[6] = blk[7] = blk[0] << 3; return; } x0 = (blk[0] << 11) + 128; //For rounding at fourth stage //First stage /*What exactly are we doing here? Do the x values have a meaning?*/ x8 = W7 * (x4 + x5); x4 = x8 + (W1 - W7) * x4; x5 = x8 - (W1 + W7) * x5; x8 = W3 * (x6 + x7); x6 = x8 - (W3 - W5) * x6; x7 = x8 - (W3 + W5) * x7; //Second stage x8 = x0 + x1; x0 -= x1; x1 = W6 * (x3 + x2); x2 = x1 - (W2 + W6) * x2; x3 = x1 + (W2 - W6) * x3; x1 = x4 + x6; x4 -= x6; x6 = x5 + x7; x5 -= x7; //Third stage x7 = x8 + x3; x8 -= x3; x3 = x0 + x2; x0 -= x2; x2 = (181 * (x4 + x5) + 128) >> 8; x4 = (181 * (x4 - x5) + 128) >> 8; //Fourth stage blk[0] = (x7 + x1) >> 8; //bit shift is to emulate 8 bit fixed point precision blk[1] = (x3 + x2) >> 8; blk[2] = (x0 + x4) >> 8; blk[3] = (x8 + x6) >> 8; blk[4] = (x8 - x6) >> 8; blk[5] = (x0 - x4) >> 8; blk[6] = (x3 - x2) >> 8; blk[7] = (x7 - x1) >> 8; } 

Answer 1

I am not an expert in DCT, but I wrote several FFT implementations at one time, so I'm going to answer that. Please do the following with a pinch of salt.

 void njRowIDCT(int* blk) 

You correctly say that the algorithm seems to be 8-bit DCX Radix-2, which uses fixed-point arithmetic with an accuracy of 24: 8. I guess the accuracy, because the last step is shifted by 8 to get the desired one (this is also a story about a fairy tale; )

Since it is 8-length, its power is 3 (2 ^ 3 = 8), which means there are 3 steps in DCT. So far, all this is very similar to the FFT. The "fourth stage", apparently, is simply scaling to restore the original accuracy after fixed-point arithmetic.

As far as I can see, the input frame is a bit-reversal from the input blk array to the local variables x0-x7. x8 seems to be a temporary variable. Sorry, I can’t be more descriptive than that.

Bit reversal stage

Update

Take a look at DSP for scientists and engineers . It provides a clear and accurate explanation of signal processing topics. This chapter is devoted to DCT (please go to p497).

Wn (twiddle coefficients) correspond to the basic functions in this chapter, although note that this is a DCT 8x8 (2D) description.

As for the 3 steps I mentioned, compare with the description of the 8-point FFT:

FFT executes butterflies on a bit-reversible input matrix (which are essentially complex multiple additions), multiplying one path by the factor Wn or twiddle on this path. FFT is performed in stages. I still don't understand what your DCT code does, but it can help decompose it into a diagram.

This or someone who knows what they are talking about the promotion ,-)

I will re-read this page and edit when I decrypt the code.

Answer 2

Both DFT and DCT are simply linear transformations that can be represented as a one-time complex matrix (sometimes shortened for strictly real input). Thus, you can simply combine the above equations to obtain a formula for each final term, which should be equivalent to one line of the linear transformation matrix (ignoring rounding problems). Then, see how the above code sequence manually performs general optimization of subexpression or refactoring between calculations and / or within lines.