离散余弦变换(DCT)详解
目录
- 引言
- 背景知识:从傅里叶变换到余弦变换
- 为什么要用余弦?
- DCT的多种类型
- DCT-II公式详解
- 逆变换(DCT-III)
- 二维DCT
- DCT的快速算法简述
- DCT的性质与应用
- DCT计算示例(详细版)
- 实现代码(Python示例)
- 代码简要解读
引言
离散余弦变换(Discrete Cosine Transform,DCT)是一种在数字信号处理和图像处理领域具有重要地位的变换方式。它能够将信号从时域(或空间域)转化到以余弦基函数为核心的频域表达形式。DCT最大的特点之一是它只包含实数部分,没有虚部,计算结果往往比傅里叶变换更紧凑,在图像和音频压缩应用中极其常见。
背景知识:从傅里叶变换到余弦变换
- 连续傅里叶变换(CFT):对一个连续的函数 f ( t ) f(t) f(t),可以用无穷多的正弦波和余弦波的线性组合来表示,这就产生了傅里叶级数或傅里叶变换。
- 离散傅里叶变换(DFT):当信号是离散的(离散采样),我们可以使用离散傅里叶变换来将其表示为有限个复数频率分量。它常被快速算法(FFT)实现。
- 离散余弦变换(DCT):DFT得到的结果是复数,而在很多实际应用(如图像处理)里,信号本身通常是实值并且我们需要的数据也常常集中在余弦分量上。DCT正是专门用“余弦基函数”来对信号做正交分解的一种特殊形式,且结果均为实数。
从“傅里叶变换”到“离散傅里叶变换(DFT)”再到“离散余弦变换(DCT)”,我们可以简单理解为:
DCT是对某些特定边界条件下的DFT进行简化得出的,且只取余弦(实数)部分的一种变换。
为什么要用余弦?
- 平滑性:许多自然信号(例如图像中的连续区域、音频中的平滑片段)在用余弦函数来进行逼近时会更容易展现出“能量集中在低频”这一特征。
- 对称边界条件:DCT可以等效地看作原信号在边界处做了对称延拓,因此在边界过渡方面相对更“平滑”,不会像离散傅里叶变换那样在信号区间首尾可能出现明显的不连续跳变。
DCT的多种类型
根据变换公式中对信号边界条件、下标(索引)定义的不同,DCT有8种变体:DCT-I、DCT-II、DCT-III、DCT-IV、DCT-V、DCT-VI、DCT-VII、DCT-VIII。
- DCT-II:应用最广泛,也被很多文献直接称为“DCT”。
- DCT-III:DCT-II的逆变换。
- DCT-IV:在某些滤波或特殊音频编码技术中会使用。
这里主要讨论最常用的DCT-II(与它的逆变换DCT-III),在JPEG、MP3等主流标准中都是对DCT-II做各种改进或变形(如MDCT)来实现高效压缩。
DCT-II公式详解
一维DCT-II
对于长度为 N N N 的离散序列 { x [ n ] } \{ x[n] \} {x[n]}, n = 0 , 1 , … , N − 1 n=0,1,\ldots,N-1 n=0,1,…,N−1,其DCT-II系数 { X [ k ] } \{ X[k] \} {X[k]} 的一般定义为:
X [ k ] = α ( k ) ∑ n = 0 N − 1 x [ n ] ⋅ cos [ π N ( n + 1 2 ) k ] , k = 0 , 1 , … , N − 1 X[k] = \alpha(k) \sum_{n=0}^{N-1} x[n] \cdot \cos \left[ \frac{\pi}{N} \left(n + \frac{1}{2}\right) k \right], \quad k=0,1,\ldots,N-1 X[k]=α(k)n=0∑N−1x[n]⋅cos[Nπ(n+21)k],k=0,1,…,N−1
其中,
α
(
k
)
=
{
1
N
,
k
=
0
,
2
N
,
k
≠
0.
\alpha(k) = \begin{cases} \sqrt{\frac{1}{N}}, & \quad k = 0, \\ \sqrt{\frac{2}{N}}, & \quad k \neq 0. \end{cases}
α(k)=⎩
⎨
⎧N1,N2,k=0,k=0.
- α ( k ) \alpha(k) α(k) 是归一化因子(或缩放因子),它能够使DCT和逆DCT保持正交性和单位能量。
- cos [ π N ( n + 1 2 ) k ] \cos \left[ \frac{\pi}{N}\left(n+\frac{1}{2}\right)k \right] cos[Nπ(n+21)k] 表示第 k k k 个余弦基函数在第 n n n 个点的值。
名词解释:
- 正交:不同频率的基函数彼此之间相互独立(在积分或求和意义上为0)。
- 正交归一化:除了正交外,还要求不同基函数在自身内积后为1,而不是一个其他常数。
逆变换(DCT-III)
一维DCT-II的逆变换被称为DCT-III,其公式如下:
x [ n ] = ∑ k = 0 N − 1 β ( k ) X [ k ] ⋅ cos [ π N ( n + 1 2 ) k ] , n = 0 , 1 , … , N − 1 x[n] = \sum_{k=0}^{N-1} \beta(k) \, X[k] \cdot \cos\left[\frac{\pi}{N}\left(n+\frac{1}{2}\right)k\right], \quad n=0,1,\ldots,N-1 x[n]=k=0∑N−1β(k)X[k]⋅cos[Nπ(n+21)k],n=0,1,…,N−1
其中
β ( k ) = { 1 N , k = 0 , 2 N , k ≠ 0. \beta(k) = \begin{cases} \sqrt{\frac{1}{N}}, & k=0, \\ \sqrt{\frac{2}{N}}, & k\neq 0. \end{cases} β(k)=⎩ ⎨ ⎧N1,N2,k=0,k=0.
也就是说,只要我们用对应的逆变换系数 β ( k ) \beta(k) β(k) 和同样的余弦基函数,就能从频域系数 { X [ k ] } \{X[k]\} {X[k]} 复原出原始信号 { x [ n ] } \{x[n]\} {x[n]}。
二维DCT
在图像处理中,更常用的是“二维DCT”,其思路与一维的完全一致,只不过我们需要分别在行和列两个维度上做变换。
对于大小为 M × N M \times N M×N 的图像块(或矩阵) { f ( m , n ) } \{f(m,n)\} {f(m,n)},其二维DCT-II可表示为:
F
(
u
,
v
)
=
α
(
u
)
α
(
v
)
∑
m
=
0
M
−
1
∑
n
=
0
N
−
1
f
(
m
,
n
)
cos
[
π
M
(
m
+
1
2
)
u
]
cos
[
π
N
(
n
+
1
2
)
v
]
F(u,v) = \alpha(u)\alpha(v) \sum_{m=0}^{M-1}\sum_{n=0}^{N-1} f(m,n)\cos\left[\frac{\pi}{M}\left(m+\frac{1}{2}\right)u\right] \cos\left[\frac{\pi}{N}\left(n+\frac{1}{2}\right)v\right]
F(u,v)=α(u)α(v)m=0∑M−1n=0∑N−1f(m,n)cos[Mπ(m+21)u]cos[Nπ(n+21)v]
其中
α
(
u
)
=
{
1
M
,
u
=
0
2
M
,
u
≠
0
α
(
v
)
=
{
1
N
,
v
=
0
2
N
,
v
≠
0
\alpha(u) = \begin{cases} \sqrt{\frac{1}{M}}, & u=0 \\ \sqrt{\frac{2}{M}}, & u\neq 0 \end{cases} \quad \alpha(v) = \begin{cases} \sqrt{\frac{1}{N}}, & v=0 \\ \sqrt{\frac{2}{N}}, & v\neq 0 \end{cases}
α(u)=⎩
⎨
⎧M1,M2,u=0u=0α(v)=⎩
⎨
⎧N1,N2,v=0v=0
- 这里的 α ( u ) \alpha(u) α(u) 和 α ( v ) \alpha(v) α(v) 与一维情况下的缩放因子相同,只是针对行、列分别做归一化处理。
- 计算二维DCT最常见的做法是先对每行做一维DCT,再对变换后的结果按列再做一次一维DCT;反之亦然也行。
DCT的快速算法简述
就像DFT可以通过FFT得到快速实现一样,DCT也有相应的快速算法(Fast DCT Algorithms),其复杂度常可以降到 O ( N log N ) O(N \log N) O(NlogN) 级别。
- 思路:将余弦基的矩阵分解为若干子块,可以使用类似“蝴蝶操作”来节省计算量。
- 实际实现:在很多编程库(如FFTW、Intel MKL、OpenCV、scipy)中,都已经提供了高效的DCT实现,直接调用即可。
DCT的性质与应用
能量集中性
对大多数“平滑”信号,低频分量往往数值较大,而高频分量往往比较小。DCT后,大部分信号的能量会集中在前几个低频系数上,这为后续的“压缩”或“滤波”提供了良好的基础——只要把高频分量进行“粗量化”或干脆“截断”就能大幅减小数据量。
与DFT的对比
- 相似性:都能将时域/空间域信号转到频域;都依赖于正交基函数。
- 差异性:
- DFT得到的是复数系数(同时携带正弦和余弦信息),而DCT仅仅使用余弦函数,因此变换结果是实数;
- DCT的边界条件是对信号做镜像扩展,而DFT默认是周期性扩展。
在图像压缩中的作用(以JPEG为例)
- 图像分块:将图像分割成 8 × 8 8\times 8 8×8 的小块。
- 对每个小块做二维DCT。
- 量化:对DCT系数按不同频率分量重要性进行量化;通常低频分量更重要,保留更多精度,而高频分量可以采用大步长量化甚至舍弃。
- 熵编码:对量化后的系数再做哈夫曼编码或其他可逆压缩算法。
在音频压缩中的作用(以MP3为例)
- MP3中使用MDCT(Modified Discrete Cosine Transform),它是一种基于DCT的改进型变换,把信号做分帧后对每帧进行变换压缩。
- 余弦变换本身能很好地配合人耳听觉模型,将不敏感的频段做更强的量化。
DCT计算示例(详细版)
下面举一个更详细的一维DCT-II计算的小例子。
假设我们有一个长度为
N
=
8
N=8
N=8 的信号序列:
x
=
[
12
,
10
,
8
,
8
,
8
,
8
,
10
,
12
]
.
x = [12, 10, 8, 8, 8, 8, 10, 12].
x=[12,10,8,8,8,8,10,12].
-
手动计算零频分量 X [ 0 ] X[0] X[0]
X [ 0 ] = 1 8 ∑ n = 0 7 x [ n ] cos [ π 8 ( n + 1 2 ) ⋅ 0 ] = 1 8 ∑ n = 0 7 x [ n ] ⋅ 1. X[0] = \sqrt{\frac{1}{8}} \sum_{n=0}^{7} x[n] \cos\left[\frac{\pi}{8}(n+\frac{1}{2}) \cdot 0\right] = \sqrt{\frac{1}{8}} \sum_{n=0}^{7} x[n] \cdot 1. X[0]=81n=0∑7x[n]cos[8π(n+21)⋅0]=81n=0∑7x[n]⋅1.
因为 cos ( 0 ) = 1 \cos(0) = 1 cos(0)=1,这里就相当于取均值再乘 1 8 \sqrt{\frac{1}{8}} 81。 -
计算一般频率分量 X [ k ] X[k] X[k], k>0
X [ k ] = 2 8 ∑ n = 0 7 x [ n ] cos [ π 8 ( n + 1 2 ) k ] . X[k] = \sqrt{\frac{2}{8}} \sum_{n=0}^{7} x[n] \cos\left[\frac{\pi}{8}\left(n+\frac{1}{2}\right)k\right]. X[k]=82n=0∑7x[n]cos[8π(n+21)k].
可以把常数 2 8 = 1 4 = 1 2 \sqrt{\frac{2}{8}} = \sqrt{\frac{1}{4}} = \frac{1}{2} 82=41=21 拿到外面,后面就是一个纯粹的余弦乘法求和。 -
在实际中,我们通常不会像上面这样一个个算,而是直接用程序库或手写的快速DCT算法完成变换。
实现代码(Python示例)
import numpy as np
from scipy.fftpack import dct, idct
# 1. 准备测试数据(1D信号)
x = np.array([12, 10, 8, 8, 8, 8, 10, 12], dtype=float)
print("原始信号 x:", x)
# 2. 一维DCT计算(使用DCT-II)
X = dct(x, type=2, norm='ortho')
print("DCT系数 X:", X)
# 3. 逆变换(DCT-III)
x_reconstructed = idct(X, type=2, norm='ortho')
print("重构后的信号 x_reconstructed:", x_reconstructed)
# 4. 二维DCT示例:对一个简单8x8矩阵做DCT
matrix_2d = np.array([
[12, 10, 8, 8, 8, 8, 10, 12],
[12, 10, 8, 8, 8, 8, 10, 12],
[12, 10, 8, 8, 8, 8, 10, 12],
[12, 10, 8, 8, 8, 8, 10, 12],
[12, 10, 8, 8, 8, 8, 10, 12],
[12, 10, 8, 8, 8, 8, 10, 12],
[12, 10, 8, 8, 8, 8, 10, 12],
[12, 10, 8, 8, 8, 8, 10, 12],
], dtype=float)
# 对行进行DCT
dct_rows = dct(matrix_2d, axis=1, type=2, norm='ortho')
# 再对列进行DCT
dct_2d = dct(dct_rows, axis=0, type=2, norm='ortho')
print("二维DCT结果:\n", dct_2d)
# 逆变换:先对列做IDCT,再对行做IDCT
idct_cols = idct(dct_2d, axis=0, type=2, norm='ortho')
idct_2d = idct(idct_cols, axis=1, type=2, norm='ortho')
print("二维DCT重构后的矩阵:\n", idct_2d)