A Fast Computational Algorithm for the Discrete Cosine Transform

【论文地址30027-5/h0085)】

离散余弦变换的一种快速计算算法

摘要

本文提出了一种快速离散余弦变换算法，与传统的使用快速傅里叶变换的离散余弦变换算法相比，计算复杂度改善了6倍，该算法以矩阵的形式导出，并用信号流图来说明，信号流图可以很容易地转换为硬件或软件实现。

引言

离散余弦变换（DCT）已经被成功应用到高分辨率图像编码中。传统的实现DCT变换的方法是利用双倍大小的快速傅里叶（FFT）算法，在计算过程中使用了复数计算。由于缺乏有效的算法，DCT在各种各样的应用中并没有像它的特性所暗示的那样广泛。本文描述了一种只涉及实数操作的更有效的计算一组N个点的快速离散余弦变换(FDCT)算法。该算法可以扩展到任意期望的值 $N=2^m,m>=2$。这种推广包括交替的余弦/正弦蝴蝶矩阵和二进制矩阵，以重新排列矩阵元素，使其在每个其他节点上保持可识别的位反转模式。泛化并不是唯一的——已经发现了几种替代方法，但这里描述的方法似乎是最简单的解释。它不一定是最有效的FDCT，但代表了一种方法扩展的技术。这种方法采用了 $(3N/2)(log_2N-1)+2$个实数加法和 $Nlog_2N-3N/2+4$​个实数乘法，这大约是使用双倍大小FFT的传统方法的6倍快。

离散余弦变换

离散函数 $f(j)，j=0,1,…,N-1$的离散余弦变换定义如下：

$$F(k)=\frac{2c(k)}{N}\sum^{N-1}_{j=0}f(j)cos[\frac{(2j+1)k\pi}{2N}],k=0,1,...,N-1$$

逆变换定义如下：

$$f(j)=\sum^{N-1}_{k=0}c(k)F(k)cos[\frac{(2j+1)k\pi}{2N}],j=0,1,...,N-1$$

其中，

$$c(k)=\frac{1}{\sqrt{2}}, k=0\\ c(k)=1,k=1,2,...,N-1$$

该变换比任何已知的快速计算算法拥有更高的能量压缩特性。该变换还具有循环卷积-乘法关系，可方便地应用于线性系统理论中。

快速计算算法

一个N*1的数据向量[f]的离散余弦变换可以用矩阵形式表达如下：

$$[F]=\frac{2}{N}[A_N][f]$$

其中，$[A_N]=[c(k)cos(2j+1)k\pi/2N];j,k=0,1,…,(N-1)$​​,[F]是N*1的转换后的向量。本文提出的快速计算算法是基于$[A_N]$矩阵的矩阵分解。如下所示，该矩阵可以写成下面的递归形式：

$$[A_N]=[P_N]\left[\begin{array}{cc}A_{N/2}&0\\0&R_{N/2}\end{array}\right][B_N]$$

其中，$[B_N]$定义如下：

$$[R_{N/2}]=[c(k)cos\frac{(2j+1)(2k+1)\pi}{2N}];j,k=0,1,...,\frac{N}{2}-1$$

$[P_N]$​是一个N x N的排列矩阵，它将变换后的向量从位逆序排列为一个自然的顺序。2*2DCT变换可以写成如下形式

$$[A_2]=\frac{1}{\sqrt(2)}\left[\begin{array}{cc}1&1\\1&-1\end{array}\right]$$

从 $[A_N]$的递归性质可以看出，只要找到一种通用的 $[R_{N/2}]$矩阵分解方法，$A_2$就可以扩展到更高阶的矩阵。

下面的讨论提出了一种分解 $[R_{N/2}]$ 矩阵的系统方法,需要强调的是这种分解方法不是唯一的，也不是最优的。另外有几种方法需要较少的计算步骤，但它们对更大的尺寸没有明显的泛化。

$[R_{N/2}]$矩阵可以用下面的方式分解成（$2log_2N-3$）个矩阵：

