修改型离散余弦变换

修改型离散余弦变换（英语：modified discrete cosine transform，下称缩写：MDCT）是一种派生自傅里叶变换的变换。MDCT以第四型离散余弦变换（DCT-IV）为基础，但具有重叠变换的性质：它在较大资料集合中，用于处理连续资料区块，其中当前资料区块的后半段与下一个资料区块的前半段相重叠。这种重叠情形有助于避免资料区块边界产生多余资料，除了具有离散余弦变换的能量压缩特性外，也使MDCT尤其有利于数据压缩方面应用。借由以上优势，MDCT在音频压缩领域中成为应用最为广泛的有损数据压缩手段，应用于当前大多数音频编码格式，如MP3、杜比数字、Ogg Vorbis、AAC等等。

MDCT是由Princen、Johnson和Bradley于1987年提出，其基本原理即时域混叠消除法（英语：time-domain aliasing cancellation，缩写：TDAC）是由Princen和Bradley于1986年提出。其他类似的变换还有如以离散正弦变换为基础的修改型离散正弦变换（英语：modified discrete sine transform，缩写：MDST）以及其他较少使用的变换，例如以其他不同类型的DCT或DCT/DST的组合为基础的MDCT。

定义

作为一种重叠变换，MDCT与其他傅里叶相关的变换比起有些不一样，因为它的输出数为输入数的一半（而非相同）。特别的是，它是一个线性函数F : R^2N -> R^N（其中R代表实数集合）。根据公式：

X_{k}=\sum _{n=0}^{2N-1}x_{n}\cos \left[{\frac {\pi }{N}}\left(n+{\frac {1}{2}}+{\frac {N}{2}}\right)\left(k+{\frac {1}{2}}\right)\right]

将2N个实数x₀, ..., x_2N-1变换成N个实数X₀, ..., X_N-1。（在这个变换前面的正规化系数，单位不为一定且根据不同的解释方法而有所改变。下面只考虑MDCT跟IMDCT的正规化乘积。）

逆变换

MDCT之逆变换为离散余弦逆变换（英语：inverse modified discrete cosine transform，缩写：IMDCT）。

MDCT由于输入数与输出数并不相同，乍看之下应不可逆。然而，透过加入造成错误并且必须要回复原始资料的子带重叠区块中重叠的IMDCT部分，MDCT可具备完全可逆性。这种技术称为时域混叠消除法。

IMDCT根据公式：

y_{n}={\frac {1}{N}}\sum _{k=0}^{N-1}X_{k}\cos \left[{\frac {\pi }{N}}\left(n+{\frac {1}{2}}+{\frac {N}{2}}\right)\left(k+{\frac {1}{2}}\right)\right]

将N个实数X₀, ..., X_N-1变换成2N个实数y₀, ..., y_2N-1。（如同第四型离散余弦变换为一正交变换，其逆变换具有与正变换一样的形式）。具有一般窗式正规化（usual window normalization,如下所示）的窗式MDCT（windowed MDCT），其逆变换的正规化系数必须乘2（即变成2/N）。

计算

虽然直接应用MDCT公式的计算复杂度为O（N²），但是透过递归的方式将计算分解化，可以将计算复杂度降到O（N log N），并且得出一样的结果，就如同快速傅里叶变换一样。另外也可以透过其他变换计算MDCT。通常一个离散傅里叶变换（或快速傅里叶变换）或是一个离散余弦变换包含的前置处理跟后置处理的复杂度为O(N)。此外，如下所述，任何第四型离散余弦变换的算法都能派生出一种方法来计算偶数项的MDCT和IMDCT。

窗函数

在一般的信号压缩应用中，变换特性可透过利用将窗函数w_n（n = 0, ..., 2N-1）乘以MDCT中的x_n及IMDCT中的y_n作进一步的改良，为了避免在n = 0 and 2N的边界点产生不连续并且使函数于这些点可以变得平滑而趋近零（即window the data before the MDCT and after the IMDCT）。原则上，x和y可以有不同的窗函数，而窗函数再资料区块变换时可以改变（尤其是当两个不同尺寸的资料区块结合在一起时）。但为简单起见，我们考虑给相同尺寸的资料区块应用的窗函数。变换仍然是可逆的（即TDAC依然可以正常作用），考虑一个对称的窗w_n = w_2N-1-n，只要w满足Princen-Bradley condition：

w_{n}^{2}+w_{n+N}^{2}=1

各种不同应用的窗函数都是一样的，例如：

w_{n}=\sin \left[{\frac {\pi }{2N}}\left(n+{\frac {1}{2}}\right)\right]

为应用于MP3和MPEG - 2的AAC的窗函数，而

w_{n}=\sin \left({\frac {\pi }{2}}\sin ^{2}\left[{\frac {\pi }{2N}}\left(n+{\frac {1}{2}}\right)\right]\right)

