The T. Agui Y. Arai and M. Nakajima DCT Algorithm (AAN DCT算法)

原文地址:http://www.stefan-kuhr.de/download/thesis2side.pdf

DCT算法是数字图像压缩中非常基本的算法,数字图像压缩时使用的是二维N*N(JPEG中N一般取8)的DCT.定义如下:

$F(\mu,\nu)=\frac{c\left(\mu\right)}{2}\frac{c\left(\nu\right)}{2}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}f(n,m)\cos\left\{\frac{\left(2n+1\right)\mu\pi}{2N}\right\}\cos\left\{\frac{\left(2m+1\right)\nu\pi}{2N}\right\}\left(\mu,\nu=0,1\dotsN-1\right)$

$c\left(i\right)=\begin{cases}\frac{1}{\sqrt2}&i=0\\1&i\ne0\end{cases}$

如果按公式直接计算,那么计算量是非常巨大的,对于N=8,变换结果有8*8个系数,每个系数都需要进行8*8次乘法和8*8-1次加法。因此总的计算量为8*8*8*8次乘法;8*8*(8*8-1)次加法.实际中的处理是把二维8*8的DCT按 $\mu,\nu$ 转换为按行和按列的两个8次一维8点DCT计算.而8点DCT只需要8*8次乘法加8*(8-1)次加法.8行(列)DCT总共需要8*8*8次乘法加8*8*(8-1)次加法.因此总的计算量为2*8*8*8次乘法,2*8*8*(8-1)次加法,仅为直接计算的1/8.

对于8点DCT $F(\mu)=\frac{c\left(\mu\right)}{2}\sum_{n=0}^{N-1}f(n)\cos\left\{\frac{\left(2n+1\right)\mu\pi}{2N}\right\}$ .,AAN算法 (T. Agui Y. Arai and M. Nakajima. A fast DCT-SQ Scheme for Images. Transactions of the IEICE, E 71:1095 – 1097, Nov 1988.)只需要13次乘法,29次加法.而且其中8次乘法都发生在计算的最后一步.对于JPEG标准DCT计算后紧跟着就是对DCT系数进行量化,这8次乘法可以合并在量化过程中.因此对于JPEG来说,AAN算法实际只需要5次乘法,29次加法. 其流程图如下

带箭头的连线表示减,不带箭头的边表示加,a1~a5为乘法的常数.

网上对于AAN算法的介绍不是很多,而且原始论文T. Agui Y. Arai and M. Nakajima. A fast DCT-SQ Scheme for Images. Transactions of the IEICE, E 71:1095 – 1097, Nov 1988 也始终找不到.对于其由来一直很不明白.很多论文现其都是直接对变换矩阵进行分解得到最终算法.最近找到一篇论文,也是用矩阵分解的方法对AAN算法进行了推导.虽然不知道为什么要进行矩阵分解,但由于其分解过程还比较详细.因此翻译记录于此,供日后慢慢领悟.

AAN算法的思路是借助16点实对称序列的WFTA计算8点DCT,原理如下:

2N点实对称序列: $f(n)=f(2N-1-n)$ 的DFT为:

$F(\mu)=\sum_{n=0}^{2N-1}f(n)e^{-j\frac{2n\mu\pi}{2N}}=\sum_{n=0}^{N-1}f(n)e^{-j\frac{2n\mu\pi}{2N}}+\sum_{n=N}^{2N-1}f(2N-1-n)e^{-j\frac{2n\mu\pi}{2N}}$

上式第二项中进行换元k=2N-1-n并利用e的周期性得:

$F(\mu)=\sum_{n=0}^{2N-1}f(n)e^{-j\frac{2n\mu\pi}{2N}}=\sum_{n=0}^{N-1}f(n)e^{-j\frac{2n\mu\pi}{2N}}+\sum_{k=0}^{N-1}f(k)e^{j\frac{2(k+1)\mu\pi}{2N}}$

将上式第二项中k带回n得:

