AAN算法描述及其实现

AAN算法描述及其实现

 

最近在学习DCT和MDCT发现网络上JPEG编码的资料蛮多, 而且DCT是JEPG压缩里头比较重要的一个环节. 单单描述DCT算法的资料比较少, 看看JEPG编码里头的DCT部分也不失为一个好的学习方法.

 

 

二维DCT

 

如果按照二维DCT公式来计算8*8像素的话, 获得一个DCT系数需要8*8=64次乘法和8*8=64次加法, 而完成整个8*8像素的DCT需要8*8*8*8=4096次乘法和8*8*8*8=4096次加法. 计算量是相当的大!

 

伪代码(N=8, PI=3.1415926):

 

int           u, v, i, j;   // u, v描述输出; i,j描述输入.

double        c_u, c_v;  //

unsigned char in_data[N][N] = { 数据 }; // 注意输入数据的类型, 色彩空间的范围是[0, 255]

double        out_data[N][N] = null;

 

for(u=0; u<N; u++)

{

    for(v=0, out_data[u][v]=0.0; v<N; v++)

{

    // 内层处理输入

       for(i=0; i<N; i++)

       {

           for(j=0; j<N; j++)

           {

              if(u==0)   c_u = sqrt(1.0/N);

              else       c_u = sqrt(2.0/N);

              if(v==0)   c_v = sqrt(1.0/N);

              else       c_v = sqrt(2.0/N);

             

    out_data[u][v] += c_u*c_v*in_data[i][j]

*cos((2*i+1)*u*PI/(2*N))*cos((2*j+1)*v*PI/(2*N)) / 4;

           }

}

}

}

这就是和二维DCT公式对应的伪代码, 未做任何优化, 运算量非常庞大, 复杂度为(N*N*N*N). 无论是软件还是硬件实现(即使是对这个算法进行优化, 如果架构不修改的化)都是非常的困难的!

 

 

一维DCT

 

二维DCT可以转换为一维DCT来计算, 那么8行8列的话, 对于一行来说需要计算的是(8*8)次乘法和(8*8)次加法, 而8行就需要计算8*(8*8)次乘法和8*(8*8)次加法. 对于列来说也是同样. 那么处理完一个block就需要计算2*(8*(8*8))=1024次乘法和2*(8*(8*8))=1024次加法, 运算量变为二维计算的1/4.

但是仍然显得够大.

 

伪代码

So tired! 以下省略代码100行...

 

 

 

 

 

AAN算法描述

a1=0.707, a2=0.541, a3=0.707, a4=1.307, a1=0.383

整个AAN算法流程如上图. 其中黑点表示加法, 箭头表示乘以-1, 方块表示乘法. 可以看出算法分为6个层.

输入为f[N] = {f0, f1, f2, f3, f4, f5, f6, f7};

输出为F[N] = {F0, F1, F2, F3, F4, F5, F6, F7};

 

那么AAN算法的伪代码描述为 :

 

// 第一层. 有一点FFT蝶型算法的味道. 第一层是以原始数据作为输入的.

tmp1_0 = f0 + f7;

tmp1_1 = f1 + f6;

tmp1_2 = f2 + f5;

tmp1_3 = f3 + f4;

tmp1_4 = f3 + (-1)*f4;   // 就是减法, 可以认为箭头是减法的.

tmp1_5 = f2 + (-1)*f5;

tmp1_6 = f1 + (-1)*f6;

tmp1_7 = f0 + (-1)*f7;

 

// 第二层, 是以第一层输出作为输入的

tmp2_0 = tmp1_0 + tmp1_3;

tmp2_1 = tmp1_1 + tmp1_2;

tmp2_2 = tmp1_1 + (-1)*tmp1_2;

tmp2_3 = tmp1_0 + (-1)*tmp1_3;

 

tmp2_4 = (-1)*tmp1_4 + (-1)*tmp1_5;

tmp2_5 = tmp1_5 + tmp1_6;

tmp2_6 = tmp1_6 + tmp1_7;

tmp2_7 = tmp1_7;

 

// 第三层, 是以第二层输出作为输入的

tmp3_0 = tmp2_0 + tmp2_1;

tmp3_1 = tmp2_0 + (-1)*tmp2_1;

tmp3_2 = tmp2_2 + tmp2_3;

tmp3_3 = tmp2_3;  // 这些可以不需要的

tmp3_4 = tmp2_4;

tmp3_5 = tmp2_5;

tmp3_6 = tmp2_6;

tmp3_7 = tmp2_7;

 

// 第四层, 是以第三层的输出作为输入的

tmp4_0 = tmp3_0;

tmp4_1 = tmp3_1;

tmp4_2 = tmp3_2*a1;

tmp4_3 = tmp3_3;

