[vorbis-dev] mdct explanation

[Segher Boessenkool]

...as promised.

This describes the mdct used in my d.m.l patch. I think it is the
same as the Lee fast-dct.

I typed it in a kind of pseudo-TeX, 'cause the ascii art would
kill me. Hope you can read TeX source; if not, ask someone who
can to make a .ps/.gif/.whatever of the TeX output, and put it
on a webpage or something. I'm to lazy to do it (and besides, I
don't have access to TeX, X, or anything else from a text terminal
with gcc and vi right now).

Hope you like it,

Segher


--

Suppose we want to calculate a MDCT (encode side):

(x_ is input (time) samples, w_ is window, b_ is output (freq) samples):

b_k = \sum_{j=0}^{4N-1} x_j w_j \cos{(2k+1)(2j+1-2N) \over 8N} \pi

After setting

a_j := w_{-j-1+N} x{-j-1+N} + w_{j+N} x{j+N},   j < N
a_j := w_{j-1-N} x{j-1-N} - w_{j+3N} x{j+3N},   otherwise

(or something similar; the cos'es should agree)

we get

b_k = \sum_{j=0}^{2N-1} a_j \cos{(2k+1)(2j+1) \over 8N} \pi

This is a DCT-IV.

For decode side, we first compute DCT-IV, and then "wrap it out"
to length 4N, with window.

We'll compute

c_k := 2 \cos{2k+1 \over 8N} \pi b_k

o

c_k = \sum_{j=0}^{2N-1} a_j ( 
        \cos{(2k+1)(j) \over 4N} \pi + \cos{(2k+1)(j+1) \over 4N} \pi
)

c_k = \sum_{j=0}^{2N-1} a_j \cos{(2k+1)(j) \over 4N} \pi

+ \sum_{j=1}^{2N-1} a_{j-1} \cos{(2k+1)(j) \over 4N} \pi

Now set

a'_j := a_j + a_{j-1}    (but a'_0 = a_0)

then

c_k = \sum_{j=0}^{2N-1} a'_j \cos{(2k+1)(j) \over 4N} \pi

This is a DCT-II.

We can split it:

c_k = \sum_{j=0}^{N-1} a'_{2j} \cos{(2k+1)(j) \over 2N} \pi

+ \sum_{j=0}^{N-1} a'_{2j+1} \cos{(2k+1)(2j+1) \over 4N} \pi

(for k < N),

and

c_(2N-1-k} = \sum_{j=0}^{N-1} a'_{2j} \cos{(2k+1)(j) \over 2N} \pi

- \sum_{j=0}^{N-1} a'_{2j+1} \cos{(2k+1)(2j+1) \over 4N} \pi

(otherwise).

Thus we calculate (for k = 0..N-1)

d_k := \sum_{j=0}^{N-1} a'_{2j} \cos{(2k+1)(j) \over 2N} \pi

(this is a half length DCT-II),

and

e_k := \sum_{j=0}^{N-1} a'_{2j+1} \cos{(2k+1)(2j+1) \over 4N} \pi

(a half length DCT-IV).

Hey, we know how to calculate those! Just recursion... The
endcase (startcase) for the recursion is the trivial summation
(N=1).

So, for calculating a MDCT:

Calculate a_j from x_j and w_j; now we need to calculate a DCT-IV.

To calculate a DCT-IV:

calc a'_j, do a DCT-II on it to get c_k, calc

b_k := c_k / ( 2 \cos{2k+1 \over 8N} \pi )

done.

To calc a DCT-II:

calc d_k is the DCT-II of the even a'_j, and e_k is the DCT_IV of
the odd,

c_k = d_k + e_k
c_{2N-1-k} = d_k - e_k

"That's all there is, and there ain't no more".

Oh, this only works for N a power of two, but you can split the
DCT-II in three or more parts as well. It will need some mul's,
though. If you need to split on three, do it _last_; there's
some nice symmetries and zeroes in there.

--- >8 ----
List archives:  http://www.xiph.org/archives/
Ogg project homepage: http://www.xiph.org/ogg/
To unsubscribe from this list, send a message to 'vorbis-dev-request at xiph.org'
containing only the word 'unsubscribe' in the body.  No subject is needed.
Unsubscribe messages sent to the list will be ignored/filtered.