【频域变换】理论与Rust实现

DCT

二维 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);
    }
}
rust

矩阵化

发现事实：

用到的 cos 值是可以预先计算的，可以做一个 cosine_table
归一化系数 alpha 也是可以预先计算的
进而发现循环可以转化为矩阵积

看 dct 公式： $F(u) = \alpha(u) \sum\limits_{x=0}^{N-1} f(x) \cos \left( \dfrac{(2x + 1)u\pi}{2N} \right)$ ，后半部分的 $\sum$ 可以写成矩阵积，而前半部分的乘法可以写成矩阵元素积（Hadamard），如下：

[F(0), F(1),..., F(N)] \\ = ([f(0), f(1), ..., f(N)] \left[ \begin{array}{l} \cos \dfrac{(2\cdot 0 + 1)\cdot 0\pi}{2N} & \cos \dfrac{(2\cdot 0 + 1)\cdot 1\pi}{2N} & ... & \cos \dfrac{(2\cdot 0 + 1)\cdot u\pi}{2N} & ... & \cos \dfrac{(2\cdot 0 + 1)\cdot N\pi}{2N} \\ \cos \dfrac{(2\cdot 1 + 1)\cdot 0\pi}{2N} & \cos \dfrac{(2\cdot 1 + 1)\cdot 1\pi}{2N} & ... & \cos \dfrac{(2\cdot 1 + 1)\cdot u\pi}{2N} & ... & \cos \dfrac{(2\cdot 1 + 1)\cdot N\pi}{2N} \\ ...&...&...&...&...&...\\ \cos \dfrac{(2x + 1)\cdot 0\pi}{2N} & \cos \dfrac{(2x + 1)\cdot 1\pi}{2N} & ... & \cos \dfrac{(2x + 1)\cdot u\pi}{2N} & ... & \cos \dfrac{(2x + 1)\cdot N\pi}{2N} \\ ...&...&...&...&...&...\\ \cos \dfrac{(2N + 1)\cdot 0\pi}{2N} & \cos \dfrac{(2N + 1)\cdot 1\pi}{2N} & ... & \cos \dfrac{(2N + 1)\cdot u\pi}{2N} & ... & \cos \dfrac{(2N + 1)\cdot N\pi}{2N} \\ \end{array}\right] ) \\ \circ [\sqrt{\frac{1}{N}},\sqrt{\frac{2}{N}},...,\sqrt{\frac{2}{N}}]

说明

$\circ$ 表示对应元素积（Hadamard 积）
因为 $\alpha$ 的特殊性，这里面的 $\circ$ 运算实际上是初等列变换，是每列乘以常数。
- 这对于 2d-dct 也成立

上面的式子可以简写为：

dct: $F_{1\times n} = (f_{1\times n}\times \cos_{n\times n}) \circ \alpha_{1\times n}$ （矩阵运算）
idct $f= (F_{1\times n} \circ \alpha_{1\times n})\times \cos_{n\times n}^T$

而把 Hadamard 积写成矩阵积，如下：

dct： $F_{1\times n} = f_{1\times n}\times \cos_{n\times n} \times \Lambda_{n\times n}$
idct $f_{1\times n}= F_{1\times n} \times \Lambda_{n\times n} \times \cos_{n\times n}^T$

说明：

其中 $\Lambda = \mathrm{diag}(\alpha(0),\alpha(1),…,\alpha(N))$ ，是个对角矩阵
然后我们发现，可以证明 $\cos_{n\times n} \times \Lambda_{n\times n}$ 是一个正交矩阵。
若记为 $E = \cos_{n\times n} \times \Lambda_{n\times n}$ ，则有 $EE^T=I$
进而转化为
- dct: $F=f\times E$
- idct: $f=F\times E^T$

对于 2d-dct，是先对每一行做 dct，然后对每一列做 dct：

对每一行做dct： $T_1=f_{k\times n} D_{n\times n}$
对每一列做dct，我们借用对每一行做 dct 的方式，先转置，对每一行做dct，然后再转置回来： $T_2=(T_1^T D_{k\times k})^T$
$T_2 = D_{k\times k}^T f_{k\times n} D_{n\times n}$
因此：
- dct-2d： $D_{k\times k}^T f_{k\times n} D_{n\times n}$
- idct-2d: $D_{k\times k} f_{k\times n} D_{n\times n}^T$
说明
- 上面的 $D_{k\times k}$ 和 $D_{n\times n}$ 是不同的矩阵。
- 如果 f 是方阵（k=n），那么 $D_{k\times k} = D_{n\times n}$