数字图像处理 --- 离散余弦变换(python实战)

离散余弦变换(python实战)

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

一维DCT正变换

已知一维DCT正变换公式如下：

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

其中，x[n]为输入离散序列，n=0，1，...，N-1。X[ $\mu$ ]为经过离散余弦变换后的结果，即离散余弦变换的系数， $\mu$ =0，1，...，N-1。

根据上述公式可知每个离散余弦系数X[ $\mu$ ]都是N个输入x[n]与N个cos函数的线性组合。因此，在计算每个离散余弦系数X[ $\mu$ ]的时候都要用到N个输入x[n]。E( $\mu$ )为归一化系数， $\mu$ =0，1，...，N-1。

当 $\mu$ =0时，归一化系数E( $\mu$ )= $\frac{1}{\sqrt{2}}$ (用于直流DC分量)，此时对应的一维变换公式为：

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

当 $\mu$ 不为0时，归一化系数E( $\mu$ )=1(用于交流AC分量)，此时对应的一维变换公式为：

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

这里的归一化系数非常有用，他能够保证用于计算DCT-II系数的计算矩阵A具有正交性。这一点在后面会看到。

假设现在要求四个点的离散余弦变换，即已知输入为x[0]，x[1]，x[2]，x[3]，要求相应的四个DCT系数X[0]，X[1]，X[2]，X[3]。

对于X[0]有(N=4)：

$X[0]=\sqrt{\frac{1}{N}}\displaystyle\sum_{n=0}^{N-1}x[n]=\sqrt{\frac{1}{4}}\displaystyle\sum_{n=0}^{3}x[n]\\=\sqrt{\frac{1}{4}}x[0]+\sqrt{\frac{1}{4}}x[1]+\sqrt{\frac{1}{4}}x[2]+\sqrt{\frac{1}{4}}x[3]$

对于X[1]有(N=4)：

$X[1]=\sqrt{\frac{2}{N}}\displaystyle\sum_{n=0}^{N-1}x[n]cos\frac{(2n+1)\mu\pi}{2N}=\sqrt{\frac{1}{2}}\displaystyle\sum_{n=0}^{3}x[n]cos\frac{(2n+1)1\pi}{8}\\=\sqrt{\frac{1}{2}}x[0]cos\frac{\pi}{8}+\sqrt{\frac{1}{2}}x[1]cos\frac{3\pi}{8}+\sqrt{\frac{1}{2}}x[2]cos\frac{5\pi}{8}+\sqrt{\frac{1}{2}}x[3]cos\frac{7\pi}{8}$

对于X[2]有(N=4)：

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

注意，在计算X[2]时，四个基础余弦波的频率是计算X[1]时的四个基波频率的两倍。

对于X[3]有(N=4)：

$X[3]=\sqrt{\frac{2}{N}}\displaystyle\sum_{n=0}^{N-1}x[n]cos\frac{(2n+1)\mu\pi}{2N}=\sqrt{\frac{1}{2}}\displaystyle\sum_{n=0}^{3}x[n]cos\frac{(2n+1)3\pi}{8}\\=\sqrt{\frac{1}{2}}x[0]cos\frac{3\pi}{8}+\sqrt{\frac{1}{2}}x[1]cos\frac{9\pi}{8}+\sqrt{\frac{1}{2}}x[2]cos\frac{15\pi}{8}+\sqrt{\frac{1}{2}}x[3]cos\frac{21\pi}{8}$

继续，在计算X[3]时，四个基础余弦波的频率是计算X[1]时的四个基波频率的三倍。

如果把X[]和x[]分别用向量来表示，则上面的四个等式可用矩阵计算表示为：

$\begin{bmatrix} X[0]\\ X[1]\\ X[2]\\ X[3] \end{bmatrix}=\begin{bmatrix} \sqrt{\frac{1}{4}}&\sqrt{\frac{1}{4}} &\sqrt{\frac{1}{4}} &\sqrt{\frac{1}{4}} \\ \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} \\\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} \\\sqrt{\frac{1}{2}}cos\frac{3\pi}{8}&\sqrt{\frac{1}{2}}cos\frac{9\pi}{8} &\sqrt{\frac{1}{2}}cos\frac{15\pi}{8} &\sqrt{\frac{1}{2}}cos\frac{21\pi}{8} \end{bmatrix}\begin{bmatrix} x[0]\\ x[1]\\ x[2]\\ x[3] \end{bmatrix}$

