数字信号处理 --- 二维离散余弦变换(python实战)

松下J27

已于 2024-08-23 22:29:47 修改

阅读量1.3k

点赞数 30

于 2024-07-30 22:12:05 首次发布

本文链接：https://blog.csdn.net/daduzimama/article/details/140780627

版权

Digital Image Processing 同时被 2 个专栏收录

13 篇文章 151 订阅

订阅专栏

Digital Signal Processing

11 篇文章 55 订阅

订阅专栏

二维离散余弦变换(python实战)

（注：文中所讨论的离散余弦变换都是DCT-II）

在上一篇文章中，我详细介绍了一维离散余弦变换的基本原理，并且以可视化的方式展示了一维DCT中用于分析输入信号的一系列基础余弦波。在这篇文章中继续关于DCT的学习，并把DCT从一维推广到二维，最后以图像的离散余弦变换为结尾。

数字图像处理 --- 一维离散余弦变换(python实战)-CSDN博客文章浏览阅读303次，点赞6次，收藏9次。本文详细介绍了基于python实现的一维离散余弦变换，是我个人的学习笔记。https://blog.csdn.net/daduzimama/article/details/140449135

二维DCT正变换

二维DCT正变换公式如下：

$X[\mu,\nu ]=E(\mu)E(\nu)\frac{2}{N}\displaystyle\sum_{m=0}^{N-1}\displaystyle\sum_{n=0}^{N-1}x[m,n]cos\frac{(2m+1)\mu\pi}{2N}cos\frac{(2n+1)\nu\pi}{2N}$

或

$X[\mu,\nu ]=E(\mu)\sqrt{\frac{2}{N}}\displaystyle\sum_{m=0}^{N-1}cos\frac{(2m+1)\mu\pi}{2N}(E(\nu)\sqrt{\frac{2}{N}}\displaystyle\sum_{n=0}^{N-1}x[m,n]cos\frac{(2n+1)\nu\pi}{2N})$

其中，x[m,n]为输入图像的像素值，m,n=0，1，...，N-1。X[ $\mu ,\nu$ ]为经过离散余弦变换后的结果，即离散余弦变换的系数， $\mu$ , $\nu$ =0，1，...，N-1。E( $\mu$ )，E( $\nu$ )为归一化系数，当 $\mu ,\nu$ =0时， E( $\mu$ ),E( $\nu$ )= $\frac{1}{\sqrt{2}}$ 。当 $\mu$ , $\nu$ =1，2，3,....,N-1时， E( $\mu$ ),E( $\nu$ )=1。

单从公式来看，我们可以知道三点。第一，频域中的每一个点X[ $\mu ,\nu$ ]都要遍历整幅图像计算。二，就某个固定的 $\mu$ , $\nu$ 而言，分析每行和每列数据用的是同一组cos函数。第三，越是往频域的右下角计算，分析图像所需的cos函数的频率就越高。

以4个点的X[0,0]的计算为例，N=4， $\mu ,\nu$ =0：

$X[0,0]=E(0)\sqrt{\frac{2}{N}}\displaystyle\sum_{m=0}^{N-1}cos\frac{(2m+1)0\pi}{2N}(E(0)\sqrt{\frac{2}{N}}\displaystyle\sum_{n=0}^{N-1}x[m,n]cos\frac{(2n+1)0\pi}{2N})$

$=1/\sqrt{N}\displaystyle\sum_{m=0}^{N-1}cos\frac{(2m+1)0\pi}{2N}(1/\sqrt{N}\displaystyle\sum_{n=0}^{N-1}x[m,n]cos\frac{(2n+1)0\pi}{2N})$

$=1/2\displaystyle\sum_{m=0}^{N-1}(1/2\displaystyle\sum_{n=0}^{N-1}x[m,n])$

$=1/2\displaystyle\sum_{m=0}^{N-1}(1/2x[m,0]+1/2x[m,1]+1/2x[m,2]+1/2x[m,3])$

$=(1/2(1/2x[0,0]+1/2x[0,1]+1/2x[0,2]+1/2x[0,3]))+(1/2(1/2x[1,0]+1/2x[1,1]+1/2x[1,2]+1/2x[1,3]))+(1/2(1/2x[2,0]+1/2x[2,1]+1/2x[2,2]+1/2x[2,3]))+(1/2(1/2x[3,0]+1/2x[3,1]+1/2x[3,2]+1/2x[3,3]))$

这个公式表明，要想求出X[0,0]。要先让原始二维信号x[m,n]的每一行乘以幅值为1/2频率为0的余弦信号[0.5,0.5,0.5,0.5]，然后再用这个结果乘以余弦信号[0.5,0.5,0.5,0.5]才能得到最终结果。下面是基于python的实现，主要是用于验证计算结果并演示计算过程：

#生成随机二维信号
x2=np.random.randn(4,4)
x2

导入opencv的dct库，计算二维DCT。后面会根据我上面推导出来的公式，计算X[0,0]，并用CV库的计算结果来验证。

import cv2
dct_matrix=cv2.dct(x2)
dct_matrix