$$[R_{N/2}]=[M1][M2][M3][M4]...[M(2log_2N-3)]$$

矩阵有四种不同的类型

类型1：[M1],第一个矩阵

类型2：$[M(2log_2N-3)]$，最后一个矩阵

类型3：[Mq],奇数矩阵，如[M3],[M5]等

类型4：[Mp],偶数矩阵，如[M2],[M4]等

在详细描述这四种矩阵之前，给出了以下定义:

$$B_N=\left[\begin{array}{cc}I_{N/2}&\overline{I}_{N/2}\\\overline{I}_{N/2} &-I_{N/2}\end{array}\right]$$ $$B_N^*=\left[\begin{array}{cc}-I_{N/2}&\overline{I}_{N/2}\\\overline{I}_{N/2} &I_{N/2}\end{array}\right]$$

$I_{N/2}$是阶数为 $\frac{N}{2}$的单位矩阵

$\overline{I}_{N/2}$是反对角单位矩阵

$[S_i^k]=sin\frac{k\pi}{i}[I_{N/2i}]$

$[\overline{s}_i^k]=sin\frac{k\pi}{i}[\overline{I}_{N/2i}]$

$[C_i^k]=cos\frac{k\pi}{i}[I_{N/2i}]$​

$[\overline{C}_i^k]=cos\frac{k\pi}{i}[\overline{I}_{N/2i}]$​

上式中，单位矩阵 $[I_{N/2i}]$表示 对角sin或者cos矩阵 $[S_i^k],[\overline{S}_i^k],[C_i^k],[\overline{C}_i^k]$的阶数，另外i>N/2 时，$[I_{N/2i}]\equiv1$。

实验

实验代码

import numpy as np
import cv2

#传入要计算的方阵
def cal_dct(X):
	N=X.shape[0]
	A=np.zeros((N,N))
	# print(X)
	for i in range(N):
		for j in range(N):
			if i == 0:
				a=np.sqrt(1/N)
			else:
				a=np.sqrt(2/N)
			A[i,j]=a*np.cos(np.pi*(j+0.5)*i/N)
	# print(A)
	#DCT变换
	#dot矩阵乘法，*号是对应元素相乘
	Y=A.dot(X).dot(A.T)
	return Y


def cal_idct(X):
	N=X.shape[0]
	A=np.zeros((N,N))
	for i in range(N):
		for j in range(N):
			if i == 0:
				a=np.sqrt(1/N)
			else:
				a=np.sqrt(2/N)
			A[i,j]=a*np.cos(np.pi*(j+0.5)*i/N)
	#DCT变换
	#dot矩阵乘法，*号是对应元素相乘
	
	Y=(A.T).dot(X).dot(A)
	return Y


