3.2. 离散余弦变换类字典¶

3.2.1. 什么是离散余弦变换¶

离散余弦变换 (Discrete Cosine Transform, DCT) 将有限长序列信号表示成不同频率的余弦函数之和, 被广泛用于音频与图像的压缩. 目前共有 DCT-I, DCT-II, DCT-III, DCT-IV, DCT-V-VIII 五种变种 2. DCT在图像压缩上的应用内容参见 1. 最长见的是 DCT-II 型, 下面进行介绍.

3.2.2. 一维离散余弦变换¶

link text

离散余弦变换 (Discrete Cosine Transform, DCT)最早由 Ahmed等人 [8] 于1974年提出, 他们发现DCT的基提供了对Toeplitz矩阵特征向量的良好近似. 序列信号 \(x[k], k=0,\cdots,N-1\) 的 DCT-II 的离散余弦变换定义如下

(3.18)¶\[y[k] = {\sum}_{n=0}^{N-1}x[n]{\rm cos}\frac{(2n+1)k\pi}{2N}, \]

其中, \(y[k]\) 为变换后的系数序列. 对 equ-DCTII-1D 中的结果分别乘以 \(1/\sqrt{2}\) (\(k=0\)) 和 \(\sqrt{2/N}\) (\(k=0,1,\cdots,N-1\)) 可以得到正交版本的DCT变换.

Definition 3.2 (离散余弦变换 (Discreate Cosine Transformation))

离散序列 \(x[n], n=0, 1,\cdots, N-1\) 的DCT变换定义为

