转载 sohu blog:编程铺路石,作者:shepherd
http://zqw100.blog.sohu.com/6110148.html
jpeg解码的核心问题是IDCT(逆离散余弦变换),这也是整个JPEG解码过程中耗时最长的步骤。本文将讨论一个8×8矩阵的IDCT简化的问题。
首先看一下IDCT的公式:
7 7 2*x+1 2*y+1
f(x,y) = sum sum alpha(u)*alpha(v)*F(u,v)*cos (------- *u*PI)* cos (------ *v*PI)
u=0 v=0 16 16
其中x,y=0,1...7
{ 1/sqrt(8) (u==0)
alpha(u) = {
{ 1/2 (u!=0)
上式中F(u,v)代表原始数据(即DCT数据)矩阵中的元素(u,v),
f(x,y)指解出的像素矩阵中的元素(x,y)
那个公式看晕了吧?好等你明白过来我继续说
……
明白了哈,我接着说。
看到了那个公式,相信很多人的第一反映就是直接展开,用嵌套的FOR循环进行运算。对于这种方法,怎么说呢,也许用主频为100G的CPU,解码的速度还可以忍受。不相信,你可以自己写一下,看那是几层循环嵌套,再计算一下解码这64个点一共需要运算多少次。现在理解了吧,直接展开的方法是不行滴。
那么如何简化呢?首先就是要把这个二维运算,分解成两个一维运算。换句话说,先对矩阵的每行进行一次一维IDCT,把结果保存在一个中间数组里,然后对这个数组的每列进行一次一维IDCT,就得到最终的结果了。如果你数学基础还可以的话,可以直接变换上面的公式;如果你的数学比较烂,可以用具体的数字进行演算,你会发现在上文那个丑陋的多重嵌套FOR循环中,有许多运算是重复进行的。
好休息一会,自己演算一下把
……
现在我们来讲解一维IDCT的简化。首先看公式
1 7 (2*x+1)*i*PI
p(x) = --- * ∑ { c(i) * DCT(i) * cos ------------ }
2 i=0 16
{ 1/sqrt(2) (i==0)
其中 c(i)={
{ 1 (i!=0)
现在假设f(t) = cos(PI*t/16)
展开上式可以得到:
P(0) = (1/4) * ( f(0)*1/sqrt(2)*DCT(0) + f(1)*1*DCT(1) + …… + f(7)*1*DCT(7) )
……
P(7) = (1/4) * ( f(0)*1/sqrt(2)*DCT(0) + f(15)*1*DCT(1) + …… + f(105)*1*DCT(7) )
暂时总结一下,将t值整理成一个矩阵:
0 1 2 3 4 5 6 7
0 3 6 9 12 15 18 21
0 5 10 15 20 25 30 35
0 7 14 21 28 35 42 49
0 9 18 27 36 45 54 63
0 11 22 33 44 55 66 77
0 13 26 39 52 65 78 91
0 15 30 45 60 75 90 105
很有规律吧?所以要接着简化啊。
1,f(t) = cos(t*PI/16),说明f(t)是余弦函数,可以利用余弦函数的周期性减少未知量的个数。
2,由于f(0)*1/sqrt(2) = 1*sqrt(2)/2 = cos(PI/4) = cos(PI*4/16) = f(4),利用这一特点可以进一步简化
假设Cn = cos(n*PI/16) = f(n),则一维IDCT的结果可以整理如下:
DCT(0)* DCT(1)* DCT(2)* DCT(3)* DCT(4)* DCT(5)* DCT(6)* DCT(7)*
P(0) = C4 C1 C2 C3 C4 C5 C6 C7
P(1) = C4 C3 C6 -C7 -C4 -C1 -C2 -C5
P(2) = C4 C5 -C6 -C1 -C4 C7 C2 C3
P(3) = C4 C7 -C2 -C5 C4 C3 -C6 -C1
P(4) = C4 -C7 -C2 C5 C4 -C3 -C6 C1
P(5) = C4 -C5 -C6 C1 -C4 -C7 C6 -C3
P(6) = C4 -C3 C6 C7 -C4 C1 -C2 C5
P(7) = C4 -C1 C2 -C3 C4 -C5 C6 -C7
这里C1-C7都已经是常数了,可以事先求出。
看花眼了吧?好休息一下接着说
……
仔细观察上面的结果,可以找到规律:
偶数列的上四行和下四行是对称的;奇数行的上四行和下四行是反对称的。
假设原始数组是x[],转化后的数组是y[],那么
a0 = x[0]*C4 + x[2]*C2 + x[4]*C4 + x[6]*C6
a1 = x[0]*C4 + x[2]*C6 - x[4]*C4 - x[6]*C2
a2 = x[0]*C4 - x[2]*C6 - x[4]*C4 + x[6]*C2
a3 = x[0]*C4 - x[2]*C2 + x[4]*C4 - x[6]*C6
b0 = x[1]*C1 + x[3]*C3 + x[5]*C5 + x[7]*C7
b1 = x[1]*C3 - x[3]*C7 - x[5]*C1 - x[7]*C5
b2 = x[1]*C5 - x[3]*C1 + x[5]*C7 + x[7]*C3
b3 = x[1]*C7 - x[3]*C5 + x[5]*C3 - x[7]*C1
y[0] = a0 + b0
y[7] = a0 - b0
y[1] = a1 + b1
y[6] = a1 - b1
y[2] = a2 + b2
y[5] = a2 - b2
y[3] = a3 + b3
y[4] = a3 - b3
以上是别人文章给出的代码,实际上还可以进一步化简a部分
t0 = x[0] * C4
t1 = x[2] * C2
t2 = x[2] * C6
t3 = x[4] * C4
t4 = x[6] * C6
t5 = x[6] * C2
a0 = t0 + t1 + t3 + t4
a1 = t0 + t2 - t3 - t5
a2 = t0 - t2 - t3 + t5
a3 = t0 - t1 + t3 - t4
后面的B部分和Y部分不用再变了。
至此化简工作结束。这样一个8*8的矩阵,只需进行两次类似的四则运算就可以了。
我统计了一下,每次一维IDCT是22次乘法,32次加减法,而且经过处理后都是整数运算。
相信你都晕乎乎了吧,好在都说完了,嘿嘿。