#快速DCT计算算法
def fast_DCT(X):
	#计算C8
	P8=np.array([[1,0,0,0,0,0,0,0],[0,0,0,0,1,0,0,0],[0,0,1,0,0,0,0,0],[0,0,0,0,0,0,1,0],[0,1,0,0,0,0,0,0],[0,0,0,0,0,1,0,0],[0,0,0,1,0,0,0,0],[0,0,0,0,0,0,0,1]])
	P4=np.array([[1,0,0,0],[0,0,1,0],[0,1,0,0],[0,0,0,1]])

	M1=np.array([[np.sin(np.pi/16),0,0,np.cos(np.pi/16)],[0,np.sin(5*np.pi/16),np.cos(5*np.pi/16),0],[0,-np.sin(3*np.pi/16),np.cos(3*np.pi/16),0],[-np.sin(7*np.pi/16),0,0,np.cos(7*np.pi/16)]])
	M2=np.sqrt(2)/2*np.array([[1,1,0,0],[1,-1,0,0],[0,0,-1,1],[0,0,1,1]])
	M3=np.array([[1,0,0,0],[0,-np.cos(np.pi/4),np.cos(np.pi/4),0],[0,np.cos(np.pi/4),np.cos(np.pi/4),0],[0,0,0,1]])

	# print(M3)
	R4=M1.dot(M2).dot(M3)
	B8=np.sqrt(2)/2*np.array([[1,0,0,0,0,0,0,1],[0,1,0,0,0,0,1,0],[0,0,1,0,0,1,0,0],[0,0,0,1,1,0,0,0],[0,0,0,1,-1,0,0,0],[0,0,1,0,0,-1,0,0],[0,1,0,0,0,0,-1,0],[1,0,0,0,0,0,0,-1]])

	#计算C4
	P2=np.array([[1,0],[0,1]])
	C2=np.sqrt(2)/2*np.array([[1,1],[1,-1]])
	R2=np.array([[np.sin(np.pi/8),np.cos(np.pi/8)],[-np.sin(3*np.pi/8),np.cos(3*np.pi/8)]])
	B4=np.sqrt(2)/2*np.array([[1,0,0,1],[0,1,1,0],[0,1,-1,0],[1,0,0,-1]])
	
	tmp1=np.concatenate((P2.T.dot(C2),np.zeros((2,2))),axis=1)
	tmp2=np.concatenate((np.zeros((2,2)),R2),axis=1)
	tmp=np.concatenate((tmp1,tmp2))
	C4=P4.dot(tmp).dot(B4)

	tmp1=np.concatenate((P4.T.dot(C4),np.zeros((4,4))),axis=1)
	tmp2=np.concatenate((np.zeros((4,4)),R4),axis=1)
	tmp=np.concatenate((tmp1,tmp2))
	#得到正交矩阵C8
	C8=P8.dot(tmp).dot(B8)

	Y=C8.dot(X).dot(C8.T)
	return Y

def fast_IDCT(X):
	#计算C8
	P8=np.array([[1,0,0,0,0,0,0,0],[0,0,0,0,1,0,0,0],[0,0,1,0,0,0,0,0],[0,0,0,0,0,0,1,0],[0,1,0,0,0,0,0,0],[0,0,0,0,0,1,0,0],[0,0,0,1,0,0,0,0],[0,0,0,0,0,0,0,1]])
	P4=np.array([[1,0,0,0],[0,0,1,0],[0,1,0,0],[0,0,0,1]])

	M1=np.array([[np.sin(np.pi/16),0,0,np.cos(np.pi/16)],[0,np.sin(5*np.pi/16),np.cos(5*np.pi/16),0],[0,-np.sin(3*np.pi/16),np.cos(3*np.pi/16),0],[-np.sin(7*np.pi/16),0,0,np.cos(7*np.pi/16)]])
	M2=np.sqrt(2)/2*np.array([[1,1,0,0],[1,-1,0,0],[0,0,-1,1],[0,0,1,1]])
	M3=np.array([[1,0,0,0],[0,-np.cos(np.pi/4),np.cos(np.pi/4),0],[0,np.cos(np.pi/4),np.cos(np.pi/4),0],[0,0,0,1]])

	# print(M3)
	R4=M1.dot(M2).dot(M3)
	B8=np.sqrt(2)/2*np.array([[1,0,0,0,0,0,0,1],[0,1,0,0,0,0,1,0],[0,0,1,0,0,1,0,0],[0,0,0,1,1,0,0,0],[0,0,0,1,-1,0,0,0],[0,0,1,0,0,-1,0,0],[0,1,0,0,0,0,-1,0],[1,0,0,0,0,0,0,-1]])

	#计算C4
	P2=np.array([[1,0],[0,1]])
	C2=np.sqrt(2)/2*np.array([[1,1],[1,-1]])
	R2=np.array([[np.sin(np.pi/8),np.cos(np.pi/8)],[-np.sin(3*np.pi/8),np.cos(3*np.pi/8)]])
	B4=np.sqrt(2)/2*np.array([[1,0,0,1],[0,1,1,0],[0,1,-1,0],[1,0,0,-1]])
	
	tmp1=np.concatenate((P2.T.dot(C2),np.zeros((2,2))),axis=1)
	tmp2=np.concatenate((np.zeros((2,2)),R2),axis=1)
	tmp=np.concatenate((tmp1,tmp2))
	C4=P4.dot(tmp).dot(B4)

	tmp1=np.concatenate((P4.T.dot(C4),np.zeros((4,4))),axis=1)
	tmp2=np.concatenate((np.zeros((4,4)),R4),axis=1)
	tmp=np.concatenate((tmp1,tmp2))
	#得到正交矩阵C8
	C8=P8.dot(tmp).dot(B8)

	Y=C8.T.dot(X).dot(C8)
	return Y