更广义的讲，对任意长度输入的N都可以用矩阵的形式计算系数向量：

$\begin{bmatrix} X[0]\\ X[1]\\ .\\.\\.\\ X[N-1] \end{bmatrix}=\begin{bmatrix} E(0)\sqrt{\frac{2}{N}}x[0]cos\frac{(2*0+1)0*\pi }{2N} & E(0)\sqrt{\frac{2}{N}}x[1]cos\frac{(2*1+1)0*\pi }{2N}&... &E(0)\sqrt{\frac{2}{N}}x[N-1]cos\frac{(2*(N-1)+1)0*\pi }{2N} \\ E(1)\sqrt{\frac{2}{N}}x[0]cos\frac{(2*0+1)1*\pi }{2N} & E(1)\sqrt{\frac{2}{N}}x[1]cos\frac{(2*1+1)1*\pi }{2N}&... &E(1)\sqrt{\frac{2}{N}}x[N-1]cos\frac{(2*(N-1)+1)1*\pi }{2N} \\ .&.& .&. \\ E(N-1)\sqrt{\frac{2}{N}}x[0]cos\frac{(2*0+1)(N-1)*\pi }{2N} & E(N-1)\sqrt{\frac{2}{N}}x[1]cos\frac{(2*1+1)(N-1)*\pi }{2N}&... &E(N-1)\sqrt{\frac{2}{N}}x[N-1]cos\frac{(2*(N-1)+1)(N-1)*\pi }{2N} \end{bmatrix}\begin{bmatrix} x[0]\\ x[1]\\ .\\.\\.\\ x[N-1] \end{bmatrix}$

基于上述结论可得到如下一维离散余弦变换的python代码：

import numpy as np
import math
from scipy.fft import dct,idct

def DCT1d_matrix(N):
    A=np.zeros((N,N))
    for miu in range(N):
        for x in range(N):
            if miu==0:
                E_miu=np.sqrt(1/2)
            else:
                E_miu=1
            A[miu,x]=E_miu*np.sqrt(2/N)*np.cos((2*x+1)*miu*np.pi/2/N)
    return A

def DCT1D(x):
    N=len(x)
    A=DCT1d_matrix(N)
    X=A@x
    return X

令N=4，可得到上面用于计算4点DCT的DCT变换矩阵A。

N=4
A=DCT1d_matrix(N)
print("DCT Transform Matrix for N = 4:")
print(A)

可能你得到的DCT变换矩阵跟我的不同，很有可能是因为你在计算DCT的时候没有代入归一化系数。

为了验证我生成的矩阵是否正确，下面我用这个矩阵计算一下四个点输入的DCT，并和scipy自带库的计算结果比较。需要注意的是，scipy库中的dct专门有一个控制norm的参数，如果这个参数不开，计算结果就和我的结果不同。只有开启了norm的参数，才和我的结果一样。

np.random.seed(123)
x=np.random.rand(4)
X=DCT1D(x)
print("Input signal:", x)
print("My DCT-II coefficients:",X)

#check result
X_check=dct(x,type=2)#不带归一化系数
X_nrom=dct(x,type=2,norm='ortho')#带归一化系数
print("\nResults are calced by lib:")
print("DCT-II coefficients(with out norm):",X_check)
print("DCT-II coefficients(with norm):",X_nrom)

可见，我根据一维dct正变换的计算公式生成变换矩阵A的“DCT1d_matrix(N)”函数没问题，且，通过矩阵计算的方式计算dct系数的“DCT1D(x)”函数也没问题。

一维DCT正变换公式的意义：变换矩阵A

回到公式，如果暂时不考虑归一化系数，只是单看红色方框中的部分。从信号的角度讲，实际上就是把输入信号x[n]与不同的频率的余弦信号做相关，准确的说应该是互相关。互相关常被用来度量两个信号之间相似性。

