背景

  DCT是离散余弦变换的缩写，由于其变换后具有较高的能量聚集度，通常作为音视频编码的变换去使用。而由于DCT的块效应，人们发明了很多方法去克服块效应。例如LOT、MDCT。在aac的编码中采用时域重叠的MDCT去实现(TDAC)。本博文仅从DFT到DCT的推导以及MDCT的编解码流程进行讲解,力求以数学的推导来阐明过程。

DFT : 离散傅立叶变换. 用于将离散的时域信号转换到频域上。
DCT : 离散余弦变换，也是正交变换。用于将离散的时域信号转换为频域上的信息
MDCT : 改进后的离散余弦变换. 通过时域重叠来消除混叠。
IMDCT : MDCT的逆变换，时域信号在经过MDCT编码以及IMDCT解码后，还原出的并不是原始信号

DFT到DCT的推导

DFT : $\large X(k) = \sum_{n=0}^{N-1} x[n].e^{\frac{-j.2\pi.kn}{N}} \text{ }$
欧拉公式 : $\large e^{-j\theta} = cos\theta + j.sin\theta$

step：

1. 虚部为0: 观察DFT变换可得，当其为实偶信号时，虚部为0。因为实偶信号的性质是x(n) = - x(n)，故在将DFT的复数部分拆开后由于其虚部为奇函数，故实偶信号的虚部将会抵消。
2. 构建实偶信号: 时域信号经抽样后皆为实数，而要满足偶函数的性质需要人为构造。
   假设抽样后具有从0到N-1的N点离散数字信号，其数学定义为 $\large x[m] = \{ {x[0],....,x[N-1]} \}$。将该序列进行偶延拓，其数学定义更改为
   $\acute{x[m]} = \begin{cases} x[m], & \text{if n belong to \{ {0,..,N-1} \}} \\ x[-m-1], & \text{if n belong to \{ {-N,..,-1} \} } \end{cases}$
   $\acute{x[m]}$信号如下图1所示：
   
   再将 $\acute{x[m]}$序列整体向右偏移 $\Large\frac{1}{2}$，令 $\large\ddot{x[m]}$为 $\large\acute{x[m-\frac{1}{2}]}$， $\large\ddot{x[m]}$如下图2所示:
   
3. 重新推导实偶信号的DFT公式: $\large X(k) = \sum_{m=-N+\frac{1}{2}}^{N-\frac{1}{2}} \ddot{x[m - \frac{1}{2}]}.e^{\frac{-j.2\pi.km}{2N}} \text{ } = 2 *\sum_{m=\frac{1}{2}}^{N-\frac{1}{2}} \ddot{x[m - \frac{1}{2}]}.e^{\frac{-j.2\pi.km}{2N}}$令n = m + $\frac{1}{2}$,则上式可化为 $\Large2*\sum_{n=0}^{N-1} \acute{x[n]}.cos(\frac{(n+\frac{1}{2}).k\pi}{N})$
4. 正交变换: 将DCT变换中与x[n]相乘的系数组织成矩阵C,如果该矩阵正交则有 $C.C^{T} = E$.故将变换核的系数2做变换可得下式: $\Large\sqrt{\frac{2}{N}}.g_k*\sum_{n=0}^{N-1} \acute{x[n]}.cos(\frac{(n+\frac{1}{2}).k\pi}{N}) \tag{1}$
   其中 $g_k$的数学定义为:
   $\large g_k = \begin{cases} 1/\sqrt{2}, & \text{ k == 0} \\ 1, & \text{ k != 0 } \end{cases}$

MDCT的编解码流程简述

  MDCT作为改进的离散余弦变换，所以编码由DCT过渡到MDCT是其本身的优势的。DCT在二维图片分量的变换中，其变换系数的高频分量集中在左上角(转换矩阵的左上角)，而由于图片的编码是将整体图片切割成一个个小方块进行编码转换，更是造成了相邻方块间在转换之后容易引入噪声，这就是方块效应，在视觉上表示为图片编码后相邻小方块间的白条。
  而诸如LOT、MDCT采用了TDAC实现的编码转换，转换后的单位抽样响应是由中间向其两边递减的，如下图3所示:

故MDCT可以很好的消除方块效应。
  在MDCT变换中，输入的离散数字信号长度为2N，但是经过IMDCT[MDCT[x[n]]]的有效信号长度实则为N，下图4能很好的表示出来:

现对上图4的编解码流程进行数学推导：

1. MDCT变换公式：
   $\Large X(k) = \frac{2}{N}*\sum_{n=0}^{N-1} x[n].cos[ \frac{2\pi}{N}.(n+\frac{1}{2}+\frac{N}{4}).(k+\frac{1}{2})] \quad \text{k $\in \{0,..,N/2-1\} $}$
   在MDCT变换中，由于X(k) == X(N+k)，所以X(k)只有N/2个独立分量，故k的范围为 $\{ 0,..,N/2-1\}$
2. IMDCT变换公式:
   $\Large x(n) = 2*\sum_{k=0}^{\frac{N}{2}-1} X[k].cos[ \frac{2\pi}{N}.(n+\frac{1}{2}+\frac{N}{4}).(k+\frac{1}{2})] \quad \text{n $\in \{0,..,N-1\} $}$
3. 如何从解码端获取原始信号:
     假设输入信号的序列为 $\large x[n]=\{ x_1,x_2 \}$,现证明经过MDCT变换以及IMDCT变换后的输出信号 $\large y[n]=\{ x_1-\acute{x_1},x_2+\acute{x_2} \}$, $\large\acute{x_1}$为 $\large x_1$的逆序序列,而 $\large\acute{x_2}$为 $\large x_2$的逆序序列。
     令输入的离散信号长度N为4， $\large x[n]=\{x_0,x_1,x_2,x_3\}$,则需证明 $\large y[n] =\{ x_0 - x_1,x_1-x_0,x_2-x_3,x_3-x_2\}$
     令长度N为4的MDCT变换矩阵为C，则C的数学定义如下:
   $C_k,_n= \begin{bmatrix} C_0,_0 & C_1,_0 \\ C_0,_1 & C_1,_1 \\ C_0,_2 & C_1,_2 \\ C_0,_3 & C_1,_3 \\ \end{bmatrix} ,\quad y[n]=x[n].C.C^T \Longrightarrow \quad y[n]=x[n].(CC^T)$
     再令 $Q = CC^T$，且Q为4*4矩阵，则上述证明转换为推导 $Q_0,_0 = 1,Q_0,_1=-1$。再N=4的情况下，C表示如下:
   $C_k,_n= \begin{bmatrix} cos\frac{3}{8}\pi & cos\frac{9}{8}\pi \\ cos\frac{5}{8}\pi & cos\frac{15}{8}\pi \\ cos\frac{7}{8}\pi & cos\frac{21}{8}\pi \\ cos\frac{9}{8}\pi & cos\frac{27}{8}\pi \\ \end{bmatrix}\quad,\quad cosa.cosb = \frac{cos(a+b) + cos(a-b)}{2}$
      $Q_0,_0=C_0,_0*C_0,_0 + C_1,_0*C_1,_0 \Longrightarrow cos\frac{3}{8}\pi.cos\frac{3}{8}\pi +cos\frac{9}{8}\pi.cos\frac{9}{8}\pi \Longrightarrow 1$
      $Q_0,_1=C_0,_1*C_0,_0 + C_1,_1*C_1,_0 \Longrightarrow cos\frac{5}{8}\pi.cos\frac{3}{8}\pi +cos\frac{15}{8}\pi.cos\frac{9}{8}\pi \Longrightarrow -1$
      $Q_0,_2=C_0,_2*C_0,_0 + C_1,_2*C_1,_0 \Longrightarrow cos\frac{7}{8}\pi.cos\frac{3}{8}\pi +cos\frac{21}{8}\pi.cos\frac{9}{8}\pi \Longrightarrow 0$
      $Q_0,_3=C_0,_3*C_0,_0 + C_1,_3*C_1,_0 \Longrightarrow cos\frac{9}{8}\pi.cos\frac{3}{8}\pi +cos\frac{27}{8}\pi.cos\frac{9}{8}\pi \Longrightarrow 0$

   故