def main():
	#测试，给定一个8*8的矩阵
	X=np.array([[42,66,68,66,42,66,68,66],[92,4,76,17,42,66,68,66],[79,85,74,71,42,66,68,66],[96,93,39,3,42,66,68,66],[42,66,68,66,42,66,68,66],[92,4,76,17,42,66,68,66],[79,85,74,71,42,66,68,66],[96,93,39,3,42,66,68,66]],dtype=np.float64)
	#计算DCT变换后的结果
	Y=cal_dct(X)
	print(Y)
	print('*************************')
	#调用cv2库中的dct函数
	YY=cv2.dct(X)
	print(YY)
	print('************************')
	#提出的快速计算算法
	YYY=fast_DCT(X)
	print(YYY)
	#IDCT
	X=cal_idct(Y)
	print(X)
	print('*************************')
	#调用cv2库中的dct函数
	X=cv2.idct(YY)
	print(X)
	print('************************')
	#提出的快速计算算法
	X=fast_IDCT(YYY)
	print(X)

if __name__ =='__main__':
	main()

实验结果

[[ 4.84750000e+02  6.41525518e+00  8.08716048e+01  1.94719777e+01
  -3.57500000e+01  1.34448255e+01  3.38807990e+01  9.57461504e+00]
 [-4.32489152e+00 -1.36497986e+01 -2.33629144e+01 -1.64769788e+01
   2.82560597e+00  1.36169047e+01  8.42538557e+00  5.23162272e-01]
 [-3.74467051e-14  2.45392075e-15 -1.42634402e-14 -7.86703667e-15
   6.20186521e-15  3.93768720e-15 -8.79218256e-15 -1.05081179e-14]
 [-1.39699475e+01 -2.88766884e+01 -3.89941365e+01 -2.50078137e+01
   6.99429145e+00  2.71861709e+01  2.25130198e+01  8.55081980e+00]
 [-6.25000000e+00 -6.21998536e-01  1.07158195e+01  4.11351653e+00
  -1.97500000e+01 -3.93065081e+01 -3.88045901e+01 -2.20780551e+01]
 [ 2.44075900e+01  2.20631412e+01  7.45093787e-02 -8.95596469e+00
  -8.24036938e+00 -1.61533515e+01 -3.05597165e+01 -2.77419121e+01]
 [-3.61631546e-13 -2.68558280e-14 -1.14408027e-13 -5.13334521e-14
   4.91046416e-14  4.25111508e-14  2.27897333e-14  5.10099810e-15]
 [ 3.17100998e+01  8.38102665e+00 -4.85264557e+01 -4.92516810e+01
  -7.86238834e+00  1.40906021e+00 -3.34341090e+01 -4.51890361e+01]]
