本文提出了一种新的快速离散余弦变换（FDCT）算法，通过矩阵分解和实数运算，降低了计算复杂度，比传统基于FFT的方法快约6倍。算法使用了交替的余弦/正弦蝴蝶矩阵和二进制矩阵，适用于图像编码和其他应用中的DCT计算。实验部分展示了算法的有效性和与标准库实现的一致性。
摘要
本文提出了一种快速离散余弦变换算法，与传统的使用快速傅里叶变换的离散余弦变换算法相比，计算复杂度改善了6倍，该算法以矩阵的形式导出，并用信号流图来说明，信号流图可以很容易地转换为硬件或软件实现。
引言
离散余弦变换（DCT）已经被成功应用到高分辨率图像编码中。传统的实现DCT变换的方法是利用双倍大小的快速傅里叶（FFT）算法，在计算过程中使用了复数计算。由于缺乏有效的算法，DCT在各种各样的应用中并没有像它的特性所暗示的那样广泛。本文描述了一种只涉及实数操作的更有效的计算一组N个点的快速离散余弦变换(FDCT)算法。该算法可以扩展到任意期望的值 N = 2 m , m > = 2 N=2^m,m>=2 N=2m,m>=2。这种推广包括交替的余弦/正弦蝴蝶矩阵和二进制矩阵，以重新排列矩阵元素，使其在每个其他节点上保持可识别的位反转模式。泛化并不是唯一的——已经发现了几种替代方法，但这里描述的方法似乎是最简单的解释。它不一定是最有效的FDCT，但代表了一种方法扩展的技术。这种方法采用了 ( 3 N / 2 ) ( l o g 2 N − 1 ) + 2 (3N/2)(log_2N-1)+2 (3N/2)(log2N−1)+2个实数加法和 N l o g 2 N − 3 N / 2 + 4 Nlog_2N-3N/2+4 Nlog2N−3N/2+4个实数乘法，这大约是使用双倍大小FFT的传统方法的6倍快。
离散余弦变换
离散函数
f
(
j
)
，
j
=
0
,
1
,
.
.
.
,
N
−
1
f(j)，j=0,1,...,N-1
f(j)，j=0,1,...,N−1的离散余弦变换定义如下：
F
(
k
)
=
2
c
(
k
)
N
∑
j
=
0
N
−
1
f
(
j
)
c
o
s
[
(
2
j
+
1
)
k
π
2
N
]
,
k
=
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(k)=N2c(k)j=0∑N−1f(j)cos[2N(2j+1)kπ],k=0,1,...,N−1
逆变换定义如下：
f
(
j
)
=
∑
k
=
0
N
−
1
c
(
k
)
F
(
k
)
c
o
s
[
(
2
j
+
1
)
k
π
2
N
]
,
j
=
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
f(j)=k=0∑N−1c(k)F(k)cos[2N(2j+1)kπ],j=0,1,...,N−1
其中，
c
(
k
)
=
1
2
,
k
=
0
c
(
k
)
=
1
,
k
=
1
,
2
,
.
.
.
,
N
−
1
c(k)=\frac{1}{\sqrt{2}}, k=0\\ c(k)=1,k=1,2,...,N-1
c(k)=21,k=0c(k)=1,k=1,2,...,N−1
该变换比任何已知的快速计算算法拥有更高的能量压缩特性。该变换还具有循环卷积-乘法关系，可方便地应用于线性系统理论中。
快速计算算法
一个N*1的数据向量[f]的离散余弦变换可以用矩阵形式表达如下：
[
F
]
=
2
N
[
A
N
]
[
f
]
[F]=\frac{2}{N}[A_N][f]
[F]=N2[AN][f]
其中，
[
A
N
]
=
[
c
(
k
)
c
o
s
(
2
j
+
1
)
k
π
/
2
N
]
;
j
,
k
=
0
,
1
,
.
.
.
,
(
N
−
1
)
[A_N]=[c(k)cos(2j+1)k\pi/2N];j,k=0,1,...,(N-1)
[AN]=[c(k)cos(2j+1)kπ/2N];j,k=0,1,...,(N−1),[F]是N*1的转换后的向量。本文提出的快速计算算法是基于
[
A
N
]
[A_N]
[AN]矩阵的矩阵分解。如下所示，该矩阵可以写成下面的递归形式：
[
A
N
]
=
[
P
N
]
[
A
N
/
2
0
0
R
N
/
2
]
[
B
N
]
[A_N]=[P_N]\left[\begin{array}{cc}A_{N/2}&0\\0&R_{N/2}\end{array}\right][B_N]
[AN]=[PN][AN/200RN/2][BN]
其中，
[
B
N
]
[B_N]
[BN]定义如下：
[
R
N
/
2
]
=
[
c
(
k
)
c
o
s
(
2
j
+
1
)
(
2
k
+
1
)
π
2
N
]
;
j
,
k
=
0
,
1
,
.
.
.
,
N
2
−
1
[R_{N/2}]=[c(k)cos\frac{(2j+1)(2k+1)\pi}{2N}];j,k=0,1,...,\frac{N}{2}-1
[RN/2]=[c(k)cos2N(2j+1)(2k+1)π];j,k=0,1,...,2N−1
[
P
N
]
[P_N]
[PN]是一个N x N的排列矩阵，它将变换后的向量从位逆序排列为一个自然的顺序。2*2DCT变换可以写成如下形式
[
A
2
]
=
1
(
2
)
[
1
1
1
−
1
]
[A_2]=\frac{1}{\sqrt(2)}\left[\begin{array}{cc}1&1\\1&-1\end{array}\right]
[A2]=(2)1[111−1]
从
[
A
N
]
[A_N]
[AN]的递归性质可以看出，只要找到一种通用的
[
R
N
/
2
]
[R_{N/2}]
[RN/2]矩阵分解方法，
A
2
A_2
A2就可以扩展到更高阶的矩阵。
下面的讨论提出了一种分解 [ R N / 2 ] [R_{N/2}] [RN/2] 矩阵的系统方法,需要强调的是这种分解方法不是唯一的，也不是最优的。另外有几种方法需要较少的计算步骤，但它们对更大的尺寸没有明显的泛化。
[
R
N
/
2
]
[R_{N/2}]
[RN/2]矩阵可以用下面的方式分解成（
2
l
o
g
2
N
−
3
2log_2N-3
2log2N−3）个矩阵：
[
R
N
/
2
]
=
[
M
1
]
[
M
2
]
[
M
3
]
[
M
4
]
.
.
.
[
M
(
2
l
o
g
2
N
−
3
)
]
[R_{N/2}]=[M1][M2][M3][M4]...[M(2log_2N-3)]
[RN/2]=[M1][M2][M3][M4]...[M(2log2N−3)]
矩阵有四种不同的类型
类型1：[M1],第一个矩阵
类型2： [ M ( 2 l o g 2 N − 3 ) ] [M(2log_2N-3)] [M(2log2N−3)]，最后一个矩阵
类型3：[Mq],奇数矩阵，如[M3],[M5]等
类型4：[Mp],偶数矩阵，如[M2],[M4]等
在详细描述这四种矩阵之前，给出了以下定义:
B
N
=
[
I
N
/
2
I
‾
N
/
2
I
‾
N
/
2
−
I
N
/
2
]
B_N=\left[\begin{array}{cc}I_{N/2}&\overline{I}_{N/2}\\\overline{I}_{N/2} &-I_{N/2}\end{array}\right]
BN=[IN/2IN/2IN/2−IN/2]
B N ∗ = [ − I N / 2 I ‾ N / 2 I ‾ N / 2 I N / 2 ] B_N^*=\left[\begin{array}{cc}-I_{N/2}&\overline{I}_{N/2}\\\overline{I}_{N/2} &I_{N/2}\end{array}\right] BN∗=[−IN/2IN/2IN/2IN/2]
I N / 2 I_{N/2} IN/2是阶数为 N 2 \frac{N}{2} 2N的单位矩阵
I ‾ N / 2 \overline{I}_{N/2} IN/2是反对角单位矩阵
[ S i k ] = s i n k π i [ I N / 2 i ] [S_i^k]=sin\frac{k\pi}{i}[I_{N/2i}] [Sik]=sinikπ[IN/2i]
[ s ‾ i k ] = s i n k π i [ I ‾ N / 2 i ] [\overline{s}_i^k]=sin\frac{k\pi}{i}[\overline{I}_{N/2i}] [sik]=sinikπ[IN/2i]
[ C i k ] = c o s k π i [ I N / 2 i ] [C_i^k]=cos\frac{k\pi}{i}[I_{N/2i}] [Cik]=cosikπ[IN/2i]
[ C ‾ i k ] = c o s k π i [ I ‾ N / 2 i ] [\overline{C}_i^k]=cos\frac{k\pi}{i}[\overline{I}_{N/2i}] [Cik]=cosikπ[IN/2i]
上式中，单位矩阵 [ I N / 2 i ] [I_{N/2i}] [IN/2i]表示 对角sin或者cos矩阵 [ S i k ] , [ S ‾ i k ] , [ C i k ] , [ C ‾ i k ] [S_i^k],[\overline{S}_i^k],[C_i^k],[\overline{C}_i^k] [Sik],[Sik],[Cik],[Cik]的阶数，另外i>N/2 时， [ I N / 2 i ] ≡ 1 [I_{N/2i}]\equiv1 [IN/2i]≡1。
实验
实验代码
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]
参考
【A-Fast-Computational-Algorithm-for-the-Discrete-Cosine-Transform】
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