tmp4_4 = tmp3_4*a2*(-1) + (tmp3_4+ tmp3_6)*a5*(-1);

tmp4_5 = tmp3_5*a3;

tmp4_6 = tmp3_6*a4 + (tmp3_4 + tmp3_6)*a5*(-1);

tmp4_7 = tmp3_7;

 

// 第五层

tmp5_0 = tmp4_0;

tmp5_1 = tmp4_1;

tmp5_2 = tmp4_2 + tmp4_3;

tmp5_3 = tmp4_2*(-1) + tmp4_3;

tmp5_4 = tmp4_4;

tmp5_5 = tmp4_5 + tmp4_7;

tmp5_6 = tmp4_6;

tmp5_7 = tmp4_5*(-1) + tmp4_7;

 

// 第六层, 也该到终点了, 寡人都累坏了 :)

tmp6_0 = tmp5_0;

tmp6_1 = tmp5_1;

tmp6_2 = tmp5_2;

tmp6_3 = tmp5_3;

tmp6_4 = tmp5_4 + tmp5_7;

tmp6_5 = tmp5_5 + tmp5_6;

tmp6_6 = tmp5_5 + (-1)*tmp5_6;

tmp6_7 = (-1)*tmp5_4 + tmp5_7;

 

// 功德圆满

F[0] = tmp6_0;

F[4] = tmp6_1;

F[2] = tmp6_2;

F[6] = tmp6_3;

F[5] = tmp6_4;

F[1] = tmp6_5;

F[7] = tmp6_6;

F[3] = tmp6_7;

 

伪代码怎么看怎么像FFT蝶形算法.

这个描述比较中规中矩, 中间有一些空心圆, 是可以直接优化的. 不需要这么多临时变量的. 为了使算法描述的更新清晰和AAN运算流更一致, 我觉得规范的AAN算法应该要写出这些”不必要”的部分.

 

下面给出JEPG编码里头的AAN部分, 这个代码是经过优化的.

优化思想有2个方面.

1. 空心点不在参与到代码中;

2. 完成第一阶后, 之后的阶都是很明显个行其事的, tmpi_0到 tmpi_3是相关的, tmpi_4到tmpi_7是相关的. 很明显的分为了两条路. (其中i表示层, i>1)

    /*

     * AAN算法第一阶变换

     */

    tmp0 = dataptr[0] + dataptr[7];

    tmp7 = dataptr[0] - dataptr[7];

    tmp1 = dataptr[1] + dataptr[6]; // 低频

    tmp6 = dataptr[1] - dataptr[6]; // 高频

    tmp2 = dataptr[2] + dataptr[5];

    tmp5 = dataptr[2] - dataptr[5];

    tmp3 = dataptr[3] + dataptr[4];

    tmp4 = dataptr[3] - dataptr[4];

   

    /* Even part 偶数部分 */

    tmp10 = tmp0 + tmp3; /* phase 2*/

    tmp13 = tmp0 - tmp3;

    tmp11 = tmp1 + tmp2;

    tmp12 = tmp1 - tmp2;

   

    dataptr[0] = tmp10 + tmp11; /* phase 3 */

    dataptr[4] = tmp10 - tmp11;

   

     z1 = (tmp12 + tmp13) * ((double) 0.707106781);

    dataptr[2] = tmp13 + z1; /* phase 5 */

    dataptr[6] = tmp13 - z1;

   

    /* Odd part 奇数部分 */

    tmp10 = tmp4 + tmp5; /* phase 2 */

    tmp11 = tmp5 + tmp6;

    tmp12 = tmp6 + tmp7;

   

    /* The rotator is modified from fig 4-8 to avoid extra negations. */

    z5 = (tmp10 - tmp12) * ((double) 0.382683433); /* c6 */

    z2 = ((double) 0.541196100) * tmp10 + z5; /* c2-c6 */

    z4 = ((double) 1.306562965) * tmp12 + z5; /* c2+c6 */

    z3 = tmp11 * ((double) 0.707106781); /* c4 */

      

    z11 = tmp7 + z3;     /* phase 5 */

    z13 = tmp7 - z3;

      

    dataptr[5] = z13 + z2;   /* phase 6 */

    dataptr[3] = z13 - z2;

    dataptr[1] = z11 + z4;

    dataptr[7] = z11 - z4;

   

   

为什么?

只知其然不知其所以燃是不对的, 对于AAN算法的理论推导我很困惑. 只是模糊的知道一点点, 很可能不对.一并记录如下 :

1. 信号是具有频率的, 有高频低频, 不管是图象还是Audio, 人类都是对低频部分是比较敏感, 那么处理JPEG图象或者Audio的时候都是希望能让低频部分保留多一些, 高频部分过滤多一些.

2.

 

 

参考文献

<<Implementation of Motion-JPEG Using an ASIC Prototype Board>>  Master’s Thesis