为Vorbis所应用。此外AC-3与MPEG - 4 AAC均采用了Kaiser-Bessel derived (KBD) window。值得注意的是，应用在MDCT的窗函数与应用于其他类型的信号分析的窗函数是不同的。因为它们必须满足Princen – Bradley condition。原因之一，MDCT应用了两次窗函数，包括MDCT（分析）和IMDCT（合成）。

第四型离散余弦变换与原始TDAC的关系

由定义可以看出，一个偶数点N的MDCT基本上就相当于一个输入位移了N/2且两个N点资料区块同时变换的第四型离散余弦变换。而再更进一步的查看这个等效关系，则重要的特性例如TDAC可以很容易地推导出来。为了更明确的定义与第四型离散余弦变换的关系，我们必须了解第四型离散余弦变换分别对应到偶数点跟奇数点的边界条件：偶数点在左边边界，（约在n=–1/2），奇数点在右边边界（约在n=N–1/2）诸如此类（而不是DFT那样的周期性边界）。这些遵从以下所述：

\cos \left[{\frac {\pi }{N}}\left(-n-1+{\frac {1}{2}}\right)\left(k+{\frac {1}{2}}\right)\right]=\cos \left[{\frac {\pi }{N}}\left(n+{\frac {1}{2}}\right)\left(k+{\frac {1}{2}}\right)\right]

和

\cos \left[{\frac {\pi }{N}}\left(2N-n-1+{\frac {1}{2}}\right)\left(k+{\frac {1}{2}}\right)\right]=-\cos \left[{\frac {\pi }{N}}\left(n+{\frac {1}{2}}\right)\left(k+{\frac {1}{2}}\right)\right]

因此，如果输入为一组长度为N的x数组，我们可以想见这个数组会扩增至（x, –x_R, –x, x_R, ...），其中x_R代表x的反向序列。考虑一个有2N个输入与N个输出的MDCT。我们将输入分成4组区块（a, b, c, d），每组大小为N/2。如果我们把这些资料位移N/2（即MDCT的定义中的+N/2项），然后（b, c, d）扩增到N点第四型离散余弦变换的输入结尾。因此，我们必须将这些资料根据上述边界条件“叠”（“fold”）回来。是故，一个具有2N点输入（a, b, c, d）的MDCT，就确实的相当于一个具有N点输入（–c_R–d, a–b_R）的第四型离散余弦变换，其中R的意义与前述相同。（以这种方式，任何用来计算第四型离散余弦变换的算法都可以很直观的套用在MDCT上）。同样地，上述IMDCT公式恰好正是第四型离散余弦变换的1/2，其中输出位移了N/2且扩增长度为2N（透过边界条件）。而第四型离散余弦逆变换可以简单的给定如上述之输入（–c_R–d, a–b_R）。当透过边界条件完成位移跟扩增的动作后，我们可以得到：

IMDCT(MDCT(a, b, c, d)) = (a–b_R, b–a_R, c+d_R, c_R+d) / 2

（因此有一半的IMDCT输出是多余的。）

我们现在明白了TDAC的运作原理。考虑计算一个MDCT的子带，有50％的重叠，2N个资料区块（c, d, e, f）。则得到的IMDCT会跟上述类似：(c–d_R, d–c_R, e+f_R, e_R+f) / 2。当加入这些到先前由重叠到一半所产生的IMDCT时，反向的项被消除了，并获得一简单的已恢复的原始资料（c, d）。

TDAC的起源

”时间域混叠消除法”的起源现在很清楚了。逻辑上第四型离散余弦变换中根据边界条件所扩增的输入资料造成资料发生混叠的现象，正如同频谱中奈奎斯特频率会与较低的频率发生混叠一样。因此，当c–d_R及其他被加进来会有准确的符号（right sign）提供作消除之用。对于奇数点N（事实上很少使用），N/2不是整数，所以MDCT并不是由第四型离散余弦变换透过一个简单的位移就可以互相对应。在这种情况下，额外的位移一半的样本，代表着MDCT/IMDCT等同于第三/二型离散余弦变换，其分析类似上述。

窗式MDCT的TDAC

MDCT window functions:
blue: Cosine, red: Sine-Cosine, green: modified Kaiser-Bessel

如上所述，TDAC性质已被证明适用于普通MDCT，显示加入子带资料区块的IMDCT于重叠的一半部分可回复原始的资料。推导这项窗式MDCT逆性质较上述稍微复杂一些，叙述如下。根据前述当 $(a,b,c,d)$ 和 $(c,d,e,f)$ 经过MDCT与IMDCT，并且加上其重叠一半时，我们得到 $(c+d_{R},c_{R}+d)/2+(c-d_{R},d-c_{R})/2=(c,d)$ ，即原始的资料。现在我们假设我们将MDCT的输入和IMDCT的输出两者均乘上一个长度为2N的窗函数。如前所述，我们假设一个对称的窗函数其形式为 $(w,z,z_{R},w_{R})$ 其中w跟z是长度N/2的向量，而R代表反向。如此则Princen – Bradley condition可以写成： $w^{2}+z_{R}^{2}=(1,1,\ldots )$ ，其中乘法和加法为点对点的运算，或等效成 $w_{R}^{2}+z^{2}=(1,1,\ldots )$ （将w跟z反向）。因此，我们对 $(wa,zb,z_{R}c,w_{R}d)$ 作MDCT取代对 $(a,b,c,d)$ 作MDCT（所有的乘法跟加法均为点对点）。当完成IMDCT再乘以窗函数（点对点相乘），最后N点的一半，成为：

