二维 dct 广泛用于图片频域相关算法，例如 JPEG 压缩
一维离散余弦变换（DCT-II）公式 $F(u) = \alpha(u) \sum\limits_{x=0}^{N-1} f(x) \cos \left( \dfrac{(2x + 1)u\pi}{2N} \right)$
其中，归一化系数 \(\alpha(u) = \begin{cases} \sqrt{\frac{1}{N}} & \text{if } u = 0 \\ \sqrt{\frac{2}{N}} & \text{otherwise} \end{cases}\)
一维离散余弦逆变换（DCT-III）公式 $f(x) = \sum\limits_{u=0}^{N-1} \alpha(u) F(u) \cos \left( \dfrac{(2x + 1)u\pi}{2N} \right)$
二维 DCT 可以由一维 DCT 得出，步骤如下：
- 步骤1:对每行应用一维DCT
- 步骤2:对每列应用一维DCT
代码
use std::f64::consts::PI;
/// 计算 DCT 的归一化系数 C(u) 或 C(v)
fn c(u: usize, n: usize) -> f64 {
if u == 0 {
(1.0 / n as f64).sqrt()
} else {
(2.0 / n as f64).sqrt()
}
}
/// 一维 DCT-II
fn dct_1d(vector: &[f64]) -> Vec<f64> {
let n = vector.len();
let mut result = vec![0.0; n];
for u in 0..n {
let mut sum = 0.0;
for x in 0..n {
sum += vector[x] * ((2 * x + 1) as f64 * u as f64 * PI / (2.0 * n as f64)).cos();
}
result[u] = c(u, n) * sum;
}
result
}
/// 一维 IDCT-III
fn idct_1d(vector: &[f64]) -> Vec<f64> {
let n = vector.len();
let mut result = vec![0.0; n];
for x in 0..n {
let mut sum = 0.0;
for u in 0..n {
sum += c(u, n) * vector[u] * ((2 * x + 1) as f64 * u as f64 * PI / (2.0 * n as f64)).cos();
}
result[x] = sum;
}
result
}
/// 二维 DCT
fn dct_2d(matrix: &[Vec<f64>]) -> Vec<Vec<f64>> {
let n = matrix.len();
let mut intermediate = vec![vec![0.0; n]; n];
let mut dct_matrix = vec![vec![0.0; n]; n];
// 对每一行应用一维 DCT
for i in 0..n {
intermediate[i] = dct_1d(&matrix[i]);
}
// 对每一列应用一维 DCT
for j in 0..n {
let column: Vec<f64> = intermediate.iter().map(|row| row[j]).collect();
let dct_col = dct_1d(&column);
for i in 0..n {
dct_matrix[i][j] = dct_col[i];
}
}
dct_matrix
}
/// 二维 IDCT
fn idct_2d(matrix: &[Vec<f64>]) -> Vec<Vec<f64>> {
let n = matrix.len();
let mut intermediate = vec![vec![0.0; n]; n];
let mut idct_matrix = vec![vec![0.0; n]; n];
// 对每一列应用一维 IDCT
for j in 0..n {
let column: Vec<f64> = matrix.iter().map(|row| row[j]).collect();
let idct_col = idct_1d(&column);
for i in 0..n {
intermediate[i][j] = idct_col[i];
}
}
// 对每一行应用一维 IDCT
for i in 0..n {
idct_matrix[i] = idct_1d(&intermediate[i]);
}
idct_matrix
}
#[test]
fn tst3() {
// 示例 4x4 矩阵（例如灰度图像块）
let matrix = vec![
vec![52.0, 55.0, 61.0, 66.0],
vec![70.0, 61.0, 64.0, 73.0],
vec![63.0, 59.0, 55.0, 90.0],
vec![67.0, 61.0, 68.0, 104.0],
];
println!("原始矩阵:");
for row in &matrix {
println!("{:?}", row);
}
// 计算二维 DCT
let dct = dct_2d(&matrix);
println!("\nDCT 矩阵:");
for row in &dct {
println!("{:?}", row);
}
// 计算二维 IDCT
let idct = idct_2d(&dct);
println!("\nIDCT 矩阵:");
for row in &idct {
println!("{:?}", row);
}
}
性能优化
- 用到的 cos 值是可以预先计算的，可以做一个
cosine_table - 归一化系数
alpha也是可以预先计算的 - 使用固定大小数组（
[f64; 4]; 4]），而非动态数组Vec<Vec<f64>>，可以提高缓存命中率 - 使用
packed_simd之类的库，可以加速向量化计算 - 使用并行计算
- 做成 table 后，发现迭代可以转化为矩阵积
迭代转化为矩阵积
- dct: $F_{1\times n} = (f_{1\times n}\times \cos_{n\times n}) .\times \alpha_{1\times n}$ （矩阵运算）
- idct $f= (F_{1\times n}.\times \alpha)\times \cos_{n\times n}^T$
关于 2d dct 的矩阵化：
- 对每一行做 dct，对应的矩阵积是 $ F_{k\times n} = (f_{k\times n}\times \cos_{n\times n}) .\times \alpha_{k\times n}$，这里的 $\alpha_{k\times n}$ 是 $\alpha_{1\times n}$ 的行复制k遍。可以用 repeat 生成，也可以左诚一个矩阵 $[1,1,…,1]^T$ 得到
- 列也一样
