今天的实验里面有个涉及到DCT的,虽然不用我们具体去写DCT的代码,不过还是来看一下DCT到底是个什么玩意吧。
DCT基本算法原理
DCT: 离散余弦变换(Discrete Cosine Transform),
一维DCT变换
基本定义
如果
{
X
(
m
)
∣
m
=
0
,
1
,
.
.
.
,
N
−
1
}
\{X(m)|m=0,1,...,N-1\}
{X(m)∣m=0,1,...,N−1}是有限带宽信号x(t)的采样序列,则一维离散余弦变换为:
Y
(
u
)
=
C
(
u
)
2
N
∑
m
=
0
N
−
1
X
(
m
)
cos
(
2
m
+
1
)
μ
π
2
N
,
u
=
1
,
2
,
.
.
.
,
N
−
1
Y(u) = C(u) \sqrt{\frac{2}{N}}\sum\limits_{m=0}^{N-1}X(m)\cos \frac{(2m+1)\mu\pi}{2N},u=1,2,...,N-1
Y(u)=C(u)N2m=0∑N−1X(m)cos2N(2m+1)μπ,u=1,2,...,N−1
其中
C ( u ) = { 1 / 2 , u = 0 1 , e l s e C(u) =\begin{cases} 1/\sqrt{2},&u=0\\1,& else\end{cases} C(u)={1/2,1,u=0else
其逆变换为
X ( m ) = 2 N ∑ m = 0 N − 1 C ( u ) Y ( u ) cos ( 2 m + 1 ) u π 2 N , m = 1 , 2 , . . . , N − 1 X(m) =\sqrt{\frac{2}{N}}\sum\limits_{m=0}^{N-1}C(u)Y(u)\cos \frac{(2m+1)u\pi}{2N},m=1,2,...,N-1 X(m)=N2m=0∑N−1C(u)Y(u)cos2N(2m+1)uπ,m=1,2,...,N−1
二者的变换核为
a ( u , m ) = C ( u ) 2 N cos ( 2 m + 1 ) u π 2 N a(u,m) = C(u)\sqrt{\frac{2}{N}}\cos \frac{(2m+1)u\pi}{2N} a(u,m)=C(u)N2cos2N(2m+1)uπ
二维DCT变换
对二维图像信号数据而言,有推广
Y ( u , v ) = C ( u ) C ( v ) 2 M N ∑ m = 0 M − 1 ∑ n = 0 N − 1 X ( m , n ) cos ( 2 m + 1 ) u π 2 M cos ( 2 n + 1 ) v π 2 N , u = 1 , 2 , . . . , M − 1 , v = 0 , 1 , . . . , N − 1 Y(u,v) = C(u) C(v){\frac{2}{\sqrt{MN}}}\sum\limits_{m=0}^{M-1} \sum\limits_{n=0}^{N-1} X(m,n)\cos \frac{(2m+1)u\pi}{2M} \cos \frac{(2n+1)v\pi}{2N} ,u=1,2,...,M-1,v = 0,1,...,N-1 Y(u,v)=C(u)C(v)MN2m=0∑M−1n=0∑N−1X(m,n)cos2M(2m+1)uπcos2N(2n+1)vπ,u=1,2,...,M−1,v=0,1,...,N−1
C ( u ) , C ( v ) = { 1 / 2 , u , v = 0 1 , e l s e C(u),C(v) =\begin{cases} 1/\sqrt{2},&u,v=0\\1,& else\end{cases} C(u),C(v)={1/2,1,u,v=0else
逆变换为
X ( m , n ) = 2 M N ∑ m = 0 M − 1 ∑ n = 0 N − 1 C ( u ) C ( v ) Y ( u , v ) cos ( 2 m + 1 ) u π 2 M cos ( 2 n + 1 ) v π 2 N , m = 1 , 2 , . . . , M − 1 , n = 0 , 1 , . . . , N − 1 X(m,n) = \frac{2}{\sqrt{MN}}\sum\limits_{m=0}^{M-1} \sum\limits_{n=0}^{N-1} C(u) C(v) Y(u,v)\cos \frac{(2m+1)u\pi}{2M} \cos \frac{(2n+1)v\pi}{2N} ,m=1,2,...,M-1,n=0,1,...,N-1 X(m,n)=MN2m=0∑M−1n=0∑N−1C(u)C(v)Y(u,v)cos2M(2m+1)uπcos2N(2n+1)vπ,m=1,2,...,M−1,n=0,1,...,N−1
变换核是两次一维DCT变换核的乘积
a ( u , v , m , n ) = a 1 ( u , m ) a 2 ( v , n ) = C ( u ) 2 M cos ( 2 m + 1 ) u π 2 M C ( v ) 2 N cos ( 2 n + 1 ) v π 2 N a(u,v,m,n) = a_1(u,m)a_2(v,n)=C(u)\sqrt{\frac{2}{M}}\cos \frac{(2m+1)u\pi}{2M}C(v)\sqrt{\frac{2}{N}}\cos \frac{(2n+1)v\pi}{2N} a(u,v,m,n)=a1(u,m)a2(v,n)=C(u)M2cos2M(2m+1)uπC(v)N2cos2N(2n+1)vπ
DCT变换在压缩中的应用
DCT的变换核构成的基向量与图像内容无关,而且变换核是可以分离的,既二维DCT可以用两次一维DCT来完成,使得数学运算难度大大简化,再配以已经发现的其它快速算法,使得DCT编码得到了广泛的应用。
特点
- DCT应用于图像数据压缩,可以减少代表图像亮度(或色度)层次数码信息,达到数据压缩的目的。
- 在编码变换过程发现图像细节的位置,以便删去或略去对视觉不敏感的部分,而更加突出视觉的敏感部分,通过选择主要数据来传输、重视图像。
信号的频谱
由图像内取出一个区块,分成8×8个像素的64格阵列,经过对逐个像素的亮度(或讨论色度)数值取样,并将像素的亮度数值列成矩阵形表格。然后利用离散余弦变换(DCT)可将各空间取样值转变为频率域的数值,即为DCT系数。
64点阵列,可得到64个DCT系数,其可以等效为一个直流平均值和63个不同频率余弦波幅值组成的64个点阵列。
这64个点代表频谱上的64个频率,相当于是64个正交基函数的系数。左上角是直流和低频部分,右下角是高频部分。(从表达式也可看出u,v越大的地方基函数频率越高)
可以按照Z字扫描顺序,将DCT系数从二维数据转变为一维数据。
该一维序列的首项是平均亮度值,后面的项代表亮度的起伏程度。
人眼睛对图像的亮度信息有较高的相对灵敏度,对图像的彩色信息不够敏感;还有,人眼睛对图像信息的低频分量具有较高的视觉灵敏度。经Z字形字扫描后所形成的数据系列,恰好与人眼睛对图像信息的敏感程度形成良好的对应关系。根据视觉生理的上述规律,可对图像数据进行压缩。
系数的量化
对每种频率分量设定不同的折算值,将转换得到的DCT系数再次进行折算,可以进一步突出视觉效果影响大的成分,而消弱或忽略视觉效果影响小的成分。这种处理方法称为量化处理,简称Q处理。
如: 将区块的64个系数分别除以量化表中对应位置量化步长,再进行四舍五入取整
对于64点阵列的64个系数来说,对应了64种不同频率,可使用64个不同的折算值。通常称这64个折算值为量化表,每个折算值称为量化步长,或称量化值。在64点阵列中,左上角的数据量化值较小,右下角的数据量化值较大。
经过量化处理后的DCT系数矩阵,可出现许多零值。一般左上角位置的数据的商数是非0,在右下角位置的数据的商数很小,经四舍五入取整值后可简写为0。
果冻君的提问:所以可以压缩成稀疏矩阵吗??
实际操作归纳
- 要将输入的图像分解为 8*8 的块
- 后对每个块进行二维离散变换把每个块转变成 64 个 DCT 系数值
- 后将变换得到的 DCT 系数进行编码和传送
- 解码时对每个块进行二维 DCT 反变换
- 将反变换后的块组合成图像
将8*8的大小代入原公式,有
Y ( u , v ) = 1 4 C ( u ) C ( v ) ∑ m = 0 7 ∑ n = 0 7 X ( m , n ) cos ( 2 m + 1 ) u π 16 cos ( 2 n + 1 ) v π 16 , u = 1 , 2 , . . . , M − 1 , v = 0 , 1 , . . . , N − 1 Y(u,v) = \frac{1}{4}C(u) C(v)\sum\limits_{m=0}^{7} \sum\limits_{n=0}^{7} X(m,n)\cos \frac{(2m+1)u\pi}{16} \cos \frac{(2n+1)v\pi}{16} ,u=1,2,...,M-1,v = 0,1,...,N-1 Y(u,v)=41C(u)C(v)m=0∑7n=0∑7X(m,n)cos16(2m+1)uπcos16(2n+1)vπ,u=1,2,...,M−1,v=0,1,...,N−1
可以等效为两个N=8的一维DCT变换。
对课程案例的分析
1维DCT变换子函数
其源代码为
void dct_1d(dct_data_t src[DCT_SIZE], dct_data_t dst[DCT_SIZE])
{
unsigned int k, n;
int tmp;
const dct_data_t dct_coeff_table[DCT_SIZE][DCT_SIZE] = {
#include "dct_coeff_table.txt"
};
DCT_Outer_Loop:
for (k = 0; k < DCT_SIZE; k++) {
DCT_Inner_Loop:
for(n = 0, tmp = 0; n < DCT_SIZE; n++) {
int coeff = (int)dct_coeff_table[k][n];
tmp += src[n] * coeff;
}
dst[k] = DESCALE(tmp, CONST_BITS);
}
}
引用的系数表为
8192, 8192, 8192, 8192, 8192, 8192, 8192, 8192,
11363, 9633, 6436, 2260, -2260, -6436, -9632,-11362,
10703, 4433, -4433,-10703,-10703, -4433, 4433, 10703,
9633, -2260,-11362, -6436, 6436, 11363, 2260, -9632,
8192, -8192, -8192, 8192, 8192, -8191, -8191, 8192,
6436,-11362, 2260, 9633, -9632, -2260, 11363, -6436,
4433,-10703, 10703, -4433, -4433, 10703,-10703, 4433,
2260, -6436, 9633,-11362, 11363, -9632, 6436, -2260
由代码可知,输入数据和输出数据都是一维数组,其首先将系数表文件导入数组中,然后对数据进行遍历。
Outer循环对k,也就是系数表的行进行遍历
Inner循环对数据和系数表的列n进行遍历
因此其每次Inner循环都是将数据src中的每个值与系数表第k行的每个数值相乘再相加,就得到了输出数据的第k个值。
从一维DCT变换的公式可知,内层循环做累加时所乘系数应该为
a ( u , m ) = C ( u ) 2 N cos ( 2 m + 1 ) μ π 2 N a(u,m)=C(u) \sqrt{\frac{2}{N}}\cos \frac{(2m+1)\mu\pi}{2N} a(u,m)=C(u)N2cos2N(2m+1)μπ
在次代码中,u为j,m为n。
最后的Descale的宏定义为
#define DESCALE(x,n) (((x) + (1 << ((n)-1))) >> n)
其效果相当于给x加一个2(n-1)的常量便宜再除以2n
很多本来就不够大的数字会被移位变成0。
二维DCT子函数
其源码为
void dct_2d(dct_data_t in_block[DCT_SIZE][DCT_SIZE],
dct_data_t out_block[DCT_SIZE][DCT_SIZE])
{
dct_data_t row_outbuf[DCT_SIZE][DCT_SIZE];
dct_data_t col_outbuf[DCT_SIZE][DCT_SIZE], col_inbuf[DCT_SIZE][DCT_SIZE];
unsigned i, j;
// DCT rows
Row_DCT_Loop:
for(i = 0; i < DCT_SIZE; i++) {
dct_1d(in_block[i], row_outbuf[i]);
}
// Transpose data in order to re-use 1D DCT code
Xpose_Row_Outer_Loop:
for (j = 0; j < DCT_SIZE; j++)
Xpose_Row_Inner_Loop:
for(i = 0; i < DCT_SIZE; i++)
col_inbuf[j][i] = row_outbuf[i][j];
// DCT columns
Col_DCT_Loop:
for (i = 0; i < DCT_SIZE; i++) {
dct_1d(col_inbuf[i], col_outbuf[i]);
}
// Transpose data back into natural order
Xpose_Col_Outer_Loop:
for (j = 0; j < DCT_SIZE; j++)
Xpose_Col_Inner_Loop:
for(i = 0; i < DCT_SIZE; i++)
out_block[j][i] = col_outbuf[i][j];
}
输入数据为二维数组,输出数据也是二维数组。
其首先进行Row_DCT_Loop,对ROWs进行DCT变换
接下来的Xpose循环实际上是对数组进行了转置,这样新的数组相当于是进行了Cols的DCT变换
然后再进行一次DCT循环。
最后再转置变回原来的方向。
参考文献
[1]孙娟. DCT原理及其在视频压缩编码中的实现[EB/OL]. 北京:中国科技论文在线 [2008-07-23].
[2] https://www.cnblogs.com/stnlcd/p/7261842.html