(z_{R},w_{R})\cdot (z_{R}c+wd_{R},zc_{R}+w_{R}d)=(z_{R}^{2}c+wz_{R}d_{R},w_{R}zc_{R}+w_{R}^{2}d)

（请注意，我们不再乘1/2，因为IMDCT的正规化与窗式的例子中的因数2有所不同）。同样的，对于产生 $(c,d,e,f)$ 的窗式MDCT和窗式IMDCT，是在其前N点的一半：

(w,z)\cdot (wc-z_{R}d_{R},zd-w_{R}c_{R})=(w^{2}c-wz_{R}d_{R},z^{2}d-w_{R}zc_{R})

当我们这两个各半的结果加总起来，我们得到：

(z_{R}^{2}c+wz_{R}d_{R},w_{R}zc_{R}+w_{R}^{2}d)+(w^{2}c-wz_{R}d_{R},z^{2}d-w_{R}zc_{R})

=\left([z_{R}^{2}+w^{2}]c+[wz_{R}-wz_{R}]d_{R},[w_{R}^{2}+z^{2}]d+[w_{R}z-w_{R}z]c_{R}\right)=(c,d)

即恢复的原始资料。

音频编码应用

就MP3格式而言，MDCT并不直接用于处理音频信号，而是用于处理32波段多相正交滤波器（英语：polyphase quadrature filter，缩写：PQF）数组输出端信号。MDCT输出结果会再经由一条混叠削减公式作后处理，以减弱PQF当中典型的混叠。这种变换与滤波器数组组合，被称作混合滤波器数组或子带MDCT。

AAC则相反，通常纯粹采用MDCT；其中唯一例外是索尼的罕见格式MPEG – 4 AAC - SSR——既采用MDCT，又在其后用到四波段PQF数组。自适应听觉变换编码（英语：adaptive transform acoustic coding，缩写：ATRAC）利用运用MDCT的堆栈型正交镜像滤波器（英语：quadrature mirror filter，缩写：QMF）。

参考

Henrique S. Malvar, Signal Processing with Lapped Transforms (Artech House: Norwood MA, 1992).

John P. Princen and Alan B. Bradley, "Analysis/synthesis filter bank design based on time domain aliasing cancellation," IEEE Trans. Acoust. Speech Sig. Proc. ASSP-34 (5), 1153-1161 (1986).（Described a precursor to the MDCT using a combination of discrete cosine and sine transforms.）

J. P. Princen and A. W. Johnson and A. B. Bradley, "Subband/transform coding using filter bank designs based on time domain aliasing cancellation," IEEE Proc. Intl. Conf. on Acoustics, Speech, and Signal Processing（ICASSP）12, 2161-2164 (1987).（Initial description of what is now called the MDCT.）

A. W. Johnson and A. B. Bradley, "Adaptive transform coding incorporating time domain aliasing cancellation," Speech Comm. 6, 299-308 (1987).

For algorithms, see e.g.:
- Chi-Min Liu and Wen-Chieh Lee, "A unified fast algorithm for cosine modulated filterbanks in current audio standards^{[永久失效链接]}", J. Audio Engineering 47 (12), 1061-1075 (1999).
- V. Britanak and K. R. Rao, "A new fast algorithm for the unified forward and inverse MDCT/MDST computation," Signal Processing 82, 433-459（2002）
- Vladimir Nikolajevic and Gerhard Fettweis, "Computation of forward and inverse MDCT using Clenshaw's recurrence formula," IEEE Trans. Sig. Proc. 51 (5), 1439-1444（2003）
- Che-Hong Chen, Bin-Da Liu, and Jar-Ferr Yang, "Recursive architectures for realizing modified discrete cosine transform and its inverse," IEEE Trans. Circuits Syst. II: Analog Dig. Sig. Proc. 50 (1), 38-45（2003）
- J.S. Wu, H.Z. Shu, L. Senhadji, and L.M. Luo, "Mixed-radix algorithm for the computation of forward and inverse MDCTs," IEEE Trans. Circuits Syst. I: Reg. Papers. 56 (4), 784-794（2009）
- V. Britanak, "A survey of efficient MDCT implementations in MP3 audio coding standard: retrospective and state-of-the-art," Signal. Process. 91 (4), 624-672(2011)
- ...and references thereof.