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

背景

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

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

DFT到DCT的推导

DFT : X ( k ) = ∑ n = 0 N − 1 x [ n ] . e − j . 2 π . k n N   \large X(k) = \sum_{n=0}^{N-1} x[n].e^{\frac{-j.2\pi.kn}{N}} \text{ } X(k)=n=0N1x[n].eNj.2π.kn 
欧拉公式 : e − j θ = c o s θ + j . s i n θ \large e^{-j\theta} = cos\theta + j.sin\theta ejθ=cosθ+j.sinθ

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[N1]}。将该序列进行偶延拓,其数学定义更改为
    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]ˊ={x[m],x[m1],if n belong to { 0,..,N-1 }if n belong to { -N,..,-1 } 
    x [ m ] ˊ \acute{x[m]} x[m]ˊ信号如下图1所示:

    再将 x [ m ] ˊ \acute{x[m]} x[m]ˊ序列整体向右偏移 1 2 \Large\frac{1}{2} 21,令 x [ m ] ¨ \large\ddot{x[m]} x[m]¨ x [ m − 1 2 ] ˊ \large\acute{x[m-\frac{1}{2}]} x[m21]ˊ x [ m ] ¨ \large\ddot{x[m]} x[m]¨如下图2所示:
    在这里插入图片描述
  3. 重新推导实偶信号的DFT公式: X ( k ) = ∑ m = − N + 1 2 N − 1 2 x [ m − 1 2 ] ¨ . e − j . 2 π . k m 2 N   = 2 ∗ ∑ m = 1 2 N − 1 2 x [ m − 1 2 ] ¨ . e − j . 2 π . k m 2 N \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}} X(k)=m=N+21N21x[m21]¨.e2Nj.2π.km =2m=21N21x[m21]¨.e2Nj.2π.km令n = m + 1 2 \frac{1}{2} 21,则上式可化为 2 ∗ ∑ n = 0 N − 1 x [ n ] ˊ . c o s ( ( n + 1 2 ) . k π N ) \Large2*\sum_{n=0}^{N-1} \acute{x[n]}.cos(\frac{(n+\frac{1}{2}).k\pi}{N}) 2n=0N1x[n]ˊ.cos(N(n+21).kπ)
  4. 正交变换: 将DCT变换中与x[n]相乘的系数组织成矩阵C,如果该矩阵正交则有 C . C T = E C.C^{T} = E C.CT=E.故将变换核的系数2做变换可得下式: 2 N . g k ∗ ∑ n = 0 N − 1 x [ n ] ˊ . c o s ( ( n + 1 2 ) . 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} N2 .gkn=0N1x[n]ˊ.cos(N(n+21).kπ)(1)
    其中 g k g_k gk的数学定义为:
    g k = { 1 / 2 ,  k == 0 1 ,  k != 0  \large g_k = \begin{cases} 1/\sqrt{2}, & \text{ k == 0} \\ 1, & \text{ k != 0 } \end{cases} gk=1/2 ,1, k == 0 k != 0 

MDCT的编解码流程简述

  MDCT作为改进的离散余弦变换,所以编码由DCT过渡到MDCT是其本身的优势的。DCT在二维图片分量的变换中,其变换系数的高频分量集中在左上角(转换矩阵的左上角),而由于图片的编码是将整体图片切割成一个个小方块进行编码转换,更是造成了相邻方块间在转换之后容易引入噪声,这就是方块效应,在视觉上表示为图片编码后相邻小方块间的白条。
  而诸如LOTMDCT采用了TDAC实现的编码转换,转换后的单位抽样响应是由中间向其两边递减的,如下图3所示:
在这里插入图片描述
MDCT可以很好的消除方块效应。
  在MDCT变换中,输入的离散数字信号长度为2N,但是经过IMDCT[MDCT[x[n]]]的有效信号长度实则为N,下图4能很好的表示出来:
在这里插入图片描述

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

  1. MDCT变换公式:
    X ( k ) = 2 N ∗ ∑ n = 0 N − 1 x [ n ] . c o s [ 2 π N . ( n + 1 2 + N 4 ) . ( k + 1 2 ) ] 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\} $} X(k)=N2n=0N1x[n].cos[N2π.(n+21+4N).(k+21)]{0,..,N/21}
    在MDCT变换中,由于X(k) == X(N+k),所以X(k)只有N/2个独立分量,故k的范围为 { 0 , . . , N / 2 − 1 } \{ 0,..,N/2-1\} {0,..,N/21}
  2. IMDCT变换公式:
    x ( n ) = 2 ∗ ∑ k = 0 N 2 − 1 X [ k ] . c o s [ 2 π N . ( n + 1 2 + N 4 ) . ( k + 1 2 ) ] 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\} $} x(n)=2k=02N1X[k].cos[N2π.(n+21+4N).(k+21)]{0,..,N1}
  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]={x1x1ˊ,x2+x2ˊ}, x 1 ˊ \large\acute{x_1} x1ˊ x 1 \large x_1 x1的逆序序列,而 x 2 ˊ \large\acute{x_2} x2ˊ x 2 \large x_2 x2的逆序序列。
      令输入的离散信号长度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]={x0x1,x1x0,x2x3,x3x2}
      令长度N为4的MDCT变换矩阵为C,则C的数学定义如下:
    C k , n = [ C 0 , 0 C 1 , 0 C 0 , 1 C 1 , 1 C 0 , 2 C 1 , 2 C 0 , 3 C 1 , 3 ] , y [ n ] = x [ n ] . C . C T ⟹ y [ n ] = x [ n ] . ( C C T ) 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) Ck,n=C0,0C0,1C0,2C0,3C1,0C1,1C1,2C1,3,y[n]=x[n].C.CTy[n]=x[n].(CCT)
      再令 Q = C C T Q = CC^T Q=CCT,且Q为4*4矩阵,则上述证明转换为推导 Q 0 , 0 = 1 , Q 0 , 1 = − 1 Q_0,_0 = 1,Q_0,_1=-1 Q0,0=1,Q0,1=1。再N=4的情况下,C表示如下:
    C k , n = [ c o s 3 8 π c o s 9 8 π c o s 5 8 π c o s 15 8 π c o s 7 8 π c o s 21 8 π c o s 9 8 π c o s 27 8 π ] , c o s a . c o s b = c o s ( a + b ) + c o s ( 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} Ck,n=cos83πcos85πcos87πcos89πcos89πcos815πcos821πcos827π,cosa.cosb=2cos(a+b)+cos(ab)
       Q 0 , 0 = C 0 , 0 ∗ C 0 , 0 + C 1 , 0 ∗ C 1 , 0 ⟹ c o s 3 8 π . c o s 3 8 π + c o s 9 8 π . c o s 9 8 π ⟹ 1 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 Q0,0=C0,0C0,0+C1,0C1,0cos83π.cos83π+cos89π.cos89π1
       Q 0 , 1 = C 0 , 1 ∗ C 0 , 0 + C 1 , 1 ∗ C 1 , 0 ⟹ c o s 5 8 π . c o s 3 8 π + c o s 15 8 π . c o s 9 8 π ⟹ − 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 Q0,1=C0,1C0,0+C1,1C1,0cos85π.cos83π+cos815π.cos89π1
       Q 0 , 2 = C 0 , 2 ∗ C 0 , 0 + C 1 , 2 ∗ C 1 , 0 ⟹ c o s 7 8 π . c o s 3 8 π + c o s 21 8 π . c o s 9 8 π ⟹ 0 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 Q0,2=C0,2C0,0+C1,2C1,0cos87π.cos83π+cos821π.cos89π0
       Q 0 , 3 = C 0 , 3 ∗ C 0 , 0 + C 1 , 3 ∗ C 1 , 0 ⟹ c o s 9 8 π . c o s 3 8 π + c o s 27 8 π . c o s 9 8 π ⟹ 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 Q0,3=C0,3C0,0+C1,3C1,0cos89π.cos83π+cos827π.cos89π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\} [x0x1x2x3]1100C1,0C1,1C1,2C1,3C2,0C2,1C2,2C2,3C3,0C3,1C3,2C3,3={y0,y1,y2,y3}
   可得 y 0 = x 0 − x 1 \large y_0 = x_0 - x_1 y0=x0x1,后续 y 1 \large y_1 y1的推导读者可以自证。
   令 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˘))={ xixiˊ,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+1xi+1ˊ,xi+2+xi+2ˊ}
   再令输出序列的 y i , y i + 1 y_i,y_{i+1} yi,yi+1进行时域50%重叠,即可还原出2 x i + 1 \large x_{i+1} xi+1

