DCT和DFT的关系以及MDCT的推导证明

背景

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

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

DFT到DCT的推导

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

step：

1. 虚部为0: 观察DFT变换可得，当其为实偶信号时，虚部为0。因为实偶信号的性质是x(n) = - x(n)，故在将DFT的复数部分拆开后由于其虚部为奇函数，故实偶信号的虚部将会抵消。
2. 构建实偶信号: 时域信号经抽样后皆为实数，而要满足偶函数的性质需要人为构造。
   假设抽样后具有从0到N-1的N点离散数字信号，其数学定义为 x [ m ] = { x [ 0 ] , . . . . , x [ N − 1 ] } \large x[m] = \{ {x[0],....,x[N-1]} \} x[m]={x[0],....,x[N−1]}。将该序列进行偶延拓，其数学定义更改为
   $$\acute{x[m]}=\{\begin{array}{ll}x[m],& \text{if n belong to { 0,..,N-1 }}\\ x[-m-1],& \text{if n belong to { -N,..,-1 } }\end{array}$$
   $\acute{x[m]}$信号如下图1所示：
   
   再将$\acute{x[m]}$序列整体向右偏移$\frac{1}{2}$，令$\ddot{x[m]}$为$\acute{x[m-\frac{1}{2}]}$，$\ddot{x[m]}$如下图2所示:
   
3. 重新推导实偶信号的DFT公式: $$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}}=2\ast \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}$,则上式可化为$2\ast {\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做变换可得下式:$$\sqrt{\frac{2}{N}}.g_k\ast \sum _{n=0}^{N-1}\acute{x[n]}.cos(\frac{(n+\frac{1}{2}).k\pi }{N})\text{\hspace{1em}(1)}$$
   其中$g_k$的数学定义为:
   $$g_k=\{\begin{array}{ll}1/\sqrt{2},& \text{ k == 0}\\ 1,& \text{ k != 0 }\end{array}$$

MDCT的编解码流程简述

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

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

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

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

   故
[ x 0 x 1 x 2 x 3 ] ∗ [ 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 ] = { y 0 , y 1 , y 2 , y 3 } \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\} [x0​​x1​​x2​​x3​​]∗⎣⎢⎢⎢⎡​1−100​C1​,0​C1​,1​C1​,2​C1​,3​​C2​,0​C2​,1​C2​,2​C2​,3​​C3​,0​C3​,1​C3​,2​C3​,3​​⎦⎥⎥⎥⎤​={y0​,y1​,y2​,y3​}
   可得$y_0=x_0-x_1$,后续$y_1$的推导读者可以自证。
   令 x i ˘ = { x i , x i + 1 } \breve{x_i} =\{ x_i,x_{i+1} \} xi​˘​={xi​,xi+1​}, x i + 1 ˘ = { x i + 1 , x i + 2 } \breve{x_{i+1}}=\{x_{i+1},x_{i+2}\} xi+1​˘​={xi+1​,xi+2​},在MDCT的输入序列中，当前序列和下个序列的时域重叠为50%.而 y i = I M D C T ( M D C T ( x i ˘ ) ) = {   x i − x i ˊ , x i + 1 + x i + 1 ˊ } y_i= IMDCT(MDCT(\breve{x_i})) = \{\ x_i - \acute{x_i},x_{i+1} + \acute{x_{i+1}} \} yi​=IMDCT(MDCT(xi​˘​))={ xi​−xi​ˊ​,xi+1​+xi+1​ˊ​} y i + 1 = I M D C T ( M D C T ( x i + 1 ˘ ) ) = {   x i + 1 − x i + 1 ˊ , x i + 2 + x i + 2 ˊ } y_{i+1}= IMDCT(MDCT(\breve{x_{i+1}})) = \{\ x_{i+1} - \acute{x_{i+1}},x_{i+2} + \acute{x_{i+2}} \} yi+1​=IMDCT(MDCT(xi+1​˘​))={ xi+1​−xi+1​ˊ​,xi+2​+xi+2​ˊ​}
   再令输出序列的$y_i,y_{i+1}$进行时域50%重叠，即可还原出2$x_{i+1}$

MDCT的快速算法

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

1. MDCT的公式如下:
   $$X(k)=\frac{2}{N}\ast \sum _{n=0}^{N-1}x[n].cos[\frac{2\pi }{N}.(n+\frac{1}{2}+\frac{N}{4}).(k+\frac{1}{2})]\text{k ∈{0,..,N/2−1}}$$
   令N=2m，可变换为:
   $$X(k)=\frac{1}{M}\ast \sum _{n=0}^{2M-1}x[n].cos[\frac{\pi }{M}.(n+\frac{1}{2}+\frac{M}{2}).(k+\frac{1}{2})]\text{k ∈{0,..,N/2−1}}$$
   同时令${\epsilon }_n=cos[\frac{\pi }{M}.(n+\frac{1}{2}).(k+\frac{1}{2})]$,而${\epsilon }_n$即是DCT变换核的系数因子。现通过${\epsilon }_n$将MDCT化简:
   $$X(k)=\frac{1}{M}\ast \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.{\epsilon }_{\frac{M}{2}+n}$$
   $$=\sum _{n=0}^{\frac{M}{2}-1}[x_n.{\epsilon }_{\frac{M}{2}+n}+x_{M-1-n}.{\epsilon }_{\frac{3M}{2}-1-n}+{\epsilon }_{2M-1-n}.x_{\frac{3M}{2}-1-n}+{\epsilon }_{2M+n}.x_{\frac{3M}{2}+n}]\text{\hspace{1em}(3)}$$又因为${\epsilon }_{2M+n}={\epsilon }_{2M-1-n}=-{\epsilon }_n$，式3可化为:
   $$\sum _{n=0}^{\frac{M}{2}-1}[x_n.{\epsilon }_{\frac{M}{2}+n}-x_{M-1-n}.{\epsilon }_{\frac{M}{2}+n}-{\epsilon }_n.x_{\frac{3M}{2}-1-n}-{\epsilon }_{2M+n}.x_{\frac{3M}{2}+n}]$$$$=\sum _{n=M/2}^{M-1}{\epsilon }_n.(x_{n-\frac{M}{2}}-x_{\frac{3M}{2}-1-n})-\sum _{n=0}^{\frac{M}{2}-1}{\epsilon }_n.(x_{\frac{3M}{2}-1-n}+x_{\frac{3M}{2}+n})$$则MDCT可以化为类似DCT变换的求和式:${\sum }_{n=0}^{M-1}\breve{x_n}{\epsilon }_n,\breve{x_n}$为:
   $$\breve{x_n}=\{\begin{array}{ll}-(x_{\frac{3M}{2}-1-n}+x_{\frac{3M}{2}+n}),& \text{ n = {0,..,2M-1 }}\\ (x_{n-\frac{M}{2}}-x_{\frac{3M}{2}-1-n}),,& \text{ n = {2M,..,M -1 }}\end{array}$$ 所以2N点的MDCT可以化为N/2的DCT进行计算，而DCT可以转化为FFT通过蝶形单元减少算法的计算时间度。