$F(\mu)=\sum_{n=0}^{2N-1}f(n)e^{-j\frac{2n\mu\pi}{2N}}=\sum_{n=0}^{N-1}f(n)e^{-j\frac{2n\mu\pi}{2N}}+\sum_{k=0}^{N-1}f(n)e^{j\frac{2(n+1)\mu\pi}{2N}}&\dots\left(1\right)$

(1)式两边同乘 $e^{-j{\frac{\mu\pi}{2N}}}$ 得

$F\left(\mu\right)e^{-j{\frac{\mu\pi}{2N}}}=\sum_{n=0}^{N-1}f(n)e^{-j\frac{\left(2n+1\right)\mu\pi}{2N}}+\sum_{k=0}^{N-1}f(n)e^{j\frac{\left(2n+1\right)\mu\pi}{2N}}=2\sum_{n=0}^{N-1}f\left(n\right)\cos\frac{\left(2n+1\right)\mu\pi}{2N}$

因此2N点实对称序列的前N个DFT系数乘上一个复数比例因子恰好就是N点实序列的DCT系数。正是利用这一关系,借助16点WFTA算法推导出了AAN算法

16点WFTA算法

DFT算法写成矩阵形式就是:

对N=16,记 $W=e^{\frac{-j\pi}{8}}$

则变换矩阵M为

其偶数行前后部分相等,奇数行数前后半部分全部相反。因此将M按如下规则重排: 所有偶数行相对位置不变放到前8行,所有奇数行相对位置不变放到后8行.由于奇数行和偶数行的前后数据关系,将重拍后的矩阵看成四个8*8的分块矩阵则: $M=P_1*M^\prime*R_1=P_1*\begin{bmatrix}x&x\\y&-y\end{bmatrix}*R_1$

$\begin{bmatrix}x&x\\y&-y\end{bmatrix}=P_1*\begin{bmatrix}x&0\\0&y\end{bmatrix}*\begin{bmatrix}I_8&I_8\\I_8&-1_8\end{bmatrix}$

再加上奇偶行的重排矩阵得:

对于 $图片$ 左上角的8*8方阵满足同样的对称关系,因此可以继续进行同样的分解,但右下角的8*8方阵不满足这样的对称关系.只能专门进行分解.暂记为 $图片$ 的结构为 $图片$ ,要想对 $图片$ 的左上角进行变换,而右下角保持不变,则变换矩阵 $图片$ 的结构应为 $图片$ 分别为变换矩阵和单位阵.经过再次分解后: $图片$ 为

同样对这次分解后剩下的不具备行对称关系4*4方阵先记为 $图片$ .这次分解后剩下的左上角4*4方阵任然具有相似的对称偶数行数据前后相等的关系.因此分解可以继续进行.并且将 $图片$ 进行分块.左上角为一个4*4方阵,右下角为12*12方阵，其它位置为0。我们只希望变换左上角4*4方阵,而保持其它部分不变,则很容易得出变换矩阵 $图片$ 的结构应该为 $图片$ .

经过再次分解后 $图片$ 为

: $\M=P_1*P_2*P_3*M^\prime^\prime^\prime*R_3*R_2*R_1$

为了使 $图片$ 左上角的2*2方阵对角化,再对 $图片$ 进行一变换,左乘一个, $图片$ 得, $图片$ 为

变换矩阵M左边的矩阵 $图片$ 都是重排列行序用的。可以把每一个 $图片$ 看成 $图片$ 的一个置换,这些 $图片$ 的乘积还是一个对行的排列,不影响计算结果.因此可以用一个 $图片$ 表示

因此最终 $M=P*\overline{M}*R_4*R_3*R_2*R_1$

其中 $图片$ 中 $图片$ ， $图片$ 分解如下(分解规则：只使用初等行变换,化成只含0，1的一些列矩阵和一个稀疏矩阵之积):

令 $\gamma=\frac{\pi}{8}$

备注:

$图片$ 的分解结果是4*4,8*8矩阵,运算的时候需要把分解的结果带入到原16*16矩阵中.方法如下.：将原方阵分块后看成一个对角阵,把 $图片$ 阶方阵也按同样分块,并把分解结果填入到对应位置中即可

例如要将 $图片$ 的分解结果带回去,已知 $图片$ 则将 $图片$ 分成8*8块的分块矩阵: $图片$

使用WFTA计算DFT所需的计算量复数加法可由 $图片$ 以及 $图片$ 的分解式中每行多余1个正负1的个数确定。复数乘法则由以上矩阵中非1的个数决定,16点WFTA只需要10次复数乘法和74次复数加法.

使用WFTA构建AAN DCT

16点WFTA变换矩阵形式为 $图片$ .其中输入 $图片$ 要求是16点的.但由于计算DCT时 $图片$ 满足实对称,因此只要其中前面8点确定了,后面8点也就确定了.可以在 $图片$ 与 $图片$ 之间增加一个映射,当取的数据在后8点时,自动映射成正确的前面8点中的某个位置.

U正是这样的映射. 因此 $图片$ 前8点数据即为 $图片$ DCT结果(相差一个比例因子) $图片$ 都为上节WFTA中的对应矩阵. $图片$ 和 $图片$ 可以合并到一起,记为 $图片$ .

因为 $图片$ 中只有 $图片$ 含有复数,因此不会出现两个复数相乘的情况.而整个结果又是实数,因此可以完全忽略 $图片$ 里复数的虚部,可以直接将对应复数虚部置0

由于只需要计算Y[0]~~Y[7]这8个数据,因为可以将以上各矩阵进一步简化以便减少计算.由矩阵P可知,Y[0],Y[2],Y[4],Y[6] 由 $图片$ 的0,4,2,5行计算出.Y[1],Y[3],Y[5],Y[7]由 $图片$ 的8,9,10,11行计算出(由矩阵P第i行中1的位置指明使用 $图片$ 的哪一行计算.)

对于0,4,2,5行.每行的后8点数据都是0,因此这4行可表示成 $图片$

将矩阵 $图片$ 分别进行适当的分块后可以把它们都看成对角阵.因此 $图片$ 从左到右的每一步之积也都是对角阵.Y[0],Y[2],Y[4],Y[6]只跟结果矩阵的0,4,2,5行有关,位于矩阵的上部分,而矩阵每一步之积顶部的8行中,后面的一半数据都全为0.因此计算时可以直接取 $图片$ 的左上角8*8小块.因此Y[0],Y[2],Y[4],Y[6]的计算可以简化为:

$图片$ 分别为 $图片$ 的左上角的8*8小块.

其中
最后两行的特殊关系还可进一步分解(交换倒数第一列和倒数第三列,再将倒数第二行的负一倍加到倒数第一行可消去最后一行的两个非零常数)

为了进一步减少运算量,还需对 $R_4^\prime,R_3^\prime,R_2^\prime$ 进行简化.

$R_1^\prime^\prime=\begin{bmatrix}1&0&0&0&0&0&0&1&\\0&1&0&0&0&0&1&0&\\0&0&1&0&0&1&0&0&\\0&0&0&1&1&0&0&0&\\0&0&1&0&0&1&0&0&\\0&1&0&0&0&0&1&0&\\0&1&0&0&0&0&1&0&\\1&0&0&0&0&0&0&1\end{bmatrix}=\begin{bmatrix}I_4&I_4^\prime\\I_4^\prime&I_4\end{bmatrix}$

其中

$I_4=\begin{bmatrix}1&0&0&0&\\0&1&0&0&\\0&0&1&0&\\0&0&0&1\end{bmatrix}I_4^\prime=\begin{bmatrix}0&0&0&1&\\0&0&1&0&\\0&1&0&0&\\1&0&0&0\end{bmatrix}$

因为有 $I_4^\prime*I_4^\prime=I_4$ 所以有:

$\begin{bmatrix}x&y&\\m&n\end{bmatrix}*R_1^\prime^\prime=\begin{bmatrix}x*I_4+y*I_4^\prime&x*I_4^\prime+y*I_4&\\m*I_4+n*I_4^\prime&m*I_4^\prime+n*I_4\end{bmatrix}=\begin{bmatrix}x*I_4+y*I_4^\prime&0&\\m*I_4+n*I_4^\prime&0\end{bmatrix}*R_1^\prime^\prime$

由以上关系得:

$图片$ 可简化为: $图片$

如此简化后 $R_1^\prime^\prime$ 下半部分可以全部用0代替

同理 $R_3^\prime,R_4^\prime$ 可分别简化为:

$R_3^\prime^\prime=\begin{bmatrix}1&1&0&0&0&0&0&0&\\1&1&0&0&0&0&0&0&\\1&-1&0&0&0&0&0&0&\\-1&1&0&0&0&0&0&0&\\0&0&0&0&1&0&0&0&\\0&0&0&0&0&1&0&0&\\0&0&0&0&0&-1&0&0&\\0&0&0&0&-1&0&0&0\end{bmatrix}$

并且 $R_2^\prime^\prime$ 的2,3,6,7行可用0代替

$R_4^\prime^\prime=\begin{bmatrix}2&0&0&0&0&0&0&0&\\0&0&0&0&0&0&0&0&\\0&0&1&0&0&0&0&0&\\0&0&-1&0&0&0&0&0&\\0&0&0&0&1&0&0&0&\\0&0&0&0&0&1&0&0&\\0&0&0&0&0&-1&0&0&\\0&0&0&0&-1&0&0&0\end{bmatrix}$

并且 $R_3^\prime^\prime$ 的1,3,6,7行可用0代替

与 $R_4^\prime^\prime$ 相邻的矩阵

$Q^\prime=\begin{bmatrix}1&0&0&0&0&0&0&0&\\0&0&1&0&0&0&0&0&\\0&0&0&0&1&0&0&0&\\0&0&0&0&0&1&0&-1\end{bmatrix}$

由于 $图片$ 中第四行和第七行相差一个-1,这个-1可以合并到 $图片$ 中去。合并后的 $图片$

最后去掉全为0的行列后得到

其流程图为:

同理,计算基数下表的系数Y[1],Y[3],Y[5],Y[7]时,可以直接采用对应矩阵的右下角8*8方阵进行计算. $R_4^\prime,R_3^\prime,R_2^\prime$ 对应右下角8*8方阵都是单位阵.因此只需要化简

即可

$图片$ 是 $图片$ 的右下角8*8分块.

令 $图片$

V可分解为:

V的分解式中最后一个矩阵

$图片$ 与 $图片$ 相乘.根据 $图片$ 中每列数据的对称关系,可以修改V的最后一个分解式.将 $图片$ 每列中两个正负1合并到V的最后一个分解式中去.{合并要求矩阵乘积的结果不变,及V最后一个分解式的修改后每行与 $图片$ 中每列对应数据相乘的结果都不能改变.记V最后一个分解式一行数据为a[i], $图片$ 每列数据为b[i],因此对V最后一个分解式每一行的修改都必须满足 $图片$ 保持不变.因此只能利用 $图片$ }的关系修改a(i) :i位为1时,如果7 - i mod8 的位置为0,则可以将i位的1移动到7 - i mod 8的位置,并改变符号.

因此修改后的V的最后一个分解式为

$\begin{bmatrix}1&0&0&0&0&0&0&0&\\0&0&1&0&0&0&0&0&\\0&1&0&0&0&0&0&0&\\1&1&0&0&0&0&0&0&\\0&0&-1&-1&0&0&0&0\end{bmatrix}$

去掉全为0的行列最后可写成

对应流程图为；

寄偶部分下标合在一起即为文章开头部分的整个流程图

...............................................................................完.............................................................................................................

.................................................................2011年3月22号凌晨1点40分初稿................................................................................/

posted on 2011-04-24 13:11 TinyCoder 阅读(1921) 评论(0) 编辑收藏举报