离散余弦变换

<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