離散餘弦變換

<image>

離散餘弦變換（discrete cosine transform, DCT）是與傅立葉變換相關的一種變換，類似於離散傅立葉變換，但是只使用實數。離散餘弦變換相當於一個長度大概是它兩倍的離散傅立葉變換，這個離散傅立葉變換是對一個實偶函數進行的（因為一個實偶函數的傅立葉變換仍然是一個實偶函數），在有些變形裡面需要將輸入或者輸出的位置移動半個單位（DCT有8種標準類型，其中4種是常見的）。

最常用的一種離散餘弦變換的類型是下面給出的第二種類型，通常我們所說的離散餘弦變換指的就是這種。它的逆，也就是下面給出的第三種類型，通常相應的被稱為「反離散餘弦變換」，「逆離散餘弦變換」或者「IDCT」。

有兩個相關的變換，一個是離散正弦變換，它相當於一個長度大概是它兩倍的實奇函數的離散傅立葉變換；另一個是改進的離散餘弦變換，它相當於對交疊的數據進行離散餘弦變換。

## 應用

離散餘弦變換，尤其是它的第二種類型，經常被信號處理和圖像處理使用，用於對信號和圖像進行有損數據壓縮。這是由於離散餘弦變換具有很強的「能量集中」特性：大多數的信號資訊（包括聲音和圖像）往往集中在離散餘弦變換後的低頻部分，而且當信號具有接近馬爾可夫過程的統計特性時，離散餘弦變換的去相關性接近於K-L變換（Karhunen-Loève變換——它具有最優的去相關性）的性能。

例如，在圖像編碼標準JPEG與視訊編碼標準MJPEG和MPEG的各個標準中都使用了離散餘弦變換。在這些標準制中都使用了二維的第二種類型離散餘弦變換，並將結果進行量化之後進行熵編碼。這時對應第二種類型離散餘弦變換中的n通常是8，並用該公式對每個8x8塊的每行進行變換，然後每列進行變換。得到的是一個8x8的變換係數矩陣。其中（0,0）位置的元素就是直流分量，矩陣中的其他元素根據其位置表示不同頻率的交流分量。

一個類似的變換，改進的離散餘弦變換被用在高級音訊編碼、Vorbis和MP3音訊壓縮當中。

離散餘弦變換也經常被用來使用譜方法來解偏微分方程，這時候離散餘弦變換的不同的變量對應著陣列兩端不同的奇/偶邊界條件。

### 常見應用

- 音訊信號處理 — 音訊編碼、音訊資料壓縮（有損和無損）、環繞聲、回音消除、音位辨識、時域混疊消除法（TDAC）
- 數位音訊
- 數位廣播 — 數位聲音廣播（DAB+）
- 語音處理 — 語音編碼、語音辨識、語音活性檢測（VAD）
- 數位電話 — VoIP、行動電話、視訊通話
- 生物辨識技術 — 指紋定向、臉部辨識系統、生物辨識浮水印、掌紋辨識
- 人臉檢測

## 在科學計算與物理模型中的應用

近年來，離散餘弦轉換已被應用於處理高維度科學大數據，以克服數據壓縮與計算複雜度的限制。特別是在核融合科學領域，TokaMind 模型利用 3D 離散餘弦轉換 (3D-DCT-II) 技術，對托卡馬克（Tokamak）反應爐產生的四維時空電漿數據進行高效特徵提取。

該方法透過以下機制優化數據處理流程：

- 能量緊縮（Energy Compaction）： 3D-DCT-II 能在空間與時間維度上將電漿演化數據中的冗餘雜訊過濾，僅保留關鍵的物理特徵（如電漿不穩定性波形），大幅提升了多模態 Transformer (MMT) 基礎模型的訓練效率。
- 維度映射：透過 3D-DCT-II 將複雜的時空數據變換至頻譜空間，使得模型能夠在較低的運算開銷下，準確識別核融合實驗中極短暫的電漿擾動現象，為人工智慧輔助能源研究提供了關鍵的數據預處理標準。TokaMind 利用 3D-DCT-II 提取的時空特徵，已被驗證可遷移應用於智慧電網（Power Grid）的負載動態預測，以DCT-II可微分最佳化建立低頻能量與模態響應等指標之Fourier-based Low-Frequency Spectral Dynamics Network，實現科學計算與物理模型架構在不同物理系統之間的廣泛應用。

## 形式化定義

形式上來看，離散餘弦變換是一個線性的可逆函數$F:R^n\rightarrow R^n$其中R是實數集，或者等價的說一個$n \times n$的方陣。離散餘弦變換有幾種變形的形式，它們都是根據下面的某一個公式把$n$個實數$x_0,\ldots ,x_{n-1}$變換到另外$n$個實數$f_0,\ldots ,f_{n-1}$的操作。

### DCT-I

$$f_m = \frac{1}{2} (x_0 + (-1)^m x_{n-1}) <pre> + \sum_{k=1}^{n-2} x_k \cos \left[\frac{\pi}{n-1} m k \right]$$

