您的位置: 首页  >  文章  >  DCT和DFT的关系以及MDCT的推导证明

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

分类: 文章 • 2024-11-15 10:54:33

目录

背景

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

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

DFT到DCT的推导

DFT : X(k)=∑n=0N−1x[n].e−j.2π.knN \large X(k) = \sum_{n=0}^{N-1} x[n].e^{\frac{-j.2\pi.kn}{N}} \text{ }
欧拉公式 : e−jθ=cosθ+j.sinθ\large 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[m],if n belong to { 0,..,N-1 }x[−m−1],if n belong to { -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}
    x[m]ˊ\acute{x[m]}信号如下图1所示：
    DCT和DFT的关系以及MDCT的推导证明
    再将x[m]ˊ\acute{x[m]}序列整体向右偏移12\Large\frac{1}{2}，令x[m]¨\large\ddot{x[m]}为x[m−12]ˊ\large\acute{x[m-\frac{1}{2}]}，x[m]¨\large\ddot{x[m]}如下图2所示:
    DCT和DFT的关系以及MDCT的推导证明
  3. 重新推导实偶信号的DFT公式: X(k)=∑m=−N+12N−12x[m−12]¨.e−j.2π.km2N =2∗∑m=12N−12x[m−12]¨.e−j.2π.km2N\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 + 12\frac{1}{2},则上式可化为2∗∑n=0N−1x[n]ˊ.cos((n+12).kπN)\Large2*\sum_{n=0}^{N-1} \acute{x[n]}.cos(\frac{(n+\frac{1}{2}).k\pi}{N})
  4. 正交变换: 将DCT变换中与x[n]相乘的系数组织成矩阵C,如果该矩阵正交则有C.CT=EC.C^{T} = E.故将变换核的系数2做变换可得下式:2N.gk∗∑n=0N−1x[n]ˊ.cos((n+12).kπN)(1)\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}
    其中gkg_k的数学定义为:
    gk={1/2, k == 01, k != 0 \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所示:
DCT和DFT的关系以及MDCT的推导证明
故MDCT可以很好的消除方块效应。
  在MDCT变换中，输入的离散数字信号长度为2N，但是经过IMDCT[MDCT[x[n]]]的有效信号长度实则为N，下图4能很好的表示出来:
DCT和DFT的关系以及MDCT的推导证明

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

  1. MDCT变换公式：
    X(k)=2N∗∑n=0N−1x[n].cos[2πN.(n+12+N4).(k+12)]k ∈{0,..,N/2−1}\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}\{ 0,..,N/2-1\}
  2. IMDCT变换公式:
    x(n)=2∗∑k=0N2−1X[k].cos[2πN.(n+12+N4).(k+12)]n ∈{0,..,N−1}\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. 如何从解码端获取原始信号:
      假设输入信号的序列为x[n]={x1,x2}\large x[n]=\{ x_1,x_2 \},现证明经过MDCT变换以及IMDCT变换后的输出信号y[n]={x1−x1ˊ,x2+x2ˊ}\large y[n]=\{ x_1-\acute{x_1},x_2+\acute{x_2} \},x1ˊ\large\acute{x_1}为x1\large x_1的逆序序列,而x2ˊ\large\acute{x_2}为x2\large x_2的逆序序列。
      令输入的离散信号长度N为4，x[n]={x0,x1,x2,x3}\large x[n]=\{x_0,x_1,x_2,x_3\},则需证明y[n]={x0−x1,x1−x0,x2−x3,x3−x2}\large y[n] =\{ x_0 - x_1,x_1-x_0,x_2-x_3,x_3-x_2\}
      令长度N为4的MDCT变换矩阵为C，则C的数学定义如下:
    Ck,n=[C0,0C1,0C0,1C1,1C0,2C1,2C0,3C1,3],y[n]=x[n].C.CT⟹y[n]=x[n].(CCT) 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=CCTQ = CC^T，且Q为4*4矩阵，则上述证明转换为推导Q0,0=1,Q0,1=−1Q_0,_0 = 1,Q_0,_1=-1。再N=4的情况下，C表示如下:
    Ck,n=[cos38πcos98πcos58πcos158πcos78πcos218πcos98πcos278π],cosa.cosb=cos(a+b)+cos(a−b)2 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}
       Q0,0=C0,0∗C0,0+C1,0∗C1,0⟹cos38π.cos38π+cos98π.cos98π⟹1Q_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
       Q0,1=C0,1∗C0,0+C1,1∗C1,0⟹cos58π.cos38π+cos158π.cos98π⟹−1Q_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
       Q0,2=C0,2∗C0,0+C1,2∗C1,0⟹cos78π.cos38π+cos218π.cos98π⟹0Q_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
       Q0,3=C0,3∗C0,0+C1,3∗C1,0⟹cos98π.cos38π+cos278π.cos98π⟹0Q_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

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

MDCT的快速算法

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

  1. MDCT的公式如下:
    X(k)=2N∗∑n=0N−1x[n].cos[2πN.(n+12+N4).(k+12)]k ∈{0,..,N/2−1}\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，可变换为:
    X(k)=1M∗∑n=02M−1x[n].cos[piM.(n+12+M2).(k+12)]k ∈{0,..,N/2−1}\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\} $}
    同时令εn=cos[πM.(n+12).(k+12)]\varepsilon_n=cos[\frac{\pi}{M}.(n+\frac{1}{2}).(k+\frac{1}{2})],而εn\varepsilon_n即是DCT变换核的系数因子。现通过εn\varepsilon_n将MDCT化简:
    X(k)=1M∗∑n=02M−1x[n].cos[2πN.(n+12+M2).(k+12)]=∑n=02M−1xn.εM2+n\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}
    =∑n=0M2−1[xn.εM2+n+xM−1−n.ε3M2−1−n+ε2M−1−n.x3M2−1−n+ε2M+n.x3M2+n](3)\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} 又因为ε2M+n=ε2M−1−n=−εn\large\varepsilon_{2M+n}=\varepsilon_{2M-1-n}=-\varepsilon_n，式3可化为:
    ∑n=0M2−1[xn.εM2+n−xM−1−n.εM2+n−εn.x3M2−1−n−ε2M+n.x3M2+n]\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} ] =∑n=M/2M−1εn.(xn−M2−x3M2−1−n)−∑n=0M2−1εn.(x3M2−1−n+x3M2+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变换的求和式:∑n=0M−1xn˘εn,xn˘\large\sum_{n=0}^{M-1}\breve{x_n}\varepsilon_n,\breve{x_n}为:
    xn˘={−(x3M2−1−n+x3M2+n), n = {0,..,M2-1 }(xn−M2−x3M2−1−n),, n = {M2,..,M -1 }\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通过蝶形单元减少算法的计算时间度。

相关推荐