二维DCT的优化算法

-11-302016-11-30姓名：学号：姓名：学号：院系：电子与信息工程学院指导教师：《音视频编码技术》课程报告《2-DDCT算法优化》摘要 离散余弦变换（DCT）是一种正交变换，它本身并不会改变信息源的熵值，因而是一种无损的变换。但是由于离散余弦变换（DCT）变换能够给量化和编码等过程创造很好的条件，因而在数字信号处理领域，特别是语音和图像压缩领域有着很广泛的应用。对于语音信号与图像信号这样相关性很强的信号而言，离散余弦变换（DCT）接近于最优的KL变换(Karhunen-LoeveTransform)。传统的2-D离散余弦变换（DCT）变换过程是行-列变换，这种方法直观简单但是效率低下，在有些实时性要求较高的图像压缩领域并不能满足要求。从离散余弦变换（DCT）的数学表达式上可以看出，其具有很高的对称性，利用这种对称性可以极大地降低算法的复杂度。本次课程报告在总结前人研究的基础上，展示并实现了一种快速变换和反变换算法，并在MATLAB平台上进行了实现与测试。结果表明，这种算法具有高效、简单、易于实现的特点。关键词：快速离散余弦变换，快速离散余弦反变换，快速算法，图像压缩

2-DDCT算法优化引言研究背景 在本次的课程报告中，DCT的优化算法主要是针对于2维输入数据，因而此算法主要针对的是图像处理。 在图像处理领域中，图像压缩编码依据原理可以分为熵编码、预测编码、变换编码和混合编码。DCT属于其中的变换编码，变换编码还包括傅里叶变换和沃尔什变换等等。变换编码是一种数学变换对，通常是将2维空间域上的图像经过正交变换映射到另一个空间域上（例如频域），并且使得变换后的系数之间的相关性降低。就变换本身而言，它是一种在数学上无损的变换，通过反变换可以恢复原信号。但是这种变换的特点在于，变换后图像的大部分能量只集中在少数几个变换系数上，这就为量化编码和熵编码带来了很大的便利，在对系数进行编码后就能实现图像的压缩。 在理论上，K-L变换是最优的正交变换，它能够完全消除图像块内像素间的线性相关性，经过K-L变换后各个系数在统计上不相关，其协方差矩阵是对角阵，因而K-L变换能够大大的减小原数据的冗余度。如果在K-L变换之后丢弃数值娇小的系数，所造成的均方误差是所有正交变换中最小的。但是K-L变换有一个很大的缺点，它的基向量是图像的协方差矩阵的特征向量。这就意味着对于不同的图像，它的基向量是不相同的，为了对图像进行K-L变换，就需要先计算每一幅图像的相关矩阵，再对相关矩阵作特征分解，取特征向量作为变换的基向量。这就使得K-L变换的使用变得很麻烦，甚至是难以实现。 DCT变换的性能仅次于K-L变换，被认为是一种准最佳变换，对于图像这种相关性较大的数据，DCT的性能接近于K-L变换。此外，DCT的变换矩阵是固定的，与图像无关。而且它是构造成对称的数据序列，避免了子图像边界处的跳跃和不连续现象。 傅里叶变换是应用最早的变换之一，并且也有快速算法。但是傅里叶变换存在一个很明显的缺陷，恢复出来的边界会存在不连续现象，这使得恢复的图像存图SEQ图\*ARABIC1离散余弦变换（DCT）效果图在明显可见的方格，影响图像质量。 沃尔什变换比DCT简单，实时性更高，但是性能比DCT略差。 基于以上的特点，DCT具有广泛的应用，在常见的JPEG静态图像编码以及MJPEG、MPEG动态图像编码等标准中都有应用。这种应用的广泛也对DCT的算法效率提出了更高的要求，对DCT快速算法也有着更加迫切的需求。前人研究 由于DCT在语音和图像压缩领域的广泛应用，已经有人提出了许多快速算法和快速计算DCT的VLSI架构。DCT的传统计算方法是行-列法，先对行计算1-DDCT，再对列计算1-DDCT。对于N×N大小的图像而言，这种算法需要计算2N个1-DDCT，算法的复杂度为O(N2)。因而为了更加高效地进行DCT运算，有人提出了直接对整个2维矩阵做DCT运算。目前，最为高效的DCT算法是由Duhamel和Guillemot在1990年的一篇论文中提出的多项式变换计算方法。在这种方法中，2-DDCT运算的乘积运算相对于传统的方法而言减小了50%。还有人提出基于FFT的算法，这种算法的复杂度在传统方法的50%到75%之间。此外，Cho和Lee提出了专门针对于矩阵大小为4×4快速DCT算法，这种算法仅仅局限于大小为4×4的变换，如果要进行扩展更大的维数 在本次的课程报告中所展示的是一种有别于以上快速算法的方法，这种方法和Cho与Lee提出的算法有异曲同工之处，本质上来说是对他们算法的一种扩展。这种算法对于大小为N×N的矩阵而言，仅仅需要计算N个1-DDCT，但是要求N=2m,m∈R。算法需要的乘积运算次数仅仅为传统算法的50%，与Duhamel和Guillemot的算法有着相同复杂度。但是相比于Duhamel和Guillemot报告的章节安排本次报告的主要内容和章节安排如下：快速DCT算法的理论分析 从数学上对算法进行理论分析，相关公式推导，总结出最终的算法结构。相应的快速IDCT算法的理论分析基于前面的理论推导，总结出快速IDCT算法的结构。与其他快速算法的复杂度比较将此算法与传统的算法、其他算法的复杂度作数值上的比较。算法的实现与测试在前面理论分析的基础上，在MATLAB平台上实现4×4的快速算法，并对图像做DCT和IDCT变换，检验算法的有效性，同时与MATLAB自带的DCT和IDCT函数的运算时间进行比较。总结与展望分析此算法的优缺点与可改进的空间。理论分析2-DDCT理论分析 对于一个给定的大小为N×N的二维矩阵x假设二维DCT矩阵为YYDCT变换公式为：Ymn=cmc其中，cm=定义ymn为Ymn对ymn=i=0N-1j=0N-1xijYmn=cm利用积化和差公式，可以进行如下的变换，cos2i+1mπ2Ncos将（2.3）带入到（2.2.a）中，得到ymn+i=0form,n=1,2,⋯,N-1(2.4)同时，我们定义两个新的变量AAmn=i=0N-1j=0N-1Bmn=i=0N-1j=0N-1x所ymnymn=12(从（2.5）与（2.6）中，可以发现，二维的DCT变换可以分解为两个分量的叠加，并且这两个分量（Amn与Bmn）的形式类似于一维DCT的形式。那么下面的目标就在于如何将Amn与Bmn转换成一维DCT的形式。从而达到对于N×N的矩阵而言，只需要计算N次一维 一维DCT的表达式为Y与（2.5.a）和（2.5.b）对比可以发现，如果{2i+1m±2j+1n}可以转换为{2i+1k, 由于{i,j=0,1,2,⋯,N-1}，那么{(2j+1)}的值要么等于{(2i+1)}整数倍除以2N的余数2j+1=p2i+1mod2N或者2j+1=2N-p2i+1mod2N其中p=1,3,5,7,⋯,N-1，如果将（2.7）进行化简，就可以得到j=(pi+p-12)mod2N或者j=(N-1-pi-p-12)mod2N其中p= 从（2.8）中可以发现，如果i取0到(N-1)的不同值，那么对应得到的j的序列是不一样的，这样就得到了N2个不同的(i,j)组合，这说明，依据（ 此外每一个i的值都对应着一个j的序列，这意味着，依据（2.8）可以对输入数据分成N组，每一个组内元素的下标(i,j)都满足（j(p|a)=(pi+p-12或者j(p|b)=(N-1-pi-p-12其中，p=1二维输入数据{xi,j,i,j=0,1,2,⋯,N-1}通过（2.9）就可以分成或者{用Rpa和RR(2.10)为了进一步进行简化，假设已知pi+p-12除以N的商为qpiRR(2.10)按照（2.10）中的分组，Amn与BAmn=pN-1{Bmn=pN-1S其中Tpam,nTpbm,n=xSpam,n=xSpbm,n=x把（2.11）带入到（2.6）中，可以得到，ymn=p为了把ymn表示为一维DCT的和的形式，需要将Tpam,n、TpbT将上面的式子拆分为n为奇数和偶数两种情况：当n为奇数时T当n为偶数是T那么以同样的方式，可以获得TpbTSS把（2.14）-（2.17）带入到（2.13）中，可以得到ymn的新表达式y可以发现，i=0i=0i=0i=0（2.19.a）和（2.19.b）分别是{xijpa+xijpb}序列的一维DCT变换序列的第m+np和m-np个元素；（2.19.c）和（2.19.d）分别是{-1qpixijpafg那么（2.19.a）和（2.19.b）就分别对应着±fpl,l=0,1,2,⋯,N-1；（2.19.c）和（2.19.d）就分别对应着±gpl,l=0,1,2,⋯,N-1。这就意味着，大小为N×N 从（2.18）可以看出，要计算出ymn，可以先计算出fpl和gpl，然后再将相对应的fpl和gpl进行线性叠加，得到最终结果cos如果把fpl表达为f那么把（2.22）带入就可以得到，± 同样的方式，当n为奇数时，g±g上面的两个式子意味着，fpl总是会以±f(N-p)l的形式出现。同样，gpl 基于以上的分析，假如输入矩阵的大小为8×8，那么这个算法就可以用以下的图来表示，图SEQ图\*ARABIC28X8DCT算法结构图，O表示相加，虚线表示取反图2（续）8X8DCT算法结构图，O表示相加，虚线表示取反2-DIDCT理论分析 对于正交变换而言，如果不考虑系数的影响，那么其相对应的反变换就直接是正变换的反过程而已。那么可以利用这个来推导出快速IDCT算法。 如果在DCT变换中不考虑（2.1）中前面的系数，那么反变换的算法结构就仅仅是正变换算法结构的取反。例如对于8X8的IDCT变换而言，如果不考虑系数，那么其算法结构如图3所示。可以看到，相对于图2而言，图3仅仅是将图2中信号的流向做了取反，相当于是反向做运算，然后做N个一维IDCT，其运算量与（一）中分析的DCT的运算量相当，相对于传统的IDCT算法而言也仅仅需要50%的一维DCT变换次数。图SEQ图\*ARABIC38X8IDCT算法结构图，O表示相加，虚线表示取反图SEQ图\*ARABIC4（续）8X8IDCT算法结构图，O表示相加，虚线表示取反 但是如果将系数考虑进来，就需要做一些额外的调整，调整的过程如下：对于初始输入数据对相对应的DCT系数做归一化处理将第一步得到的数据乘以系数1每个一维IDCT的第一个输入数据需要再乘以系数1经过以上三个步骤，处理后就可以得到相对应的IDCT数据。算法的复杂度分析 从（一）和（二）部分的内容中可以看到对于N×N为的DCT变换，仅仅需要计算N个一维的DCT。我们在本小节中以需要的乘法运算的次数和加法运算次数来作为算法复杂度比较的标准 由于N×N大小的二维DCT变换仅仅需要N次一维的DCT变换，那么相对于传统的二维DCT算法而言，本报告中的算法所需要的乘法运算量就仅仅是传统方法的50%。而对于一维DCT的算法复杂度，目前比较高效的算法需要的N如果在二维DCT变换中使用这个一维DCT算法，那么所需乘法运算量为，N表格SEQ表格\*ARABIC1所需乘法运算次数比较矩阵大小传统算法[1][2][3]本报告算法4X432241616168X8192144104969616X16102476856851251232X3251203840274025602560表格SEQ表格\*ARABIC2所需加法运算次数比较矩阵大小传统算法[1][2][3]本报告算法4X472727068748X846446446248446616X162592259225582531253032X321337613376129501257812738 至于所需的加法运算量，从前面的分析中可以知道，除了DCT运算的加法运算量，其他地方加法运算量为N21+log2N-2N-(N-2)。而对于一维DCT变换而言5 依据（2.25）和（2.26）和其他参考文献中的快速算法，列出了表格1和表格2，分别展示了各个算法所需的乘法运算次数和加法运算次数。从表格中可以发现，本报告中的算法所需的乘法运算量与文献[3]中一样，并且比其他算法都要低。同时，加法运算量与其他算法基本一致。实现与测试 在本章节中，将根据第二章节中的DCT和IDCT快速算法的理论分析，在MATLAB平台上实现相应的DCT变换算法和IDCT变换算法。然后利用这两个算法对一个实际图像做变换与反变换处理，观察处理的效果。算法的实现 在本小节里面，将会依据第二章节的分析，在MATLAB平台上实现4X4的快速DCT与IDCT算法。 从第二章中的分析中可以知道，快速算法可以用流向图很简洁地表示出来。那么依据前面的公式，对于4X4的DCT和IDCT变化可以得到以下的流向图，图SEQ图\*ARABIC54X4算法结构图，（a）DCT（b）IDCT那么依据图5，算法的流程图如下：图SEQ图\*ARABIC6快速DCT算法流程图图SEQ图\*ARABIC7快速IDCT算法流程图算法的测试 测试环境：Dell商用台式机，MATLAB2013a，测试图片大小为614x614.图SEQ图\*ARABIC8快速DCT与IDCT测试效果图 从图8的实际测试结果来看，原图经过了4X4的DCT变换之后，在每个4X4大小的方格内，图像的能量都集中分布于少数几个系数上，这符合正交变换的特点。同时经过了IDCT变换之后，可以完整的恢复图像原来的模样，且不存在失真与格纹现象。这证明了算法的正确性。表格SEQ表格\*ARABIC3程序运行时间比较自定义函数MATLAB已有函数函数名称调用次数运行时间函数名称调用次数运行时间main13.969smain16.804sfdct2237161.594sdct2237163.228sfidct2237161.396sidct2237162.674s 表格3是MATLAB利用CPU的运行时间计算出来的结果，虽然结果的具体数值并不具有很高的可靠性，但是在相同条件下得出来的结果，相互之间具有可比性。通过比较可以发现，在相同的调用次数的条件下，无论是fdct2函数还是fidct2函数，其运行时间都比MATLAB已有函数的运行时间减少了将近50%。总结与展望 本次课程报告首先分析了离散余弦变换（DCT）在图像压缩领域的重要性，它本身虽然不具有压缩能力，但是变换给后续的编码和量化操作带来了便利，在多个图像、视频标准中都有应用。DCT的广泛应用也为人们对它的研究增添了动力，由于传统计算方法的效率低下，许多文献提出了相关的DCT快速算法。本次报告就研究、实现和分析了其中一种性能比较优越的快速算法。这种快速算法的乘法运算量仅仅只有传统算法的50%，而加法运算量和其他方法相当，具有很高的效率，同时，这种算法可以能够用流向图来直观的表示，简单而直观。而且由于这种算法具有比较好的对称性，它能够很容易地应用到多线程并行处理中，这样效率又可以得到很大的提升，同时这种算法还可以直接用VLSI结构来实现，这说明这种算的应用相对简单。在进行了理论分析之后，基于理论分析在MATLAB平台上实现了此算法，并用实际图片进行了测试，测试结果表明，此算法可以快速地进行DCT变换和IDCT变换，能够完整的恢复出原图。此外，还将此算法与MATLAB自带的DCT变换和IDCT变换函数在相同的条件下进行了运行时间比较，比较结果表明，此算法的运行时间相对于MATLAB自带的函数减少了将近50%，具有较高的效率。 虽然这种的算法具有很高的效率，但是它依然有它的局限性。这种算法的最大的局限性在于，如果二维矩阵的维数很大，那么这种算法将会难以实现。因为这种算法需要对输入数据进行排列，在进行了一维DCT后还需要进行多次的线性运算，这种排列与线性运算的组合并不具有很明显的规律性。在算法的实现过程中，对于16X16大小的矩阵，需要列出256个式子；但是对于32X32的矩阵，需要的式子达到了1024个。所以对于大型的矩阵，这种方法是不适合的，比较适合于16X16以下的矩阵。参考文献[1] C.Ma,“Afastrecursivetwodimensioncosinetransform”,IntelligentRobotsand ComputerVisio:SeventhinaSeries,DavidP.Casasent,Ed.,inProc.SPIE1002, pp.541-548,1988.[2] M.Vetterli,“Fast2-Ddiscretecosinetransform,”inProc.ICASSP’85,Mar.1985. P.DuhamelandC.Guillemot,“Polynomialtransformcomputationof2-DDCT,”in Proc.ICASSP’90,pp.1515-1518,Apr.1990.[3] N.I.ChoandS.U.Lee,“Afast4X4DCTalgorithmfortherecursive2-DDCT,”IEEE Trans.Acoust.,Speech,SignalProcessing,submittedforpublication.[4] F.A.KamangarandK.R.Rao,“Fastalgorithmforthe2-Ddiscretecosinetansform,” IEEETrans,Comput.,vol.C-31,pp.899-906,Spet.1982.[5] B.G.Lee,“Anewalgorithmtocomputethediscretecosinetransform,”IEEETrans. Acoust.,Speech,SignalProcessing,vol.ASSP-32,pp.1243-1245,Dec.1987.[7] 杨帆.数字图像处理与分析.北京：北京航空航天大学出版社,2007[8] 全红艳，曹桂涛.数字图像处理原理与实现方法.北京：机械工业出版社，2014

附件函数fdct2源码function[Y]=fdct2(X)%快速4X4DCT变换函数%%By黄跃辉2016/11/30%%%检测输入量[a,b]=size(X);ifa~=4||b~=4Y=zeros(4,4);return;end%%A=zeros(4,4);B=zeros(4,4);A(1,1)=X(1,1)+X(1,4);A(2,1)=X(2,2)+X(2,3);A(3,1)=X(3,3)+X(3,2);A(4,1)=X(4,4)+X(4,1);A(1,2)=X(1,1)-X(1,4);A(2,2)=X(2,2)-X(2,3);A(3,2)=X(3,3)-X(3,2);A(4,2)=X(4,4)-X(4,1);A(1,3)=X(1,2)+X(1,3);A(2,3)=X(2,1)+X(2,4);A(3,3)=X(3,4)+X(3,1);A(4,3)=X(4,3)+X(4,2);A(1,4)=X(1,2)-X(1,3);A(2,4)=X(2,4)-X(2,1);A(3,4)=X(3,1)-X(3,4);A(4,4)=X(4,3)-X(4,2);%%%获取变换矩阵T=zeros(4);fori=0:3forj=0:3T(i+1,j+1)=cos(pi*(j+0.5)*i/4);endend%对A中的每一列做一维DCT,得到Bfori=1:4B(:,i)=T*A(:,i);endclearAT%%%由B元素的线性组合得到结果Y=zeros(4,4);Y(1,1)=(B(1,1)+B(1,3))*0.25;Y(2,1)=(B(2,1)+B(2,3))*sqrt(2)/4;Y(3,1)=(B(3,1)+B(3,3))*sqrt(2)/4;Y(4,1)=(B(4,1)+B(4,3))*sqrt(2)/4;Y(2,2)=0.5*(B(3,2)+B(3,4)+B(1,2))*0.5;Y(1,2)=(B(2,2)+B(4,4))*sqrt(2)/4;Y(4,4)=0.5*(B(1,2)-B(3,2)-B(3,4))*0.5;Y(3,4)=0.5*(B(2,2)-B(4,4)-B(4,2)-B(2,4))*0.5;Y(3,3)=0.5*(B(1,1)-B(1,3))*0.5;Y(2,3)=0.5*(B(2,1)-B(2,3)+B(4,1)-B(4,3))*0.5;Y(1,3)=(B(3,1)-B(3,3))*sqrt(2)/4;Y(4,3)=0.5*(B(2,1)+B(4,3)-B(2,3)-B(4,1))*0.5;Y(2,4)=0.5*(B(3,2)-B(1,4)-B(3,4))*0.5;Y(1,4)=(B(4,2)-B(2,4))*sqrt(2)/4;Y(4,2)=0.5*(B(1,4)+B(3,2)-B(3,4))*0.5;Y(3,2)=0.5*(B(2,2)-B(4,4)+B(4,2)+B(2,4))*0.5;end函数fidct2源码function[Y]=fidct2(X)%快速4X4DCT反变换%%By黄跃辉2016/11/30%%%检测输入量[a,b]=size(X);ifa~=4||b~=4Y=zeros(4,4);return;end%%%对X做归一化处理X(1,1)=X(1,1)*4;X(2,1)=X(2,1)*4/sqrt(2);X(3,1)=X(3,1)*4/sqrt(2);X(4,1)=X(4,1)*4/sqrt(2);X(2,2)=X(2,2)*2;X(1,2)=X(1,2)*4/sqrt(2);X(4,4)=X(4,4)*2;X(3,4)=X(3,4)*2;X(3,3)=X(3,3)*2;X(2,3)=X(2,3)*2;X(1,3)=X(1,3)*4/sqrt(2);X(4,3)=X(4,3)*2;X(2,4)=X(2,4)*2;X(1,4)=X(1,4)*4/sqrt(2);X(4,2)=X(4,2)*2;X(3,2)=X(3,2)*2;X=X./8;%%%对归一化X做线性组合B=zeros(4,4);B(1,1)=0.5*(X(1,1)+2*X(3,3));B(2,1)=X(2,1)+X(2,3)+X(4,3);B(3,1)=X(3,1)+X(1,3);B(4,1)=X(4,1)+X(2,3)-X(4,3);B(1,2)=X(2,2)+X(4,4);B(2,2)=X(1,2)+X(3,2)+X(3,4);B(3,2)=X(2,2)-X(4,4)+X(2,4)+X(4,2);B(4,2)=X(3,2)-X(3,4)+X(1,4);B(1,3)=0.5*(X(1,1)-2*X(3,3));B(2,