有些人認為應該將$x_0$和$x_{n-1}$乘以$\sqrt{2}$，相應的將$f_0$和$f_{n-1}$乘以$\frac{1}{\sqrt{2}}$。這樣做的結果是這種DCT-I矩陣變為了正交矩陣（再乘一個係數的話），但是這樣就不能直接和一個實偶離散傅立葉變換對應了。

一個$n=5$的對實數abcde的DCT-I型變換等價於一個8點的對實數abcdedcb（偶對稱）的DFT變換，結果再除以2（對應的，DCT-II~DCT-IV相對等價的DFT有一個半個抽樣的位移）。需要指出的是，DCT-I不適用於$n<2$的情況（其它的DCT類型都適用於所有的整數n）。

所以，DCT-I暗示的邊界條件是：$x_k$相對於$k=0$點偶對稱，並且相對於$k=n-1$點偶對稱；對$f_m$的情況也類似。

### DCT-II

$$f_m = <pre> \sum_{k=0}^{n-1} x_k \cos \left[\frac{\pi}{n} m \left(k+\frac{1}{2}\right) \right]$$

DCT-II大概是最常用的一種形式，通常直接被稱為DCT。

有些人更進一步的將$f_0$再乘以$\frac{1}{\sqrt{2}}$（參見下面的DCT-III型的對應修改）。這將使得DCT-II成為正交矩陣（再乘一個係數的話），但是這樣就不能直接和一個有半個抽樣位移的實偶離散傅立葉變換對應了。

所以，DCT-II暗示的邊界條件是：$x_k$相對於$k=-\frac{1}{2}$點偶對稱，並且相對於$k=n-\frac{1}{2}$點奇對稱；對$f_m$相對於$m=0$點偶對稱，並且相對於$m=n$點奇對稱。

### DCT-III

$$f_m = \frac{1}{2} x_0 + <pre> \sum_{k=1}^{n-1} x_k \cos \left[\frac{\pi}{n} \left(m+\frac{1}{2}\right) k \right]$$

因為這是DCT-II的逆變換（再乘一個係數的話），這種變形通常被簡單的稱為逆離散餘弦變換。

有些人更進一步的將$x_0$再乘以$\sqrt{2}$（參見上面的DCT-II型的對應修改），這將使得DCT-III成為正交矩陣（再乘一個係數的話），但是這樣就不能直接和一個結果有半個抽樣位移的實偶離散傅立葉變換對應了。

所以，DCT-III暗示的邊界條件是：$x_k$相對於$k=0$點偶對稱，並且相對於$k=n$點奇對稱；對$f_m$相對於$m=-\frac{1}{2}$點偶對稱，並且相對於$m=n-\frac{1}{2}$點偶對稱。

### DCT-IV

$$f_m = <pre> \sum_{k=0}^{n-1} x_k \cos \left[\frac{\pi}{n} \left(m+\frac{1}{2}\right) \left(k+\frac{1}{2}\right) \right]$$

DCT-IV對應的矩陣是正交矩陣（再乘一個係數的話）。

一種DCT-IV的變形，將不同的變換的數據重疊起來，被稱為改進的離散餘弦變換。

DCT-IV暗示的邊界條件是：$x_k$相對於$k=-\frac{1}{2}$點偶對稱，並且相對於$k=n-\frac{1}{2}$點奇對稱；對$f_m$類似。

### DCT V~VIII

上面提到的DCT I~IV是和偶數階的實偶DFT對應的。原則上，還有四種DCT變換（Martucci, 1994）是和奇數階的實偶DFT對應的，它們在分母中都有一個$n+\frac{1}{2}$的係數。但是在實際應用中，這幾種變型很少被用到。

最平凡的和奇數階的實偶DFT對應的DCT是1階的DCT（1也是奇數），可以說變換隻是乘上一個係數$a$而已，對應於DCT-V的長度為1的狀況。

## 反變換

DCT-I的反變換是把DCT-I乘以係數$\frac{2}{n-1}$。

DCT-IV的反變換是把DCT-IV乘以係數$\frac{2}{n}$。

DCT-II的反變換是把DCT-III乘以係數$\frac{2}{n}$，反之亦然。

和離散傅立葉變換類似，變化前面的歸一化係數僅僅是常規而已，改變這個係數並不改變變換的性質。例如，有些人喜歡在DCT-II變換的前面乘以$\sqrt{\frac{2}{n}}$，這樣反變換從形式上就和變換更相似，而不需要另外的歸一化係數。

## 計算

儘管直接使用公式進行變換需要進行$O(n^2)$次操作，但是和快速傅立葉變換類似，我們有複雜度為$O(n \log(n))$的快速算法，這就是常常被稱做蝶形變換的一種分解演算法。另外一種方法是通過快速傅立葉變換來計算DCT，這時候需要$O(n)$的預操作和後操作。

以下簡單介紹兩種利用DFT來計算DCT-II的方法

### 方法一