以下是基于上面的推导结果求得的X[0，0]，整个计算过程是先用幅值为1/2频率为0的cos函数乘以二维数组的每一行并得到一个中间结果，最后用这个中间结果再次乘以cos函数得到最终结果，这一计算过程和推导一致：

A=DCT1d_matrix(4)
print("Forward transform matrix A:\n",A)

#原始信号的每一行乘以余弦信号，得到一个中间结果
step1=x2@A[0,:]
print("result of step1:",step1)

#中间结果再度乘以余弦信号，得到最终结果
step2=step1@A[0,:]
print("result of step2:",step2)

根据计算结果来看，1，计算时所使用的cos函数为一维DCT矩阵的第一行。毕竟DCT矩阵中的每一行都是一个固定频率的cos函数。2，按照这种方法计算出来的最终结果0.8999和用cv库算出来的左上角第一个值一致。

0,0点的例子比较特殊，现在我们再看一个点X[1,2]。还是4个点的DCT，N=4， $\mu =1$ ， $\nu =2$ ：

$X[1,2]=E(1)\sqrt{\frac{2}{N}}\displaystyle\sum_{m=0}^{N-1}cos\frac{(2m+1)\pi}{2N}(E(2)\sqrt{\frac{2}{N}}\displaystyle\sum_{n=0}^{N-1}x[m,n]cos\frac{(2n+1)2\pi}{2N})$

$=\sqrt{\frac{1}{2}}\displaystyle\sum_{m=0}^{N-1}cos\frac{(2m+1)\pi}{8}(\sqrt{\frac{1}{2}}\displaystyle\sum_{n=0}^{N-1}x[m,n]cos\frac{(2n+1)2\pi}{8})$

$=\sqrt{\frac{1}{2}}\displaystyle\sum_{m=0}^{N-1}cos\frac{(2m+1)\pi}{8}(\sqrt{\frac{1}{2}}x[m,0]cos\frac{2\pi}{8}+\sqrt{\frac{1}{2}}x[m,1]cos\frac{6\pi}{8}+\sqrt{\frac{1}{2}}x[m,2]cos\frac{10\pi}{8}+\sqrt{\frac{1}{2}}x[m,3]cos\frac{14\pi}{8})$

$=(\sqrt{\frac{1}{2}}cos\frac{\pi}{8}(\sqrt{\frac{1}{2}}x[0,0]cos\frac{2\pi}{8}+\sqrt{\frac{1}{2}}x[0,1]cos\frac{6\pi}{8}+\sqrt{\frac{1}{2}}x[0,2]cos\frac{10\pi}{8}+\sqrt{\frac{1}{2}}x[0,3]cos\frac{14\pi}{8}))+(\sqrt{\frac{1}{2}}cos\frac{3\pi}{8}(\sqrt{\frac{1}{2}}x[1,0]cos\frac{2\pi}{8}+\sqrt{\frac{1}{2}}x[1,1]cos\frac{6\pi}{8}+\sqrt{\frac{1}{2}}x[1,2]cos\frac{10\pi}{8}+\sqrt{\frac{1}{2}}x[1,3]cos\frac{14\pi}{8}))+(\sqrt{\frac{1}{2}}cos\frac{5\pi}{8}(\sqrt{\frac{1}{2}}x[2,0]cos\frac{2\pi}{8}+\sqrt{\frac{1}{2}}x[2,1]cos\frac{6\pi}{8}+\sqrt{\frac{1}{2}}x[2,2]cos\frac{10\pi}{8}+\sqrt{\frac{1}{2}}x[2,3]cos\frac{14\pi}{8}))+(\sqrt{\frac{1}{2}}cos\frac{7\pi}{8}(\sqrt{\frac{1}{2}}x[3,0]cos\frac{2\pi}{8}+\sqrt{\frac{1}{2}}x[3,1]cos\frac{6\pi}{8}+\sqrt{\frac{1}{2}}x[3,2]cos\frac{10\pi}{8}+\sqrt{\frac{1}{2}}x[3,3]cos\frac{14\pi}{8}))$

结合推导结果来看，要想得到X[1,2]的结果。先是要让原始二维信号的每一行都乘以频率为2/4pi的cos函数得到中间结果，这一步所使用的cos函数为：

$[\sqrt{\frac{1}{2}}cos\frac{2\pi}{8},\sqrt{\frac{1}{2}}cos\frac{6\pi}{8},\sqrt{\frac{1}{2}}cos\frac{10\pi}{8},\sqrt{\frac{1}{2}}cos\frac{14\pi}{8}]$

此时 $\nu =2$ ，然后再用这个中间结果乘以频率为1/4pi的cos函数得到最终结果。这一步所使用的cos函数为：

$[\sqrt{\frac{1}{2}}cos\frac{\pi}{8},\sqrt{\frac{1}{2}}cos\frac{3\pi}{8},\sqrt{\frac{1}{2}}cos\frac{5\pi}{8},\sqrt{\frac{1}{2}}cos\frac{7\pi}{8}]$

此时 $\mu =1$ ，实际上这两组cos函数就是一维4点DCT变换矩阵其中的两行。