MDCT的快速算法

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

  1. MDCT的公式如下:
    X ( k ) = 2 N ∗ ∑ n = 0 N − 1 x [ n ] . c o s [ 2 π N . ( n + 1 2 + N 4 ) . ( k + 1 2 ) ] 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\} $} X(k)=N2n=0N1x[n].cos[N2π.(n+21+4N).(k+21)]{0,..,N/21}
    令N=2m,可变换为:
    X ( k ) = 1 M ∗ ∑ n = 0 2 M − 1 x [ n ] . c o s [ π M . ( n + 1 2 + M 2 ) . ( k + 1 2 ) ] 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\} $} X(k)=M1n=02M1x[n].cos[Mπ.(n+21+2M).(k+21)]{0,..,N/21}
    同时令 ε n = c o s [ π M . ( n + 1 2 ) . ( k + 1 2 ) ] \varepsilon_n=cos[\frac{\pi}{M}.(n+\frac{1}{2}).(k+\frac{1}{2})] εn=cos[Mπ.(n+21).(k+21)],而 ε n \varepsilon_n εn即是DCT变换核的系数因子。现通过 ε n \varepsilon_n εn将MDCT化简:
    X ( k ) = 1 M ∗ ∑ n = 0 2 M − 1 x [ n ] . c o s [ 2 π N . ( n + 1 2 + M 2 ) . ( k + 1 2 ) ] = ∑ n = 0 2 M − 1 x n . ε M 2 + 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} X(k)=M1n=02M1x[n].cos[N2π.(n+21+2M).(k+21)]=n=02M1xn.ε2M+n
    = ∑ n = 0 M 2 − 1 [ x n . ε M 2 + n + x M − 1 − n . ε 3 M 2 − 1 − n + ε 2 M − 1 − n . x 3 M 2 − 1 − n + ε 2 M + n . x 3 M 2 + 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} =n=02M1[xn.ε2M+n+xM1n.ε23M1n+ε2M1n.x23M1n+ε2M+n.x23M+n](3)又因为 ε 2 M + n = ε 2 M − 1 − n = − ε n \large\varepsilon_{2M+n}=\varepsilon_{2M-1-n}=-\varepsilon_n ε2M+n=ε2M1n=εn,式3可化为:
    ∑ n = 0 M 2 − 1 [ x n . ε M 2 + n − x M − 1 − n . ε M 2 + n − ε n . x 3 M 2 − 1 − n − ε 2 M + n . x 3 M 2 + 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=02M1[xn.ε2M+nxM1n.ε2M+nεn.x23M1nε2M+n.x23M+n] = ∑ n = M / 2 M − 1 ε n . ( x n − M 2 − x 3 M 2 − 1 − n ) − ∑ n = 0 M 2 − 1 ε n . ( x 3 M 2 − 1 − n + x 3 M 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} ) =n=M/2M1εn.(xn2Mx23M1n)n=02M1εn.(x23M1n+x23M+n)则MDCT可以化为类似DCT变换的求和式: ∑ n = 0 M − 1 x n ˘ ε n , x n ˘ \large\sum_{n=0}^{M-1}\breve{x_n}\varepsilon_n,\breve{x_n} n=0M1xn˘εn,xn˘为:
    x n ˘ = { − ( x 3 M 2 − 1 − n + x 3 M 2 + n ) ,  n = {0,.., M 2 -1 } ( x n − M 2 − x 3 M 2 − 1 − n ) , ,  n = { M 2 ,.., 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} xn˘=(x23M1n+x23M+n),(xn2Mx23M1n),, n = {0,..,2M-1 } n = {2M,..,M -1 } 所以2N点的MDCT可以化为N/2的DCT进行计算,而DCT可以转化为FFT通过蝶形单元减少算法的计算时间度。
已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页