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
注:本文转自 http://blog.csdn.net/sshcx/article/details/1674224