修改型離散餘弦轉換

修改型離散餘弦轉換（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 : R2N -> RN（其中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個實數x0, ..., x2N-1轉換成N個實數X0, ..., XN-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個實數X0, ..., XN-1轉換成2N個實數y0, ..., y2N-1。（如同第四型離散餘弦轉換為一正交轉換，其逆轉換具有與正轉換一樣的形式）。具有一般窗式正規化（usual window normalization,如下所示）的窗式MDCT（windowed MDCT），其逆轉換的正規化係數必須乘2（即變成2/N）。

## 計算

雖然直接應用MDCT公式的計算複雜度為O（N2），但是透過遞迴的方式將計算分解化，可以將計算複雜度降到O（N log N），並且得出一樣的結果，就如同快速傅立葉轉換一樣。另外也可以透過其他轉換計算MDCT。通常一個離散傅立葉轉換（或快速傅立葉轉換）或是一個離散餘弦轉換包含的前置處理跟後置處理的複雜度為O(N)。此外，如下所述，任何第四型離散餘弦轉換的演演算法都能衍生出一種方法來計算偶數項的MDCT和IMDCT。

## 窗函數

在一般的信號壓縮應用中，轉換特性可透過利用將窗函數wn（n = 0, ..., 2N-1）乘以MDCT中的xn及IMDCT中的yn作進一步的改良，為了避免在n = 0 and 2N的邊界點產生不連續並且使函數於這些點可以變得平滑而趨近零（即window the data before the MDCT and after the IMDCT）。原則上，x和y可以有不同的窗函數，而窗函數在資料區塊轉換時可以改變（尤其是當兩個不同尺寸的資料區塊結合在一起時）。但為簡單起見，我們考慮給相同尺寸的資料區塊應用的窗函數。

轉換仍然是可逆的（即TDAC依然可以正常作用），考慮一個對稱的窗wn = w2N-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, –xR, –x, xR, ...），其中xR代表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點輸入（–cR–d, a–bR）的第四型離散餘弦轉換，其中R的意義與前述相同。（以這種方式，任何用來計算第四型離散餘弦轉換的演演算法都可以很直觀的套用在MDCT上）。同樣地，上述IMDCT公式恰好正是第四型離散餘弦轉換的1/2，其中輸出位移了N/2且擴增長度為2N（透過邊界條件）。而第四型離散餘弦逆轉換可以簡單的給定如上述之輸入（–cR–d, a–bR）。當透過邊界條件完成位移跟擴增的動作後，我們可以得到：

IMDCT(MDCT(a, b, c, d)) = (a–bR, b–aR, c+dR, cR+d) / 2

（因此有一半的IMDCT輸出是多餘的。）

我們現在明白了TDAC的運作原理。考慮計算一個MDCT的子帶，有50％的重疊，2N個資料區塊（c, d, e, f）。則得到的IMDCT會跟上述類似：(c–dR, d–cR, e+fR, eR+f) / 2。當加入這些到先前由重疊到一半所產生的IMDCT時，反向的項被消除了，並獲得一簡單的已回復的原始資料（c, d）。

## TDAC的起源

」時間域混疊消除法」的起源現在很清楚了。邏輯上第四型離散餘弦轉換中根據邊界條件所擴增的輸入資料造成資料發生混疊的現象，正如同頻譜中奈奎斯特頻率會與較低的頻率發生混疊一樣。因此，當c–dR及其他被加進來會有準確的符號（right sign）提供作消除之用。對於奇數點N（事實上很少使用），N/2不是整數，所以MDCT並不是由第四型離散餘弦轉換透過一個簡單的位移就可以互相對應。在這種情況下，額外的位移一半的樣本，代表著MDCT/IMDCT等同於第三/二型離散餘弦轉換，其分析類似上述。

## 窗式MDCT的TDAC

<image>

如上所述，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反向）。

因此，我們對$(w a,z b,z_R c,w_R d)$作MDCT取代對$(a,b,c,d)$作MDCT（所有的乘法跟加法均為點對點）。當完成IMDCT再乘以窗函數（點對點相乘），最後N點的一半，成為：

$$(z_R, w_R) \cdot (z_R c+w d_R, z c_R + w_R d) = (z_R^2 c+wz_R d_R, w_R z c_R + w_R^2 d)$$

（請注意，我們不再乘1/2，因為IMDCT的正規化與窗式的例子中的因數2有所不同）。同樣的，對於產生$(c,d,e,f)$的窗式MDCT和窗式IMDCT，是在其前N點的一半：

$$(w,z) \cdot (w c - z_R d_R, z d - w_R c_R) = (w^2 c - w z_R d_R, z^2 d - w_R z c_R)$$

當我們這兩個各半的結果加總起來，我們得到：

$$(z_R^2 c+wz_R d_R, w_R z c_R + w_R^2 d) + (w^2 c - w z_R d_R, z^2 d - w_R z c_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.