理论上，两个信号在互相关的时候一个信号动，另一个信号不动。移动的信号每移动一次得到一个值。只不过这里的互相关是正好当两个信号移动到完全重合的时候得到的值，即，两个信号逐一相乘的结果。

现在，我把上式展开，“互相关”这一点就能看的更清楚。

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

如果分别用向量x和c表示输入信号x和与之相关的信号cos函数：

$x=\begin{bmatrix} x[0]\\ x[1]\\ ...\\ x[N-1] \end{bmatrix}$ $c_{\mu }=\begin{bmatrix} cos\frac{(2*0+1)\mu\pi}{2N}\\cos\frac{(2*1+1)\mu\pi}{2N} \\... \\ cos\frac{(2*(N-1)+1)\mu\pi}{2N}\end{bmatrix}$

则上式可改写为：

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

对于不同的 $\mu$ =0，1，...，N-1有：

$X[0]=E(0)\sqrt{\frac{2}{N}}(x\cdot c_{0})$

$X[1]=E(1)\sqrt{\frac{2}{N}}(x\cdot c_{1})$

。。。

$X[N-1]=E(N-1)\sqrt{\frac{2}{N}}(x\cdot c_{N-1})$

这就是说，离散余弦变换实际上是用输入信号x和不同频率的cos信号做互相关，每个频率所对应的dct系数，实际上就是互相关的结果。更准确的说，就是两个信号的乘积。而这些不同频率的cos信号，实际上就是dct变换矩阵中的每一行，cos信号的长度和输入x相同。

如果我们把DCT公式中的cos函数展开后，并用熟知的幅值+频率+相位的公式来表示的话：

$y=Acos(B(x+C))+D$

$E(\mu)\sqrt{\frac{2}{N}}cos\frac{(2n+1)\mu\pi}{2N}=E(\mu)\sqrt{\frac{2}{N}}cos(\frac{\mu \pi }{N}(n+1/2))$

其中：

$A=E(\mu)\sqrt{\frac{2}{N}}$

$B=\frac{\mu \pi }{N}$

$C=\frac{1}{2}$

$D=0$

这就是说，实际上离散余弦变换就是用上面的这些幅值为A，周期为B，相位为C的余弦信号与输入信号做相关后得到变换结果的(也就是频域中不同频率cos函数的权重)。下面我们分别以N=4，N=8和N=16为例子，看看这些余弦信号究竟长什么样。这里再次重复一下，这些余弦信号就是变换矩阵A中的每一个行向量。

4点DCT的基础余弦信号：