cos0=np.array([np.sqrt(1/2)*np.cos(2*np.pi/8),np.sqrt(1/2)*np.cos(6*np.pi/8),np.sqrt(1/2)*np.cos(10*np.pi/8),np.sqrt(1/2)*np.cos(14*np.pi/8)])
print("cos with freq of 2/4pi:",cos0)

cos1=np.array([np.sqrt(1/2)*np.cos(1*np.pi/8),np.sqrt(1/2)*np.cos(3*np.pi/8),np.sqrt(1/2)*np.cos(5*np.pi/8),np.sqrt(1/2)*np.cos(7*np.pi/8)])
print("cos with freq of 1/4pi:",cos1)

最后再次基于推导结果计算X[1,2]，先用 $\nu =2$ 时的cos函数乘以二维信号的每一行，得到中间结果。然后再用 $\mu =1$ 时的cos函数乘以中间结果，得到最终结果：

#原始信号的每一行乘以nu=2时的余弦信号，得到一个中间结果
step1=x2@A[2,:]
print("result of step1:",step1)

#中间结果再乘以mu=1时的余弦信号，得到最终结果
step2=step1@A[1,:]
print("result of step2:",step2)

二维DCT正变换的矩阵计算

在上面的两个例子中，我分别列举了频域中两个点的计算，事实上，和一维DCT一样，二维信号所有点的DCT可以通过矩阵运算的方式一次性求出来，具体来说就是先对二维信号的每一行进行一维DCT，得到一个中间结果矩阵。然后对中间结果的每一列进行一维DCT，得到最终的二维DCT结果。

1，首先，一维DCT的矩阵计算公式为：

$X[\mu ]=Ax[n]$

2, 基于这一公式，我们可以先用 $Ax[m,n]^{T}$ 实现对x的每行数据进行一维DCT变换。其中， $Ax[m,n]^{T}$ 可拆分成下图的计算方式，其中每个Acol是矩阵中的一列。

$Ax^{T}=\begin{bmatrix} A\cdot col0 &A\cdot col1 & ... &A\cdot colN \end{bmatrix}$

上面公式中的每一列 $A\cdot col$ 又是按照下图的方式计算的：

$\begin{bmatrix} .& . & .\\ . &A &. \\ . & . & . \end{bmatrix}\begin{bmatrix} .\\ col\\ . \end{bmatrix}=\begin{bmatrix} .\\ .\\ . \end{bmatrix}$

上图的意思是矩阵A的每一行都是一组cos基础函数，每行cos函数乘以col中的数据就是col的DCT系数。又因为col中所装的都来自于 $x^{T}$ ，即x的转置。因此，矩阵A所乘的每一列col就是原始二维数据中的每一行。A乘col的结果是一个列向量。

3，上一步中 $Ax[m,n]^{T}$ 的结果是一个与x尺寸相同的矩阵，是中间结果。该矩阵中的每一列都是原始二维矩阵x每行一维DCT的结果。首先，我们要对中间结果转置，即， $(Ax[m,n]^{T})^{T}$ ，转置后矩阵中每一行都是x每一行DCT的结果。然后再对转置后的中间结果的每一列做一维DCT，直接套用一维DCT的矩阵计算公式即可，得到：

$A(Ax[m,n]^{T})^{T}=Ax[m,n]A^{T}$

最终得到二维DCT的矩阵运算公式：

$X[\mu ,\nu ]=Ax[m,n]A^{T}$

用python代码实现如下：

X2=(A@x2)@A.T
X2

这一计算结果与之前用cv2库计算所得的结果一致！

二维DCT反变换

已知二维DCT正变换的矩阵计算公式为：

$X[\mu ,\nu ]=Ax[m,n]A^{T}$

现在要求二维DCT反变换，即已知X，要求x，需要分两步计算：

1，等式两边的前面同时乘以A^T，得到：

$A^TX[\mu,\nu]=A^TAx[x,y]A^T \Rightarrow A^TX[\mu,\nu]=x[x,y]A^T$

2，等式两边的后面同时乘以A：

$A^TX[\mu,\nu]A=x[x,y]A^TA\Rightarrow A^TX[\mu,\nu]A=x[x,y]$

最终得到二维DCT反变换的公式：

$x[x,y]=A^TX[\mu,\nu]A$

print("Original data:\n",x2)
print("\nForward 2d DCT:\n",X2)
x2_inv=A.T@X2@A
print("\nInv 2d DCT:\n",x2_inv)

用相应的库函数验证DCT结果：

idct_matrix = idct(idct(dct_matrix.T, norm='ortho').T, norm='ortho')
idct_matrix

（全文完）

作者 --- 松下J27

参考文献：

1，《数字图像处理技术详解与Visual C++实践》---左飞

2，dct — SciPy v1.14.0 Manual

3，dct

4，https://en.wikipedia.org/wiki/Discrete_cosine_transform

5，https://en.wikipedia.org/wiki/Cross-correlation

6，Amplitude, Period, Phase Shift and Frequency

（配图与本文无关）