探讨JPEG：IDCT

 

转载 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次加减法，而且经过处理后都是整数运算。
相信你都晕乎乎了吧，好在都说完了，嘿嘿。