用Python手写DCT变换：从傅里叶公式到图像压缩实战
用Python手写DCT变换：从傅里叶公式到图像压缩实战
如果你曾经好奇过，为什么一张几兆的JPEG图片在放大后会出现模糊的色块，而一张同样大小的PNG图片却能保持清晰，那么你其实已经触摸到了离散余弦变换（DCT）的门槛。在数字图像和视频压缩的世界里，DCT扮演着一位沉默的“能量搬运工”，它不动声色地将图像中肉眼不敏感的高频细节转化为可被大幅压缩的数据，从而让我们的手机能存储成千上万张照片，让在线视频流畅播放成为可能。
对于算法工程师和计算机视觉的初学者而言，仅仅调用cv2.dct()或scipy.fftpack.dct是远远不够的。真正理解DCT，意味着你能从最基础的傅里叶级数出发，亲手推导出变换矩阵，并用代码实现从空域到频域的完整映射。这个过程不仅能让你深刻理解JPEG、MPEG等主流压缩标准的底层逻辑，更能让你在遇到复杂的信号处理问题时，拥有从原理层面进行优化和创新的能力。今天，我们就抛开现成的库函数，在Jupyter Notebook中，一步步从离散傅里叶变换（DFT）的公式出发，手写实现DCT-II，并最终将其应用于一个8x8图像块的压缩实战，亲眼见证“能量集中”这一神奇特性。
1. 从离散傅里叶变换到离散余弦变换：原理推导
要理解离散余弦变换，我们必须先回到它的源头——离散傅里叶变换。DFT是信号从时域（或空域）转换到频域的核心工具，但其复数运算和周期性边界条件在处理实信号（如图像像素值）时，会引入不必要的计算复杂度和冗余信息。DCT正是为了解决这些问题而诞生的。
1.1 DFT的局限性与DCT的诞生
一个长度为N的实数序列x[n]，其N点DFT定义为：
import numpy as np
def naive_dft(x):
"""
手写实现DFT公式，用于理解原理。
注意：此实现计算复杂度为O(N^2)，仅用于教学。
"""
N = len(x)
X = np.zeros(N, dtype=complex)
for k in range(N):
for n in range(N):
X[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
return X
这个公式计算出的X[k]是一个复数序列。对于实信号x[n]，其DFT具有共轭对称性，即X[k] = conj(X[N-k])。这意味着近一半的频谱信息是冗余的。此外，DFT隐含的周期性假设（将有限长序列视为周期信号的一个周期）在处理非周期性图像块边界时，可能会在边界处引入不连续的高频分量，即所谓的“边界效应”。
提示：DFT的周期性意味着它默认信号的首尾是相连的。想象一下，如果一张图片的左边和右边像素值差异巨大，DFT会认为这是一个剧烈的跳变，从而产生大量不必要的高频能量，这不利于压缩。
DCT的核心思想是：将一个长度为N的实数序列，通过对称延拓，构造成一个长度为2N的偶对称序列，然后对这个新序列做2N点的DFT。由于构造后的序列是实偶对称的，其DFT结果将是一个实数序列，并且所有虚部为零，同时消除了由非周期性边界引起的虚假高频分量。最常用的DCT-II，其延拓方式可以直观理解为将原序列像镜子一样反射一遍。
1.2 DCT-II公式的矩阵化推导
DCT-II的定义式如下，对于长度为N的输入序列x[n]，输出系数X[k]为：
$$ X[k] = \sum_{n=0}^{N-1} x[n] \cdot \cos\left[\frac{\pi}{N} k (n + \frac{1}{2}) \right] \quad \text{for } k = 0, 1, ..., N-1 $$
注意，这里的k=0对应的是直流分量（DC coefficient），k>0对应的是交流分量（AC coefficients）。为了将其转化为矩阵运算，我们可以构造一个变换矩阵C，使得X = C * x。矩阵C的每个元素C[k, n]由下式给出：
$$ C[k, n] = \alpha(k) \cdot \cos\left[\frac{\pi}{N} k (n + \frac{1}{2}) \right] $$
其中，归一化因子α(k)通常有两种选择，对应正交变换和正交归一变换：
- 正交变换：
α(0) = sqrt(1/N),α(k) = sqrt(2/N)for k>0 - 正交归一变换（使变换矩阵是酉矩阵）：
α(k) = sqrt(2/N)for all k，但k=0时公式略有调整。
在JPEG标准中，通常使用正交变换形式。下面我们用Python来构建这个变换矩阵：
def dct_matrix(N, orthonormal=False):
"""
生成N点DCT-II的变换矩阵C。
Args:
N: 变换长度。
orthonormal: 如果为True，生成正交归一化矩阵（酉矩阵）。
Returns:
C: N x N的DCT变换矩阵。
"""
C = np.zeros((N, N))
for k in range(N):
for n in range(N):
# 计算余弦核
val = np.cos(np.pi * k * (2*n + 1) / (2 * N))
# 应用归一化因子
if orthonormal:
alpha = np.sqrt(2.0 / N)
if k == 0:
alpha = np.sqrt(1.0 / N)
C[k, n] = alpha * val
else:
# JPEG常用形式
if k == 0:
alpha = np.sqrt(1.0 / N)
else:
alpha = np.sqrt(2.0 / N)
C[k, n] = alpha * val
return C
# 生成一个8x8的DCT变换矩阵（正交归一化版本）
C8 = dct_matrix(8, orthonormal=True)
print("DCT变换矩阵C8的形状:", C8.shape)
print("验证正交性 (C8 * C8.T 应接近单位矩阵):")
print(np.round(C8 @ C8.T, 10)) # 忽略微小浮点误差
运行这段代码，你会得到一个8x8的矩阵。这个矩阵的每一行可以看作是一个基向量，对应一个特定频率的余弦波采样。DCT变换的本质，就是将原始信号x投影到这组完备的余弦基上，得到一组系数X，这组系数代表了原始信号在不同频率余弦分量上的“强度”。
2. 二维DCT与图像块的频域表示
图像是二维信号，因此我们需要将一维DCT扩展到二维。幸运的是，二维DCT是可分离的，这意味着我们可以先对图像的每一行做一维DCT，再对结果的每一列做一维DCT（或者先列后行，结果相同）。对于一个M x N的图像块f，其二维DCT变换F可以通过矩阵运算高效完成：
def dct2d_block(block):
"""
对单个图像块进行二维DCT变换。
使用可分离性质：先对行变换，再对列变换。
Args:
block: 一个M x N的numpy数组（图像块）。
Returns:
coeff: DCT系数矩阵。
"""
M, N = block.shape
# 生成变换矩阵
C_M = dct_matrix(M, orthonormal=True)
C_N = dct_matrix(N, orthonormal=True)
# 二维DCT: F = C_M * block * C_N^T
coeff = C_M @ block @ C_N.T
return coeff
def idct2d_block(coeff):
"""
对DCT系数矩阵进行二维逆变换，重建图像块。
Args:
coeff: DCT系数矩阵。
Returns:
block: 重建的图像块。
"""
M, N = coeff.shape
C_M = dct_matrix(M, orthonormal=True)
C_N = dct_matrix(N, orthonormal=True)
# 二维IDCT: block = C_M^T * coeff * C_N
block = C_M.T @ coeff @ C_N
return block
为了直观感受DCT的效果，我们找一个8x8的图像块来试试。这里我们不用真实图片，而是构造一个包含简单边缘和纹理的合成块：
import matplotlib.pyplot as plt
# 创建一个8x8的测试图像块，模拟一个从黑到白的渐变和一条边缘
test_block = np.zeros((8, 8))
for i in range(8):
test_block[i, :] = i * 32 # 垂直渐变
test_block[3:5, 2:6] = 200 # 中间一个亮块，模拟一个边缘
# 进行DCT变换
dct_coeff = dct2d_block(test_block)
# 可视化
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
im0 = axes[0].imshow(test_block, cmap='gray', vmin=0, vmax=255)
axes[0].set_title('原始8x8图像块')
plt.colorbar(im0, ax=axes[0])
# 显示DCT系数（取对数显示以看清细节）
dct_coeff_abs_log = np.log10(np.abs(dct_coeff) + 1e-10) # 避免log(0)
im1 = axes[1].imshow(dct_coeff_abs_log, cmap='hot')
axes[1].set_title('DCT系数幅值（对数刻度）')
plt.colorbar(im1, ax=axes[1])
# 重建图像
reconstructed_block = idct2d_block(dct_coeff)
im2 = axes[2].imshow(reconstructed_block, cmap='gray', vmin=0, vmax=255)
axes[2].set_title('重建图像块 (应无失真)')
plt.colorbar(im2, ax=axes[2])
# 检查重建误差
print(f"重建图像与原始图像的最大绝对误差: {np.max(np.abs(reconstructed_block - test_block)):.2e}")
print(f"均方误差 (MSE): {np.mean((reconstructed_block - test_block)**2):.2e}")
plt.tight_layout()
plt.show()
运行这段代码，你会看到三个图。第一个是原始块，第二个是其DCT系数。注意观察DCT系数矩阵的左上角，那里的值（尤其是F[0,0]，即DC系数）通常最大，代表了图像块的平均亮度。越往右下角，系数代表的频率越高，其幅值通常越小。这就是DCT的能量集中特性：对于典型的自然图像，其大部分视觉信息（能量）都集中在低频部分，即DCT系数矩阵的左上区域。
3. 能量集中与Zigzag扫描：JPEG压缩的核心
为什么能量集中特性对压缩如此重要？因为我们可以利用人类视觉系统对高频细节不敏感的特性，在几乎不影响主观质量的前提下，丢弃或粗量化这些高频系数，从而实现数据量的锐减。
3.1 量化：有损压缩的关键步骤
量化是JPEG压缩中唯一有损的步骤。其原理非常简单：将每个DCT系数除以一个对应的量化步长，然后取整。步长越大，量化越粗糙，压缩率越高，但失真也越大。JPEG标准为亮度（Y）分量和色度（Cb, Cr）分量分别定义了推荐的量化表。下面是一个标准的亮度量化表（8x8）：
# JPEG标准亮度量化表 (Quality = 50)
Q_luminance = np.array([
[16, 11, 10, 16, 24, 40, 51, 61],
[12, 12, 14, 19, 26, 58, 60, 55],
[14, 13, 16, 24, 40, 57, 69, 56],
[14, 17, 22, 29, 51, 87, 80, 62],
[18, 22, 37, 56, 68, 109, 103, 77],
[24, 35, 55, 64, 81, 104, 113, 92],
[49, 64, 78, 87, 103, 121, 120, 101],
[72, 92, 95, 98, 112, 100, 103, 99]
], dtype=np.float32)
def quantize(dct_coeff, quantization_table):
"""
对DCT系数进行量化。
Args:
dct_coeff: DCT系数矩阵。
quantization_table: 量化表。
Returns:
quantized_coeff: 量化后的系数（整数）。
"""
# 量化公式：Q = round(DCT / Q_table)
quantized = np.round(dct_coeff / quantization_table)
return quantized.astype(np.int32)
def dequantize(quantized_coeff, quantization_table):
"""
对量化后的系数进行反量化。
Args:
quantized_coeff: 量化后的系数（整数）。
quantization_table: 量化表。
Returns:
approx_dct_coeff: 反量化后的DCT系数（有损，是原始系数的近似）。
"""
return quantized_coeff * quantization_table
# 对我们的测试块进行量化和反量化
quantized = quantize(dct_coeff, Q_luminance)
restored_dct = dequantize(quantized, Q_luminance)
restored_block = idct2d_block(restored_dct)
# 计算量化带来的误差
quant_error = np.mean((restored_block - test_block) ** 2)
print(f"经过量化/反量化后的重建MSE: {quant_error:.2f}")
量化后，许多高频系数（量化表右下角数值大的区域）会变成0。这些连续的0为后续的熵编码（如霍夫曼编码）创造了极佳的压缩条件。
3.2 Zigzag扫描：将二维系数转换为一维序列
为了更高效地对量化后系数进行编码，JPEG采用了一种称为“Zigzag扫描”的技术，将8x8的二维系数矩阵按照“之”字形顺序重新排列成一维序列。这样做的目的是将非零系数（主要集中左上角）聚集在序列的前部，而大量的零系数则集中在序列尾部，便于行程编码（Run-Length Encoding, RLE）。
| 扫描顺序索引 | 对应矩阵位置 (行, 列) | 频率特点 |
|---|---|---|
| 0 | (0,0) | DC系数，最低频 |
| 1 | (0,1) | 低频 |
| 2 | (1,0) | 低频 |
| 3 | (2,0) | 中低频 |
| ... | ... | ... |
| 63 | (7,7) | 最高频 |
下面我们实现这个扫描过程：
def zigzag_scan(matrix):
"""
对8x8矩阵进行Zigzag扫描，返回一维序列。
"""
rows, cols = matrix.shape
assert rows == cols == 8, "目前只支持8x8矩阵"
order = []
for d in range(rows + cols - 1):
if d % 2 == 0: # 偶数对角线，方向向上
row = min(d, rows - 1)
col = d - row
while row >= 0 and col < cols:
order.append(matrix[row, col])
row -= 1
col += 1
else: # 奇数对角线，方向向下
col = min(d, cols - 1)
row = d - col
while col >= 0 and row < rows:
order.append(matrix[row, col])
row += 1
col -= 1
return np.array(order)
def inverse_zigzag_scan(vector):
"""
将Zigzag扫描后的一维序列还原为8x8矩阵。
"""
matrix = np.zeros((8, 8), dtype=vector.dtype)
rows, cols = 8, 8
idx = 0
for d in range(rows + cols - 1):
if d % 2 == 0:
row = min(d, rows - 1)
col = d - row
while row >= 0 and col < cols and idx < len(vector):
matrix[row, col] = vector[idx]
row -= 1
col += 1
idx += 1
else:
col = min(d, cols - 1)
row = d - col
while col >= 0 and row < rows and idx < len(vector):
matrix[row, col] = vector[idx]
row += 1
col -= 1
idx += 1
return matrix
# 对量化后的系数进行Zigzag扫描
zigzag_seq = zigzag_scan(quantized)
print("量化后系数的Zigzag序列 (前20个):")
print(zigzag_seq[:20])
print(f"序列中零值的数量: {np.sum(zigzag_seq == 0)} / {len(zigzag_seq)}")
运行后你会发现，序列末尾出现了大量的0。在实际的JPEG编码中，这个一维序列会经过差分脉冲编码调制（DPCM）处理DC系数（因为相邻块的DC系数通常很接近），并对AC系数进行行程编码（RLE），最后再使用霍夫曼编码或算术编码进行熵编码，从而获得极高的压缩比。
4. 实战：构建一个完整的图像块压缩与可视化工具
现在，让我们把所有环节串联起来，构建一个可以交互式探索DCT压缩效果的小工具。我们将使用OpenCV读取一张灰度图像，将其分割成8x8的块，对每个块进行DCT、量化、反量化、IDCT，并对比原始图像与压缩后图像的差异。
import cv2
def jpeg_compress_block(block, q_table, quality_factor=50):
"""
模拟JPEG对一个8x8块的处理流程。
Args:
block: 8x8图像块，像素值范围建议为[0, 255]。
q_table: 基础量化表。
quality_factor: 质量因子 (1-100)，用于调整量化强度。
Returns:
compressed_block: 压缩后重建的块。
quantized_coeff: 量化后的系数（可用于计算压缩率）。
"""
# 1. 电平偏移：将像素值从[0,255]平移到[-128, 127]，使数据围绕0对称，有利于DCT。
block_shifted = block.astype(np.float32) - 128.0
# 2. 计算二维DCT
dct_coeff = dct2d_block(block_shifted)
# 3. 根据质量因子调整量化表
# 质量因子越高，量化步长越小，质量越好
if quality_factor < 50:
scale = 5000.0 / quality_factor if quality_factor > 0 else 5000.0
else:
scale = 200.0 - 2.0 * quality_factor
scaled_q_table = np.floor((q_table * scale + 50) / 100.0)
scaled_q_table = np.clip(scaled_q_table, 1, 255) # 步长至少为1
# 4. 量化
quantized = quantize(dct_coeff, scaled_q_table)
# 5. 反量化
restored_dct = dequantize(quantized, scaled_q_table)
# 6. 逆DCT
restored_block_shifted = idct2d_block(restored_dct)
# 7. 反向电平偏移
compressed_block = np.clip(restored_block_shifted + 128.0, 0, 255).astype(np.uint8)
return compressed_block, quantized
def process_full_image(image_path, quality=75):
"""
处理整张图像，展示压缩效果。
"""
# 读取图像并转为灰度图
img = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
if img is None:
print(f"无法读取图像: {image_path}")
return
h, w = img.shape
# 调整图像尺寸为8的倍数
h_new = (h // 8) * 8
w_new = (w // 8) * 8
img = img[:h_new, :w_new]
# 创建输出图像
compressed_img = np.zeros_like(img, dtype=np.uint8)
total_nonzero = 0
total_coeffs = 0
# 分块处理
for i in range(0, h_new, 8):
for j in range(0, w_new, 8):
block = img[i:i+8, j:j+8]
compressed_block, quantized_coeff = jpeg_compress_block(block, Q_luminance, quality)
compressed_img[i:i+8, j:j+8] = compressed_block
# 统计非零系数，粗略估计压缩率
total_nonzero += np.count_nonzero(quantized_coeff)
total_coeffs += 64
# 计算指标
mse = np.mean((img.astype(float) - compressed_img.astype(float)) ** 2)
psnr = 10 * np.log10(255.0**2 / mse) if mse > 0 else float('inf')
compression_ratio_est = total_coeffs / max(total_nonzero, 1) # 非常粗略的估计
# 可视化
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
axes[0].imshow(img, cmap='gray')
axes[0].set_title(f'原始图像 ({w_new}x{h_new})')
axes[0].axis('off')
axes[1].imshow(compressed_img, cmap='gray')
axes[1].set_title(f'压缩后图像 (质量因子={quality})')
axes[1].axis('off')
diff = np.abs(img.astype(int) - compressed_img.astype(int))
im_diff = axes[2].imshow(diff, cmap='hot', vmin=0, vmax=50)
axes[2].set_title('差异图 (放大显示)')
axes[2].axis('off')
plt.colorbar(im_diff, ax=axes[2], fraction=0.046, pad=0.04)
plt.suptitle(f'PSNR: {psnr:.2f} dB | 非零系数占比: {total_nonzero/total_coeffs*100:.1f}% | 粗略压缩比: {compression_ratio_est:.2f}:1')
plt.tight_layout()
plt.show()
# 使用示例：你需要准备一张名为'test_image.jpg'的图片，或者使用下面的代码生成一个测试图
# 生成一个简单的测试图
test_img = np.zeros((256, 256), dtype=np.uint8)
cv2.putText(test_img, 'DCT Test', (50, 120), cv2.FONT_HERSHEY_SIMPLEX, 2, 255, 3)
cv2.imwrite('dct_test_image.jpg', test_img)
# 处理图像，尝试不同的质量因子
process_full_image('dct_test_image.jpg', quality=90)
process_full_image('dct_test_image.jpg', quality=50)
process_full_image('dct_test_image.jpg', quality=10)
这个工具可以让你清晰地看到，随着质量因子的降低（量化更粗糙），图像的细节逐渐丢失，出现典型的“块效应”和模糊，但文件大小会显著减小。PSNR（峰值信噪比）指标客观地反映了失真程度，而非零系数的占比则直观地展示了DCT的能量集中特性——即使在高质量设置下，大部分系数经过量化后也变成了0。
通过这个从公式推导到代码实现，再到完整应用案例的旅程，你应该已经对DCT在图像压缩中的核心作用有了深刻的理解。手写实现这些算法，虽然效率远不及高度优化的FFT库，但其教育意义是无可替代的。它让你穿透了“黑箱”，看到了从连续的余弦波到离散的像素矩阵，再到高效的二进制编码之间那条清晰而优美的路径。下次当你保存一张JPEG图片时，或许会会心一笑，因为你已经知道，是成千上万个这样的8x8小块，经过DCT的“提纯”和量化的“取舍”，最终安静地躺在你的硬盘里。
更多推荐
12
0- 0
已为社区贡献3条内容
所有评论(0)