N=4
A=DCT1d_matrix(N)
print(A)
fig,axs=plt.subplots(N//2,2)
plt.subplots_adjust(hspace=0.5)
for u in range(N):
    row=u//2
    col=u%2
    axs[row, col].set_ylim((-1, 1))
    axs[row, col].set_title(f"cos wave with freq of {u}*pi/{N}")
    axs[row, col].plot(A[u, :])
    axs[row, col].plot(A[u, :],'ro')

注意：图像中的数值与变换矩阵A中每行的数值一一对应！

8点DCT的基础余弦信号：

N=8
A=DCT1d_matrix(N)
print(A)
fig,axs=plt.subplots(N//2,2,figsize=(6,10))
plt.subplots_adjust(hspace=0.5)
for u in range(N):
    row=u//2
    col=u%2
    axs[row, col].set_ylim((-1, 1))
    axs[row, col].set_title(f"cos wave with freq of {u}*pi/{N}")
    axs[row, col].plot(A[u, :])
    axs[row, col].plot(A[u, :],'ro')

注意：图像中的数值与变换矩阵A中每行的数值一一对应！

16点DCT的基础余弦信号：

N=16
A=DCT1d_matrix(N)
print(A)
fig,axs=plt.subplots(N//2,2,figsize=(7,15))
plt.subplots_adjust(hspace=0.5)
for u in range(N):
    row=u//2
    col=u%2
    axs[row, col].plot(A[u, :])
    axs[row, col].plot(A[u, :],'ro')
    axs[row, col].set_ylim((-1, 1))
    axs[row, col].set_title(f"cos wave with freq of {u}*pi/{N}")

注意：图像中的数值与变换矩阵A中每行的数值一一对应！

如果比较不同长度dct的前四个基础cos函数可知：

1，第一个也是最容易被观察到的就是，随着N的增加，需要越来越多的高频cos函数合成原始信号x。

2，就前四个共有的基础cos函数而言，随着N的增加，4个cos函数的频率都是相同的，样本数越来越稠密，但由于归一化系数的存在幅值越来越小。

小结：

简而言之，所谓的DCT离散余弦变换就是要用不同频率的cos函数去合成输入信号。而每个cos函数的幅值，是通过cos函数与输入信号的互相关/乘积/点积决定的。互相关的结果的和就是输入信号DCT的结果，即，DCT系数。

以本文开始处4个点的随机数的DCT为例：

np.random.seed(123)
x=np.random.rand(4)
print("Input signal x:",x)
X=DCT1D(x)
print("DCT of x:",X)

A=DCT1d_matrix(4)
print("DCT transfrom matrix:\n",A)

cos1=A[1,:]
ccr1=x*cos1
print("Result of cross correlation:",ccr1)
print("DCT coef=sum of ccr=",sum(ccr1))

fig,axs=plt.subplots(3,1,figsize=(5,7))
plt.subplots_adjust(hspace=0.5)
axs[0].plot(x)
axs[0].plot(x,'ro')
axs[0].set_title('input signal x')

axs[1].plot(cos1)
axs[1].plot(cos1,'ro')
axs[1].set_title('$cos_1(x)$')

axs[2].plot(ccr1)
axs[2].plot(ccr1,'ro')
axs[2].set_title('DCT coef of x')

下图展示了输入信号x与变换矩阵A中第二个cos函数互相关的结果：

DCT系数等于互相关结果的和，这一点和DCT1D函数计算的结果一致：

一维DCT反变换

前面提到过1维DCT正变换函数带有归一化系数，是因为这使得DCT正变换矩阵A是一个标准正交矩阵。又因为，DCT正变换可以通过矩阵计算得到结果。

因此，我在计算DCT反变换时，可以抛弃DCT反变换的计算公式，直接通过DCT正变换矩阵A的逆求得原始输入x。

Tips: 若矩阵Q为标准正交矩阵，则有

$Q \cdot Q^T=Q^T \cdot Q=I$

这说明，正交矩阵Q的逆矩阵是该矩阵的转置，即：

$Q^{-1}=Q^T$

因此，计算DCT反变换时只需用DCT正变换矩阵A的转置乘以系数X即可：

$x[n]=A^TX[\mu]$

np.random.seed(123)
x=np.random.rand(4)
print("Input signal x:",x)
X=DCT1D(x)
print("\nDCT of x:",X)

A=DCT1d_matrix(4)
AT=A.T
print("\nForward transform matrix A:\n",A)
print("\nInverse transform matrix AT:\n",AT)
x_inv=AT@X
print("\nIDCT of X:",x_inv)

这里稍微提醒一下，在使用上述方法进行DCT反变换之前。最好先验证一下DCT正变换矩阵A是否是标准正交矩阵。(以N=4为例)

1，看看每个向量的长度是否为1，即，单位向量。

2，看看两两向量之间是否彼此正交。

3，看看矩阵A的转置是否为A的逆。

#检查所有行向量是否为单位向量
print("row check:")
for i in range(4):
    print(A[i,:]@A[i,:].T)

#检查所有行向量是否彼此正交
for i in range(4):
    for j in range(i+1,4):
        print(A[i,:]@A[j,:].T)
        
    
#检查所有列向量是否为单位向量
print("\ncol check:")
for i in range(4):
    print(A[:,i]@A[:,i].T)
    

#检查所有列向量是否彼此正交
for i in range(4):
    for j in range(i+1,4):
        print(A[:,i]@A[:,j].T)
        
#检查矩阵A的转置是否为他的逆矩阵
print("\nA*AT:\n",A@A.T)
print("\nAT*A:\n",A.T@A)