$\large \begin{bmatrix} x_0 & x_1 &x_2 &x_3 \\ \end{bmatrix} * \begin{bmatrix} 1 & C_1,_0 &C_2,_0 &C_3,_0 \\ -1 &C_1,_1 &C_2,_1 &C_3,_1 \\ 0 &C_1,_2 &C_2,_2 &C_3,_2 \\ 0 &C_1,_3 &C_2,_3 &C_3,_3 \\ \end{bmatrix} = \{y_0,y_1,y_2,y_3\}$
   可得 $\large y_0 = x_0 - x_1$,后续 $\large y_1$的推导读者可以自证。
   令 $\breve{x_i} =\{ x_i,x_{i+1} \}$, $\breve{x_{i+1}}=\{x_{i+1},x_{i+2}\}$,在MDCT的输入序列中，当前序列和下个序列的时域重叠为50%.而 $y_i= IMDCT(MDCT(\breve{x_i})) = \{\ x_i - \acute{x_i},x_{i+1} + \acute{x_{i+1}} \}$ $y_{i+1}= IMDCT(MDCT(\breve{x_{i+1}})) = \{\ x_{i+1} - \acute{x_{i+1}},x_{i+2} + \acute{x_{i+2}} \}$
   再令输出序列的 $y_i,y_{i+1}$进行时域50%重叠，即可还原出2 $\large x_{i+1}$

MDCT的快速算法

N点的MDCT可以转化为N/2点的DCT-IV进行计算。

1. MDCT的公式如下:
   $\Large X(k) = \frac{2}{N}*\sum_{n=0}^{N-1} x[n].cos[ \frac{2\pi}{N}.(n+\frac{1}{2}+\frac{N}{4}).(k+\frac{1}{2})] \quad \text{k $\in \{0,..,N/2-1\} $}$
   令N=2m，可变换为:
   $\Large X(k) = \frac{1}{M}*\sum_{n=0}^{2M-1} x[n].cos[ \frac{pi}{M}.(n+\frac{1}{2}+\frac{M}{2}).(k+\frac{1}{2})] \quad \text{k $\in \{0,..,N/2-1\} $}$
   同时令 $\varepsilon_n=cos[\frac{\pi}{M}.(n+\frac{1}{2}).(k+\frac{1}{2})]$,而 $\varepsilon_n$即是DCT变换核的系数因子。现通过 $\varepsilon_n$将MDCT化简:
   $\Large X(k) = \frac{1}{M}*\sum_{n=0}^{2M-1} x[n].cos[ \frac{2\pi}{N}.(n+\frac{1}{2}+\frac{M}{2}).(k+\frac{1}{2})] = \sum_{n=0}^{2M-1} x_n.\varepsilon_{\frac{M}{2}+n}$
   $\Large = \sum_{n=0}^{\frac{M}{2}-1} [ x_n.\varepsilon_{\frac{M}{2}+n} + x_{M-1-n}.\varepsilon_{\frac{3M}{2}-1-n}+ \varepsilon_{2M-1-n}.x_{\frac{3M}{2}-1-n} + \varepsilon_{2M+n}.x_{\frac{3M}{2}+n} ]\tag{3}$又因为 $\large\varepsilon_{2M+n}=\varepsilon_{2M-1-n}=-\varepsilon_n$，式3可化为:
   $\Large \quad \sum_{n=0}^{\frac{M}{2}-1} [ x_n.\varepsilon_{\frac{M}{2}+n} - x_{M-1-n}.\varepsilon_{\frac{M}{2}+n} - \varepsilon_{n}.x_{\frac{3M}{2}-1-n} - \varepsilon_{2M+n}.x_{\frac{3M}{2}+n} ]$ $\Large = \sum_{n=M/2}^{M-1} \varepsilon_{n}.(x_{n-\frac{M}{2}} - x_{\frac{3M}{2}-1-n} ) - \sum_{n=0}^{\frac{M}{2}-1} \varepsilon_{n}.( x_{\frac{3M}{2}-1-n} + x_{\frac{3M}{2}+n} )$则MDCT可以化为类似DCT变换的求和式: $\large\sum_{n=0}^{M-1}\breve{x_n}\varepsilon_n,\breve{x_n}$为:
   $\large \breve{x_n} = \begin{cases} -(x_{\frac{3M}{2}-1-n}+x_{\frac{3M}{2}+n}), & \text{ n = \{0,..,$\frac{M}{2}$-1 \}} \\ (x_{n-\frac{M}{2}} - x_{\frac{3M}{2}-1-n}),, & \text{ n = \{$\frac{M}{2}$,..,$\small M$ -1 \}} \end{cases}$ 所以2N点的MDCT可以化为N/2的DCT进行计算，而DCT可以转化为FFT通过蝶形单元减少算法的计算时间度。