(3.19)¶\[y[k] = {\rm DCT}(x[n]) = \left\{ {\begin{array}{lll} {\sqrt{\frac{2}{N}}\sum_{n=0}^{N-1}x[n]\frac{1}{\sqrt 2}, \quad k=0}\\ {\sqrt{\frac{2}{N}}\sum_{n=0}^{N-1}x[n]{\rm cos}\frac{(2n + 1)k\pi}{2N}, \quad k=1, \cdots, N-1} \end{array}} \right. \]

用矩阵表示为

(3.20)¶\[{\bm y} = {\bm D}{\bm x} \]

其中, \({\bm x} = (x_n)_{N\times 1}, x_n = x[n]\), \({\bm D} = (d_{ij})_{N\times N}\) 表示为

(3.21)¶\[{\bm D} = \sqrt{\frac{2}{N}}\left[\begin{array}{ccccc}{1/\sqrt{2}} & {1/\sqrt{2}} & {1/\sqrt{2}} & {\cdots} & {1/\sqrt{2}} \\ {\cos \frac{\pi}{2 N}} & {\cos \frac{3 \pi}{2 N}} & {\cos \frac{5 \pi}{2 N}} & {\cdots} & {\cos \frac{(2 N-1) \pi}{2 N}} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ {\cos \frac{(N-1) \pi}{2 N}} & {\cos \frac{3(N-1) \pi}{2 N}} & {\cos \frac{5(N-1) \pi}{2 N}} & {\cdots} & {\cos \frac{(2 N-1)(N-1) \pi}{2 N}}\end{array}\right] \]

向量 \({\bm x}\) 为信号 \({\bm y} = (y_k)_{N\times 1}, y_k = y[k]\) 在基 \(\{1/{\sqrt 2}, {\rm cos}\frac{(2n + 1)k\pi}{2N}\}\) 下的坐标, 亦即离散余弦变换系数, 代表频率成分(与傅立叶变换一致). 变换前后的数据维度不变.

易知, 一维离散余弦变换矩阵 \(\bm D\) 为正交矩阵 (\({\bm D}{\bm D}^T = {\bm D}^T {\bm D} = {\bm I}\)), 故逆一维离散余弦变换可以通过下式计算

(3.22)¶\[{\bm x} = {\bm D}^{-1}{\bm y} = {\bm D}^T{\bm y} \]

3.2.3. 多维离散余弦变换¶

根据一维离散余弦变换的定义, 很容易将其扩展成多维形式, 记一个 \(L\) 维序列信号为 \(x_{k_1,k_2,\cdots,k_L}\), 其中 \(k_l = 0,1,\cdots,N_l-1\), \(l=1,2,\cdots,L\), \(N_l\) 为信号在第 \(l\) 为的样点数. 则扩展后的多维DCT-II变换可表示为

(3.23)¶\[\begin{aligned} y_{k_1,k_2,\cdots,k_L} = {\sum}_{n_1=0}^{N_1-1}\cdots \left({\sum}_{n_L=0}^{N_L-1} x_{n_1,\cdots,n_L}{\rm cos}\frac{(2n_L+1)k_L\pi}{2N_L}\right)\cdots{\rm cos}\frac{(2n_1+1)k_1\pi}{2N_1}\\ = {\sum}_{n_1=0}^{N_1-1}\cdots {\sum}_{n_L=0}^{N_L-1} x_{n_1,\cdots,n_L} {\rm cos}\frac{(2n_1+1)k_1\pi}{2N_1}\cdots {\rm cos}\frac{(2n_L+1)k_L\pi}{2N_L}. \end{aligned} \]

如二维的DCT变换可表示为

(3.24)¶\[\begin{aligned} y_{k_1,k_2} &= {\sum}_{n_1=0}^{N_1-1}\left({\sum}_{n_2=0}^{N_2-1}x_{n_1,n_2}{\rm cos}\frac{(2n_2+1)k_2\pi}{2N_2}\right){\rm cos}\frac{(2n_1+1)k_1\pi}{2N_1}\\ &= {\sum}_{n_1=0}^{N_1-1}{\sum}_{n_2=0}^{N_2-1}{\rm cos}\frac{(2n_1+1)k_1\pi}{2N_1}{\rm cos}\frac{(2n_2+1)k_2\pi}{2N_2}\\ &= {\sum}_{n_2=0}^{N_2-1}\left({\sum}_{n_1=0}^{N_1-1}x_{n_1,n_2}{\rm cos}\frac{(2n_1+1)k_1\pi}{2N_1}\right){\rm cos}\frac{(2n_2+1)k_2\pi}{2N_2}\\ &= {\sum}_{n_2=0}^{N_2-1}{\sum}_{n_1=0}^{N_1-1}{\rm cos}\frac{(2n_2+1)k_2\pi}{2N_2}{\rm cos}\frac{(2n_1+1)k_1\pi}{2N_1}\\ \end{aligned} \]

记 \({\bm X}\in{\mathbb R}^{N_1\times N_2}\), \({\bm D}_1\in{\mathbb R}^{N_1\times N_1}\) 为第一维上的DCT变换, \({\bm D}_2\in{\mathbb R}^{N_2\times N_2}\) 为第二维上的DCT变换, \({\bm D}_1,{\bm D}_2\) 均由式. 生成, \({\bm Y}\in{\mathbb R}^{N_1\times N_2}\) 为二维DCT变换后的系数矩阵, 则

(3.25)¶\[{\bm Y} = {\bm D}_1{\bm X}{\bm D}_2^T \]

设 \({\bm X}\) 为一二维矩阵, 记 \({\rm DCT}({\bm X}, 1)\) 表示对 \({\bm X}\) 的第一维(列)进行DCT变换, \({\rm DCT}({\bm X}, 2)\) 表示对 \({\bm X}\) 的第二维(行)进行DCT变换, 则矩阵 \({\bm X}\) 的二维DCT变换可以表示为

(3.26)¶\[{\bm Y} = {\rm DCT2}({\bm X}) = {\rm DCT}({\rm DCT}({\bm X}, 1),2) = {\rm DCT}({\rm DCT}({\bm X}, 2),1) \]

逆二维离散余弦变换可以表示为

(3.27)¶\[{\bm X} = {\rm IDCT2}({\bm Y}) = {\rm IDCT}({\rm IDCT}({\bm Y}, 1),2) = {\rm IDCT}({\rm IDCT}({\bm Y}, 2),1) \]

3.2.4. 过完备DCT字典¶

在过完备字典(参见 Section-DictionaryType)中, 原子(列)间是线性相关的, 存在冗余原子, 原子的数目 \(N\) 大于原子的维度 \(M\). 式. 所表示的DCT基矩阵是一个完备字典, 那么如何得到过完备的DCT字典呢? 可以通过对其进行上采样和调制得到.

对于给定的信号 \({\bm y}\in {\mathbb R}^{M\times 1}\), 期望找到一组过完备基 \({\bm D} = [{\bm d}_0, {\bm d}_1,\cdots, {\bm d}_{N-1}]\) 和一组系数 \({\bm x} = [x_0, x_1, \cdots, x_{N-1}]\), 满足

\[{\bm y} = {\bm D}{\bm x} = x_0 {\bm d}_0 + x_1 {\bm d}_1 + \cdots, x_{N-1}{\bm d}_{N-1} \]

其中, \({\bm d}_i = \{1/{\sqrt 2}, \cdots, {\rm cos}\frac{(2i + 1)k\pi}{2N}, \cdots, {\rm cos}\frac{(2i + 1)(M-1)\pi}{2N}\} \in {\mathbb R}^{M\times 1}\) 为DCT基向量, \(i=0, 1, \cdots, N-1\), \(k=1, \cdots, M-1\), 过完备DCT字典 \({\bm D} \in {\mathbb R}^{M\times N}\), 可表示为

(3.28)¶\[{\bm D} = \sqrt{\frac{2}{N}}\left[\begin{array}{ccccc}{1/\sqrt{2}} & {1/\sqrt{2}} & {1/\sqrt{2}} & {\cdots} & {1/\sqrt{2}} \\ {\cos \frac{\pi}{2 N}} & {\cos \frac{3 \pi}{2 N}} & {\cos \frac{5 \pi}{2 N}} & {\cdots} & {\cos \frac{(2 N-1) \pi}{2 N}} \\ {\vdots} & {\vdots} & {\vdots} & {\vdots} \\ {\cos \frac{(M-1) \pi}{2 N}} & {\cos \frac{3(M-1) \pi}{2 N}} & {\cos \frac{5(M-1) \pi}{2 N}} & {\cdots} & {\cos \frac{(2 N-1)(M-1) \pi}{2 N}}\end{array}\right] \]

通常, 会对字典的列进行归一化, 归一化后的字典可以表示为

(3.29)¶\[{\bm D} = \left[\frac{{\bm d}_0}{\|{\bm d}_0\|_2}, \frac{{\bm d}_1}{\|{\bm d}_1\|_2},\cdots, \frac{{\bm d}_{N-1}}{\|{\bm d}_{N-1}\|_2}\right] \]

一维过完备DCT字典¶

一维过完备DCT字典的构建较为简单, 只需按照式. 生成矩阵, 再按照式. 对字典的每一列(原子)进行归一化即可.

多维过完备DCT字典¶

多维过完备DCT字典的构建以一维过完备DCT字典为基础, 设矩阵 \({\bm D}_{1d} \in {\mathbb R}^{M\times N}\) 为一维归一化过完备DCT字典, 则二维DCT字典 \({\bm D}_{2d}\) 可以通过下式计算

(3.30)¶\[{\bm D}_{2d} = {\bm D}_{1d} \times {\bm D}_{1d}, \]

其中, \(\times\) 表示 Kronecker积 (Kronecker Product). 依次类推可以得到多维过完备DCT字典

(3.31)¶\[{\bm D}_{nd} = {\bm D}_{(n-1)d} \times {\bm D}_{(n-1)d}. \]

提示

按此规则生成的字典 \({\bm D}_{nd}\) 大小为 \(M^n \times N^n\), 每一列代表一个原子, 将原子重新排成 \(n\) 维的矩阵即可.

3.2.5. 实验与分析¶

一维离散余弦变换¶

原始定义¶

Python实现

def dct1(x):
    x = np.array(x)
    N = np.size(x)

    y = np.zeros(N)

    y[0] = np.sum(x) / np.sqrt(N)

    for k in range(1, N):
        d = np.cos(
            (np.linspace(0, N, N, endpoint=False) * 2 + 1) * k * np.pi / (2.0 * N))
        y[k] = np.sqrt(2.0 / N) * np.sum(x * d)

    return y

设置序列 \(x[n] = n, n=0,1,2,\cdots,5\), 调用上述函数以及Matlab中的 dct 函数, 得实验结果

Matlab结果

>> x=[0,1,2,3,4,5]

x =

     0     1     2     3     4     5

>> dct(x)

ans =

    6.1237   -4.1626         0   -0.4082         0   -0.0801

Python结果

x = [0, 1, 2, 3, 4, 5]
y = dct(x)
print(y)

>> [ 6.12372436e+00 -4.16256180e+00 -1.53837015e-15 -4.08248290e-01 -1.53837015e-15 -8.00788912e-02]

矩阵法实现¶

Matlab实现构造DCT变换矩阵

x = [0, 1, 2, 3, 4, 5]

n = length(x)

[cc,rr] = meshgrid(0:n-1);

T = sqrt(2 / n) * cos(pi * (2*cc + 1) .* rr / (2 * n));
T(1,:) = T(1,:) / sqrt(2);

y = T * x;
disp(y)

Python实现构造DCT变换矩阵

x = [0, 1, 2, 3, 4, 5]

N = np.size(x)

r, c = np.mgrid[0: N, 0: N]

T = np.sqrt(2 / N) * np.cos(np.pi * (2 * c + 1) * r / (2 * N))
T[0, :] = T[0, :] / np.sqrt(2)

y = np.matmul(T, x)
print(y)

图像的1维DCT变换¶

图 3.14 lena 图像的1维DCT变换¶

二维离散余弦变换¶

仿真数据¶

Matlab结果

a =

     0     1     2
     3     4     5

>> dct2(a)

ans =

    6.1237   -2.0000   -0.0000
   -3.6742         0         0

def dct1(x, axis=0):
    x = np.array(x)
    if np.ndim(x) > 1:
        N = np.size(x, axis=axis)
        T = dctmat(N)
        if axis is 0:
            return np.matmul(T, x)
        if axis is 1:
            return np.matmul(T, x.transpose()).transpose()
    if np.ndim(x) is 1:
        N = np.size(x)
        T = dctmat(N)
        return np.matmul(T, x)

def dct2(X):
   return dct1(dct1(x, axis=0), axis=1)

x = [[0, 1, 2], [3, 4, 5]]
y = dct2(x)
print(y)

>> [[ 6.12372436e+00 -2.00000000e+00  5.32490308e-16]
    [-3.67423461e+00 -3.44458101e-16 -2.51493529e-16]]

真实图像数据¶

图 3.15 lena 图像二维DCT变换及IDCT变换¶

DCT过完备字典¶

Python 实现参见文件 demo_odctdict.py

生成大小分别为 \(16\times 8, 16\times 16, 16\times 32, 16\times 64\) 的一维过完备DCT字典, 基于这些字典, 根据式. 生成的二维过完备字典.

图 3.16 1-dimension overcomplete DCT dictionary. These four dictionary are generated from 1-dimension DCT dictionary with size of \(16\times 8, 16\times 16, 16\times 32, 16\times 64\).¶

依次将由式. 生成的二维过完备字典的每一列(原子) \({\bm d}_i\in {\mathbb R}^{256\times 1}\) reshape 成 \(16\times 16\) 大小的矩阵, 并按顺序从左至右从上至下排列成方阵

图 3.17 2-dimension overcomplete DCT dictionary. These four dictionary are generated from 1-dimension DCT dictionary with size of \(16\times 8, 16\times 16, 16\times 32, 16\times 64\) using formular 式. respectively.¶

Footnotes

1: http://c.csie.org/~itct/slide/DCT_larry.pdf
2: http://en.volupedia.org/wiki/Discrete_cosine_transform