*************************
[[ 4.84750000e+02  6.41525518e+00  8.08716048e+01  1.94719777e+01
  -3.57500000e+01  1.34448255e+01  3.38807990e+01  9.57461504e+00]
 [-4.32489152e+00 -1.36497986e+01 -2.33629144e+01 -1.64769788e+01
   2.82560597e+00  1.36169047e+01  8.42538557e+00  5.23162272e-01]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]
 [-1.39699475e+01 -2.88766884e+01 -3.89941365e+01 -2.50078137e+01
   6.99429145e+00  2.71861709e+01  2.25130198e+01  8.55081980e+00]
 [-6.25000000e+00 -6.21998536e-01  1.07158195e+01  4.11351653e+00
  -1.97500000e+01 -3.93065081e+01 -3.88045901e+01 -2.20780551e+01]
 [ 2.44075900e+01  2.20631412e+01  7.45093787e-02 -8.95596469e+00
  -8.24036938e+00 -1.61533515e+01 -3.05597165e+01 -2.77419121e+01]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
   0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 3.17100998e+01  8.38102665e+00 -4.85264557e+01 -4.92516810e+01
  -7.86238834e+00  1.40906021e+00 -3.34341090e+01 -4.51890361e+01]]
************************
[[ 4.84750000e+02  6.41525518e+00  8.08716048e+01  1.94719777e+01
  -3.57500000e+01  1.34448255e+01  3.38807990e+01  9.57461504e+00]
 [-4.32489152e+00 -1.36497986e+01 -2.33629144e+01 -1.64769788e+01
   2.82560597e+00  1.36169047e+01  8.42538557e+00  5.23162272e-01]
 [-6.61401448e-15  3.17122151e-16 -3.92887044e-15  2.51146777e-15
   8.22335925e-15  3.09337706e-15  3.90041359e-15  2.22506330e-15]
 [-1.39699475e+01 -2.88766884e+01 -3.89941365e+01 -2.50078137e+01
   6.99429145e+00  2.71861709e+01  2.25130198e+01  8.55081980e+00]
 [-6.25000000e+00 -6.21998536e-01  1.07158195e+01  4.11351653e+00
  -1.97500000e+01 -3.93065081e+01 -3.88045901e+01 -2.20780551e+01]
 [ 2.44075900e+01  2.20631412e+01  7.45093787e-02 -8.95596469e+00
  -8.24036938e+00 -1.61533515e+01 -3.05597165e+01 -2.77419121e+01]
 [ 4.71027738e-16  1.52198036e-15 -1.05487038e-15 -1.60330487e-15
  -1.57009246e-16 -4.26576458e-16  1.96645187e-15  3.96516478e-16]
 [ 3.17100998e+01  8.38102665e+00 -4.85264557e+01 -4.92516810e+01
  -7.86238834e+00  1.40906021e+00 -3.34341090e+01 -4.51890361e+01]]
[[42. 66. 68. 66. 42. 66. 68. 66.]
 [92.  4. 76. 17. 42. 66. 68. 66.]
 [79. 85. 74. 71. 42. 66. 68. 66.]
 [96. 93. 39.  3. 42. 66. 68. 66.]
 [42. 66. 68. 66. 42. 66. 68. 66.]
 [92.  4. 76. 17. 42. 66. 68. 66.]
 [79. 85. 74. 71. 42. 66. 68. 66.]
 [96. 93. 39.  3. 42. 66. 68. 66.]]
*************************
[[42. 66. 68. 66. 42. 66. 68. 66.]
 [92.  4. 76. 17. 42. 66. 68. 66.]
 [79. 85. 74. 71. 42. 66. 68. 66.]
 [96. 93. 39.  3. 42. 66. 68. 66.]
 [42. 66. 68. 66. 42. 66. 68. 66.]
 [92.  4. 76. 17. 42. 66. 68. 66.]
 [79. 85. 74. 71. 42. 66. 68. 66.]
 [96. 93. 39.  3. 42. 66. 68. 66.]]
************************
[[42. 66. 68. 66. 42. 66. 68. 66.]
 [92.  4. 76. 17. 42. 66. 68. 66.]
 [79. 85. 74. 71. 42. 66. 68. 66.]
 [96. 93. 39.  3. 42. 66. 68. 66.]
 [42. 66. 68. 66. 42. 66. 68. 66.]
 [92.  4. 76. 17. 42. 66. 68. 66.]
 [79. 85. 74. 71. 42. 66. 68. 66.]
 [96. 93. 39.  3. 42. 66. 68. 66.]]
[Finished in 0.4s]