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. 步骤1:对每行应用一维DCT
2. 步骤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);
    }
}

矩阵化

发现事实：

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

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

说明

• $\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}$

DWT（离散小波变换）

DWT-1D

1. 基本概念：1D DWT 将信号（或数据向量）分解为近似系数（低频部分, ca）和细节系数（高频部分, cd）。
2. 数学原理：通过一对低通和高通滤波器对信号进行卷积，然后下采样（去除部分数据）来实现信号的多分辨率分析。具体步骤如下：
   • 输入信号 $x(t)$ 与小波函数进行卷积，得到一组系数。
   • 将信号分解为两部分：低频部分（近似部分）和高频部分（细节部分）。
3. 结果：最终的结果是一组小波系数，表示信号在不同尺度或分辨率下的特征。

公式：

其中，$A_k$ 为低频部分，$D_k$ 为高频部分，
$\phi_k[n]$ 和 $\psi_k[n]$ 分别为低通和高通滤波器函数：

对于 Haar:

其中， \(h = \left[\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}\right] \\ g = \left[-\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}}\right]\)

DWT2（2D 离散小波变换）：

1. 基本概念：2D DWT 是在二维数据（如图像）上应用 1D DWT，
2. 算法
   1. 首先对图像的每一行进行 1D DWT，得到两个矩阵：ca_row 和 cd_row
   2. 再对每个矩阵的每一列进行 1D DWT，得到4个矩阵 ca, ch, cv, cd
3. 这样将信号分解为 4 个子带（各为原图 1/4 像素数）
   • LL（cA, 低频近似）：表示图像在水平和垂直方向上的低频部分。
   • LH（cH, 水平细节部分）：表示图像在水平方向上的高频部分，也就是垂直边缘
   • HL（cV, 垂直细节部分）：表示图像在垂直方向上的高频部分，也就是水平边缘。
   • HH（cD, 对角细节部分）：表示图像在水平和垂直方向上的高频部分，也就是纹理/噪声。

import numpy as np


# 保证偶数个元素。如果不是，用最后一个填充
def pad_1d(signal):
    if signal.shape[0] % 2 == 1:
        signal = np.append(signal, signal[-1])
    return signal


def dwt1d(signal):
    signal = pad_1d(signal)

    signal = signal / np.sqrt(2)
    signal1 = signal[::2]
    signal2 = signal[1::2]

    ca = signal1 + signal2
    cd = signal1 - signal2

    return ca, cd


def idwt1d(ca, cd):
    signal = np.zeros(ca.shape[0] * 2)
    sig1 = ca + cd
    sig2 = ca - cd
    signal[::2] = sig1
    signal[1::2] = sig2
    signal = signal / np.sqrt(2)
    return signal


def pad_2d(arr: np.ndarray) -> np.ndarray:
    # 保证偶数行、偶数列
    res = arr
    # 补行
    if res.shape[0] % 2 != 0:
        last_row = res[-1:, :]  # shape (1, n)
        res = np.vstack([res, last_row])
    # 补列
    if res.shape[1] % 2 != 0:
        last_col = res[:, -1:]  # shape (n, 1)
        res = np.hstack([res, last_col])
    return res


def dwt2d(img: np.ndarray):
    img = pad_2d(img)
    # 对行进行 dwt
    img = img / 2  # 本应是 sqrt(2)，但下面还有一个 sqrt(2)
    ca_row = img[:, ::2] + img[:, 1::2]
    cd_row = img[:, ::2] - img[:, 1::2]

    # 对ca_row、cd_row 的列进行 dwt
    ca = ca_row[::2, :] + ca_row[1::2, :]
    ch = ca_row[::2, :] - ca_row[1::2, :]
    cv = cd_row[::2, :] + cd_row[1::2, :]
    cd = cd_row[::2, :] - cd_row[1::2, :]

    return ca, ch, cv, cd


def idwt2d(ca, ch, cv, cd):
    img = np.zeros((ca.shape[0] * 2, ca.shape[1] * 2))
    ca_row = np.zeros((ca.shape[0] * 2, ca.shape[1]))
    cd_row = np.zeros((ca.shape[0] * 2, ca.shape[1]))

    ca_row[::2, :] = (ca + ch) / 2
    ca_row[1::2, :] = (ca - ch) / 2
    cd_row[::2, :] = (cv + cd) / 2
    cd_row[1::2, :] = (cv - cd) / 2

    img[:, ::2] = (ca_row + cd_row) / 2
    img[:, 1::2] = (ca_row - cd_row) / 2
    return img * 2


# %% 验证与 PyWavelets 的 dwt 一致
import pywt

for sig_len in range(4, 20):
    for _ in range(10):
        signal = np.random.rand(sig_len) * 10
        signal1 = signal.copy()

        ca, cd = dwt1d(signal)
        pywt_ca, pywt_cd = pywt.dwt(signal, 'haar')

        signal_raw = idwt1d(ca, cd)

        assert np.abs(pywt_ca - ca).max() < 1e-7
        assert np.abs(pywt_cd - cd).max() < 1e-7
        assert np.abs(signal_raw[:sig_len] - signal).max() < 1e-7

# %% 验证与 PyWavelets 的 dwt2 一致
for sig_row in range(3, 10):
    for sig_line in range(3, 10):
        img = np.random.rand(sig_row, sig_line)

        ca, ch, cv, cd = dwt2d(img)

        p_cA, (p_cH, p_cV, p_cD) = pywt.dwt2(img, 'haar')

        assert np.abs(ca - p_cA).max() < 1e-7
        assert np.abs(ch - p_cH).max() < 1e-7
        assert np.abs(cv - p_cV).max() < 1e-7
        assert np.abs(cd - p_cD).max() < 1e-7

        img_recover = idwt2d(ca, ch, cv, cd)

        assert np.abs(img - img_recover[:sig_row, :sig_line]).max() < 1e-7