[vorbis-dev] mdct explanation
Segher Boessenkool
segher at wanadoo.nl
Tue Jan 2 07:03:55 PST 2001
...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.
More information about the Vorbis-dev
mailing list