令輸入信號為$x(n) \, , n=0,1,...,N-1$

並將$y(n)$以$x(n)$在$(2N-1)/2$處對稱表示

即$y(n)=\left\{\begin{matrix} x(n), & \mbox{if }n=0,1,...,N-1 \\ x(2N-n-1), & \mbox{if }n=N,N+1,...2N-1 \end{matrix}\right.$

此時令$W_{2N}$ 表示 $e^{\frac{-j2\pi}{2N}}$

則$y(n)$之DFT為

$$Y(m) = \Sigma_{n=0}^{2N-1}y(n)W_{2N}^{nm}$$

將$Y(m)$ 做以下化簡

$$\begin{align} Y(m) &= \sum_{n=0}^{N-1}y(n)W_{2N}^{nm} + \sum_{n=N}^{2N-1}y(n)W_{2N}^{nm} \\ & = \sum_{n=0}^{N-1}y(n)W_{2N}^{nm} + \sum_{n=N}^{2N-1}x(2N-n-1)W_{2N}^{nm} \\ & = \sum_{n=0}^{N-1}y(n)W_{2N}^{nm} + \sum_{n=0}^{N-1}x(n)W_{2N}^{(2N-n-1)m} \\ & = \sum_{n=0}^{N-1}x(n)[W_{2N}^{nm}+W_{2N}^{-(n+1)m}], \, \,\,\, m=0,1,...,2N-1 \end{align}$$

此時兩側同乘 $\frac{1}{2}W_{2N}^{m/2}$

可得$\frac{1}{2} W_{2N}^{m/2} Y(m) = \sum_{n=0}^{N-1} x(n)\cos{[(2n+1)\frac{m\pi}{2N}]}, \,\,\,\,\,\, m=0,1,...,N-1$

此時右式即為欲求之DCT轉換，而左式可藉由2N點數的DFT來計算，使用快速演演算法的情況下，運算之時間複雜度為$O(NlogN)$

### 方法二

第二個方法由Narasimha與Peterson在1978年提出，此方法係藉由巧妙的編排$y(n)$來達成，首先令

$y(n)=x(2n)$ 並且 $y(N-1-n)=x(2n+1), \,\, \,\,\,\, n=0,1,...,\frac{N}{2}-1$

此時X(m)可化簡為

$$X(m)= \sum_{n=0}^{N/2-1}y(n)\cos{[\frac{(4n+1)m\pi}{2N}]}+ \sum_{n=0}^{N/2-1}y(N-n-1)\cos{[\frac{(4n+3)m\pi}{2N}]}, \,\,\, \,\,\,\, m=0,1,...,N-1$$

令第二項之$n$改為 $n'=N-1-n$，則兩式可合併為

$$X(m)=\sum_{n=0}^{N-1}y(n)\cos{[\frac{(4n+1)m\pi}{2N}]}, \,\,\,\,\,\, m=0,1,...,N-1$$

右側為對$y(n)$之N點的scaled DFT

因此，$X(m)=Re[Z(m)]$，其中

$$Z(m)=W_{4N}^{m}Y(m)=W_{4N}^{m}\sum_{n=0}^{N-1}y(n)W_{N}^{nm},\,\,\, \,\,\,\, m=0,1,...,N-1$$

其中$Y(m)$ 是對$y(n)$之N點的DFT，並且可以簡單的驗證$Z(m)$具有如下性質

$$Z(N-m)=-jZ(m)^{*}$$

而因$y(n)$ 為實數輸入。

因此欲求之$X(m) = Re[Z(m)]$，$X(N-m) = -Im[Z(m)],\,\,\, \, \, \, \, m=0,1,...,\frac{N}{2}$

在使用FFT快速演演算法的情況下，運算之時間複雜度同樣為$O(NlogN)$

但此方法較方法一直接運算2N點數的DFT快上約2倍。

## 參考

- K. R. Rao and P. Yip, 離散餘弦變換：演算法、優點和應用（Discrete Cosine Transform: Algorithms, Advantages, Applications） (Academic Press, Boston, 1990).
- A. V. Oppenheim, R. W. Schafer, and J. R. Buck, 時間離散信號處理 (Discrete-Time Signal Processing), second edition (Prentice-Hall, New Jersey, 1999).
- S. A. Martucci, 對稱卷積和離散正弦餘弦變換 (Symmetric convolution and the discrete sine and cosine transforms), IEEE Trans. Sig. Processing SP-42, 1038-1051 (1994).
- Matteo Frigo and Steven G. Johnson: FFTW, http://www.fftw.org/ . 一個免費的C語言庫GPL，可以計算DCT-I~IV的1維到多維的任意大小的變換
- M. Frigo and S. G. Johnson, " FFTW3的設計和實現 ," Proceedings of the IEEE 93 (2), 216–231 (2005).
- On the Computation of the Discrete Cosine Transform. (1978, June 1). IEEE Journals & Magazine | IEEE Xplore. https://ieeexplore.ieee.org/document/1094144