二维DCT的优化算法

,《音视频编码技术》课程报告《2-DDCT算法优化》姓名：学号：院系：电子与信息工程学院指导教师：摘要离散余弦变换（DCT）是一种正交变换，它本身并不会改变信息源的燧值，因而是一种无损的变换。但是由于•离散余弦变换（DCT）变换能够给量化和编码等过程创造很好的条件，因而在数字信号处理领域，特别是语音和图像压缩领域有着很广泛的应用。对于语音信号与图像信号这样相关性很强的信号而言，离散余弦变换（DCT）接近于最优的KL变换（Karhunen-LoeveTransform）＜＞传统的2-D离散余弦变换（DCT）变换过程是行-列变换，这种方法直观简单但是效率低3在仃些实时性要求较高的图像压缩领域并不能满足要求。从离散余弦变换（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的变换矩阵是固定的，与图像无关。而且它是构造成对称的数据序列，避免了子图像边界处的跳跃和不连续现象。傅里叶变换是应用最早的变换之一，并且也有快速算法。但是傅里叶变换存在一个很明显的缺陷，恢复出来的边界会存在不连续现象，这使得恢复的图像存图1离散余弦变换（DCT）效果图在明显可见的方格，影响图像质量。沃尔什变换比DCT简单，实时性更高，但是性能比DCT略差。基于以上的特点，DCT具有广泛的应用，在常见的JPEG静态图像编码以及MJPEG、MPEG动态图像编码等标准中都有应用。这种应用的广泛也对DCT的算法效率提出了更高的要求，对DCT快速算法也有着更加迫切的需求。（二）前人研究由于DCT在语音和图像压缩领域的广泛应用，已经有人提出了许多快速算法和快速计算DCT的VLSI架构。DCT的传统计算方法是行-列法，先对行计算1-DDCT,再对列计完1-DDCTo对于NXN大小的图像而言，这种算法需要计算2N个1-DDCT,算法的复杂度为0（N2）。因而为了更加高效地进行DCT运算，有人提出了直接对整个2维矩阵做DCT运算。目前，最为高效的DCT算法是由Duhamel和Guillemot在1990年的一篇论文中提出的多项式变换计算方法。在这种方法中，2-DDCT运算的乘积运算相对于传统的方法而言减小了50%。还有人提出基于FFT的算法，这种算法的复杂度在传统方法的50%到75%之间。此外，Cho和Lee提出了专门针对于矩阵大小为4X4快速DCT算法，这种算法仅仅局限于大小为4X4的变换，如果要进行扩展更大的维数，则需要与迭代算法相结合，但是如果维数很大，这种迭代算法这会大大增加整体尊法的复杂度。在本次的课程报告中所展示的是一种有别于以上快速算法的方法，这种方法和Ch。与Lee提出的算法有异曲同工之处，本质上来说是对他们算法的一种扩展。这种算法对于大小为NXN的矩阵而言，仅仅需要计算N个l-DDCT,但是要求N=2，mCR。算法需要的乘积运算次数仅仅为传统算法的50%,与Duhamel和Guillemot的算法有着相同复杂度。但是相比于Duhamel和Guillemot的算法而言，这种算法的结构有着更高的对称性，更加简单直观，同时能够更加容易的应用到硬件的并行处理中，也更加适合VLSI的实现。(=)报告的章节安排本次报告的主要内容和章节安排如下：.快速DCT算法的理论分析从数学上对算法进行理论分析，相关公式推导，总结出最终的算法结构。.相应的快速IDCT算法的理论分析基于前面的理论推导，总结出快速IDCT算法的结构。.与其他快速算法的复杂度比较将此算法与传统的算法、其他算法的复杂度作数值上的比较。.算法的实现与测试在前面理论分析的基础上，在MATLAB平台上实现4X4的快速算法，井对图像做DCT和IDCT变换，检验算法的有效性，同时与MATLAB自带的DCT和IDCT函数的运算时间进行比较。.总结与展望分析此算法的优缺点与可改进的空间。

二、理论分析(―)2-DDCT理论分析对于一个给定的大小为NxN的二维矩阵X殉£X,i,j=0,1,2,•••—1假设二维DCT矩阵为yymnW匕m,n=0,1,2,…，N—1DCT变换公式为：K加=c(m)c(7i)Z篇7才为8sa^^l)c。生(2.1)其中，jm=0J^，n=0c(m)=彳厂c(n)=,厂舟…[>n*°定义将n为4"对c(m)c(tz)做归一化得到的变量ymn=2M包cos；烂萨)c。罐嗜巧(22a)Ymn=c(rn)c(n)ymn(2.2.b)利用枳化和差公式，可以进行如卜的变换，/(2i+l)nwr\/(2;+l)n7r\1f/(2i+l)m+(2;-l)n/(2i+l)m-(2;+l)n'、/38s8S(号丁)=5(8S(7T)+8S(汗))(2.3)将(2.3)带入到(22a)中，得到2121/.、\l=UJ=UV，(2i+l)m-⑵+l)n+22/八叫2N\l=UJ=UV，(2i+l)m-⑵+l)n+22/八叫2N1i=0J=0'formfn=L；同时，我们定义两个新的变量4rm,8恤4m=SfLo1xtjcosr)2,…，N-l(2.4)Bnn=X..灯勺8S.+D：；⑵+%)Q5b)所将71可以表示为y)nn=5(4mn+^mn)(2.6)从(2.5)与(2.6)中，可以发现，二维的DCT变换可以分解为两个分量的叠加，并且这两个分量(4皿与为皿)的形式类似于一维DCT的形式。那么下面的目标就在于如何将4ml与5M转换成一维DCT的形式。从而达到对于NXN的矩阵而言，只需要计算N次一维DCT就可以获得二维DCT。一维DCT的表达式为NT/.、T.(2i+l)mnYm=c(m)工xtjcosE^-)与(2.5.a)和(2.5.b)对比可以发现，如果｛(2i+l)m±(2j+1)叫可以转换为｛(2i+l)k,kGR｝的形式，那么就可以将4ml与B7nn转换成一维DCT的形式。由于｛ij=0,12,…,N-1｝,那么｛⑵+1)｝的值要么等于｛(2i+1)｝整数倍除以2N的余数，要么等于2N减去(2i+1)的整数倍除以2N的余数。(274-1)=p(2i+1)mod2N(2.7.a)或者(2;+1)=2N—p(2i+1)mod2N(2.7.a)其中p=L3,5,7,...，N-1,如果p超过这个范围，这会得到相同的(2/+1)的值。将(2.7)进行化简，就可以得到j=(pi+mod2N(2.8.a)或者j=(N-1-pi-mod2N(2.8.a)其中p=l,3,5,7,…,N-l。从(2.8)中可以发现，如果i取0至U(N—1)的不同值，那么对应得至IJ的j的序列是不一样的，这样就得到了N2个不同的(ij)组合，这说明，依据(2.8)的公式可以遍历NXN矩阵内的所有元素。此外每一个i的值都对应着一个j的序列，这意味着，依据(2.8)可以对输入数据分成N组，每一个组内元素的下标&力都满足(2.8),每一个分组表示为/(p|a)=(pi+、一)mod2N(2.9.a)或者j(p\b)=(N-1—pi--■)mod2N(2.9.a)其中，p=1,3,57…,N—l,i=0,1,2,-,/V-l二维输入数据口m"=0,12...,/7-1｝通过(2.9)就可以分成N个一维数据i=p=1,3,5,-,N-1)或者{4/(p|a)，i=1,2,…，N—l,p=1,3,5,…，N—1)用用和畤来表示每一组数据Rp={苞j(p⑷,i=1,2,…，N-lj(p|a)=(pi+^—)mod2Mp=135,7,…，N-l}畤={々j(p|b)，i=1,2，…，N—lj(p|b)=(N-1-pi-mod2N,p=135,7「・・，N—1}(210)为了进一步进行简化，假设己知(pi+?)除以N的商为qp「于是(2.10)可以化简为P—1Rp={%j(p|a)，i=12…，N—lj(p|a)=pt+———Nqphp乙=1,357,…,N-l}.p—1R；==L2,…,N-lj(plb)=N-1-pi--+Nqpimod2Mp=1,357,…，N-1}(MO)按照(2.10)中的分组，4刖与巴的可以表达为，4M=X；T(V(m,〃)+若(犯")}(2.11.a)Bmn=Ep15,p(jn,n)+(m,n)(2.11.b)其中邛(叫几)=X”叫/8显2，+*；⑵+1)[)q.12.3)

Tp(m,n)=eRb%i；co显2f2/+i”e)(2.12.b)S£(m,n)=Xfw堞与co."叱2…n)(2.12©S*m,n)=%的.85r+吗；(2尸1川r(2.12.d)把(2.11)带入到(2.6)中，可以得到，yrnn=ErT(^(m,n)+*(m,7i)+S£(m,n)+S?(mm)}(2.13)为了把必加表示为一维DCT的和的形式，需要将Tp(rnfn).Tp(m,n),S'(m,n)、S夕(m,n)转换成一维DCT的形式。把(2.10)带入到(2.12.a)中,可以得到,九)2N校,巾£殉…舒2-P⑵九)2N将上面的式子拆分为n为奇数和偶数两种情况:当n为奇数时N-lZi%(p|a)CN-lZi%(p|a)C。迎i=0(2i+l)m+(m+pn)2N几)(2.14.a)当n为偶数是Tp(m,n)=当n为偶数是Tp(m,n)=N-i…、Vr1、qq⑵+1)加+⑺+pn)/.JI)%与J(p|a)CO武-TT)(2.14.b)那么以同样的方式，可以获得4(犯八)、S;(m.n)、S?()n,7i)的表达式,松n,n)=,fNT‘V$2i+l)m+(m-pn)2松n,n)=,fNT‘V$2i+l)m+(m-pn)2句91犷。5甲i=0N-l一£(一i)qp』(pib)8si=02Nn),niseven(2.15.a)(21+l)m+(m—pn)2NnLnisodd(2.15.b)Sf(m,n)=vh(2i+Sf(m,n)=vh(2i+l)m+(m-p72)./xtj(pia)cos^-n)fniseveni=0IN-l•2(-1A旬3a)8迎⑵+l)L(m阿九)'nisoddJo(2.16.a)(2.16.b)S；g,n)=Z_lxij(p\b)cos^i=0N-l(2i+1)m+(m+pn)2Nit),niseven(2A7.d)把(2.14)-(2.17)NT工Pymn=W（-l）qp%3b）cS；g,n)=Z_lxij(p\b)cos^i=0N-l(2i+1)m+(m+pn)2Nit),niseven(2A7.d)把(2.14)-(2.17)NT工Pymn=W（-l）qp%3b）c。典(21+1)m+(m+pn)2Ntt),nisodd(2」7.b)7=0带入到（2.13）中，可以得到%刖的新表达式,.、(2i+l)(m+np)/a(pla)(pli))cos2jV"+jV-1Vr,\(2i+l)(m-np)/(xt/(p|a)+x(/(p|b)J8S—117=0,nIseven(2.18.a)/Vr…\(2i+l)(m+np)\_2(-1)乂为31。)一沏(p|b))8s-n+“4Vr(\(2i+l)(m-np)2/)(和*373㈤)8s2N7r/jV-1,nisodd(2.18.b)可以发现,N-lW(4(p|a)+xtj(p|b))cosi=0N-l2(xi；(p|a)+xij(p|b))cosN-lW(4(p|a)+xtj(p|b))cosi=0N-l2(xi；(p|a)+xij(p|b))cosi=0(2i+l)(m+np)2Nzr(2.19.a)(2i+l)(m—np)2N7r(2.19.b)N-lW(—(殉(p|a)—殉(p|b))1=0NT/.(一1加(/(p|a)一殉(p|b))i=0coscos(21+l)(m+np)2Nn(2.19.c)(21+l)(m—np)2Nn(2.19.d)（2.19.a）和（2.19.b）分别是｛%⑦⑷+叼⑺⑶｝序列的一维DCT变换序列的第（m+叩）和（m-np）个元素；（2.19.€：）和（2.19.£1）分别是｛（一1尸卬（狗（肉）一罚9田）｝序列的一维DCT变换序列的第（m+他）和（m-zip）个元素。定义两个量:fpi和外］(2i+1)1

(2i+1)1

2N-7r(2.20)fpl=W(Xi/(P|a)+X4(p|b))cosN-lf、(2i+1)/9pi=L(一1)①(殉(pla)+殉(p|b)Jcos—忝—nQ21)i=0那么(2.19.3)和(2.19.b)就分别对应着土加,1=0,1,2,…，N-l：(2.19.C)和（Z19.d）就分别对应着土gw」=0,l,2,・>N—l0这就意味着，大小为NxN矩

阵的二维DCT只需要先计算出N个一维DCTo从（2.18）可以看出，要计算出K所，可以先计算出加和。加，然后再将相对应的加和外1进行线性叠加，得到最终结果。已知，n为偶数时,/(2i+l)(m-?i(N-p))、/(2i+l)g+叩)\cos(n)=±cos(町(2.22)+xij(p\b))COS⑵+l)(m+xij(p\b))COS⑵+l)(m+np)2N那么把（2.22）带入就可以得到,NT士|N-p)1=/NT士|N-p)1=/&/(N-p|a)+啊(N-p|b))8s(2i+l)(m+n(n—p))2N7T(2.23)同样的方式，当n为奇数时，同样的方式，当n为奇数时，NTg“=W(T)%(殉(p|a)+xij(plb))cost=0N-l土g(N-p)(NT)=，(-1片(勺3-p\a)i=0(2i+l)(m+np)2N+xij(N-plb))C0S+xij(N-plb))C0S(2i+l)(m+(N-p)n)2N(2.24)同样，gpi会以上面的两个式子意味着，加总是会以±/^_加的形式出现。同样，gpi会以土g（N-p）（NT）的形式出现。基于以上的分析，假如输入矩阵的大小为8X8,那么这个算法就可以用以下的图来表示，

1-DDCT1-DDCTT171-DDCT1-DDCT7/Z1-DDCT1-DDCT—fao-fsi1-DDCT1-DDCTT171-DDCT1-DDCT7/Z1-DDCT1-DDCT—fao-fsi—(53—f54——f»6-工—fel建r—fee-tte•*-f7O一f71——tn一(73二篇一m1-DDCT—6n图2（续）8X8DCT算法结构图，。表示相加，虚线表示取反图28X8DCT算法结构图，0表示相加，虚线表示取反（二）2-DIDCT理论分析对于正交变换而言，如果不考虑系数的影响，那么其相对应的反变换就直接是正变换的反过程而己。那么可以利用这个来推导出快速IDCT算法。如果在DCT变换中不考虑（2.1）中前面的系数，那么反变换的算法结构就仅仅是正变换算法结构的取反。例如对于8X8的IDCT变换而言，如果不考虑系数，那么其算法结构如图3所示。可以看到，相对于图2而言，图3仅仅是将图2中信号的流向做了取反，相当于是反向做运算，然后做N个一维IDCT,其运算鼠与（一）中分析的DCT的运算鼠相当，相对于传统的IDCT算法而言也仅仅需要50%的一维DCT变换次数。图38X8IDCT算法结构图,。表示相加，虚线表示取反fio川抬抬118盘耨爵盘能名fio川抬抬118盘耨爵盘能名433言名rg&g轲012O4687I2G466ro1734be734s4or6CS5SK5-&M6A(s&&GfiAJQNT777777^V"r图4（续）8X8IDCT算法结构图，。表示相加，虚线表示取反但是如果将系数考虑进来，就需要做一些额外的调整，调整的过程如下：1）对于初始输入数据对相对应的DCT系数做归一化处理2）将第一步得到的数据乘以系数JO3）每个一维IDCT的第一个输入数据需要再乘以系数：经过以上三个步骤，处理后就可以得到相对应的IDCT数据。（三）算法的复杂度分析从（一）和（二）部分的内容中可以看到对于NXN为的DCT变换，仅仅需要计算N个一维的DCT。我们在本小节中以需要的乘法运算的次数和加法运算次数来作为算法复杂度比较的标准。由于NXN大小的二维DCT变换仅仅需要N次一维的DCT变换，那么相对于传统的二维DCT算法而言，本报告中的算法所需要的乘法运算量就仅仅是传统方法的50%。而对于一维DCT的算法复杂度，目前比较高效的算法需要的乘法运算量为N510g2N如果在二维DCT变换中使用这个一维DCT算法，那么所需乘法运算量为，N2■—log2N(2.25)表格1所需乘法运算次数比较矩阵大小传统算法[1][2][3]本报告算法4X432241616168X8192144104969616X16102476856851251232X3251203840274025602560表格2所需加法运算次数比较矩阵大小传统算法[1][2][3]本报告算法4X472727068748X846446446248446616X162592259225582531253032X321337613376129501257812738至于所需的加法运算量，从前面的分析中可以知道，除了DCT运算的加法运算量，其他地方加法运算量为N2(l+log2N)-2N-(N-2)。而对于一维DCT变换而言，所需的加法运算量为号log2N-N+l0所以二维DCT算法所需的加法运算量为，5N2—log2/V-2A+2(2.26)依据(2.25)和(2.26)和其他参考文献中的快速算法，列出了表格1和表格2,分别展示了各个算法所需的乘法运算次数和加法运算次数。从表格中可以发现，本报告中的算法所需的乘法运算量与文献⑶中一样，并且比其他算法都要低。同时，加法运算量与其他算法基本一致。三、实现与测试在本章节中，将根据第二章节中的DCT和IDCT快速算法的理论分析，在MATLAB平台上实现相应的DCT变换算法和IDCT变换算法。然后利用这两个算法对一个实际图像做变换与反变换处理，观察处理的效果。(一)算法的实现在本小节里面，将会依据第二章节的分析，在MATLAB平台上实现4X4的快速DCT与IDCT算法。从第二章中的分析中可以知道，快速算法可以用流向图很简洁地表示出来。那么依据前面的公式，对于4X4的DCT和IDCT变化可以得到以下的流向图，图54X4算法结构图，(a)DCT(b)IDCT那么依据图5,算法的流程图如下:根据(210)分组图6快速DCT算法流程图

图7快速IDCT算法流程图（二）算法的测试测试环境：Dell商用台式机，MATLAB2013a,测试图片大小为614x614.朦图DCT变帙亚图朦图DCT变帙亚图图8快速DCI与IDCT测试效果图从图8的实际测试结果来看，原图经过了4X4的DCT变换之后，在每个4X4大小的方格内，图像的能量都集中分布于少数几个系数上，这符合正交变换的特点。同时经过了IDCT变换之后，可以完整的恢复图像原来的模样，且不存在失真与格纹现象。这证明了算法的正确性。表格3程序运行时间比较自定义函数MATLAB已有函数函数名称调用次数运行时间函数名称调用次数运行时间main13.969smain16.804sfdct2237161.594sdct2237163.228sfidct2237161.396sidct2237162.674s表格3是MATLAB利用CPU的运行时间计算出来的结果，虽然结果的具体数值并不具有很高的可靠性，但是在相同条件下得出来的结果，相互之间具有可比性。通过比较可以发现，在相同的调用次数的条件下，无论是fdct2函数还是fidct2函数，其运行时间都比MATLAB己有函数的运行时间减少了将近50%o!1!总结与展望!1!总结与展望本次课程报告首先分析了离散余弦变换（DCT）在图像压缩领域的重要性，它本身虽然不具有压缩能力，但是变换给后续的编码和量化操作带来了便利，在多个图像、视频标准中都有应用。DCT的广泛应用也为人们对它的研究增添了动力，由于传统计算方法的效率低下，许多文献提出了相关的DCT快速算法。本次报告就研究、实现和分析了其中一种性能比较优越的快速算法。这种快速算法的乘法运算量仅仅只有传统算法的50%,而加法运算量和其他方法相当，具有很高的效率，同时，这种算法可以能够用流向图来直观的表示，简单而直观。而且由于这种算法具有比较好的对称性，它能够很容易地应用到多线程并行处理中，这样效率乂可以得到很大的提升，1司时这种算法还可以直接用VLSI结构来实现，这说明这种算的应用相对简单。在进行了理论分析之后，基于理论分析在MATLAB平台上实现了此算法，并用实际图片进行了测拭，测试结果表明，此算法可以快速地进行DCT变换和IDCT变换，能够完整的恢复出原图。此外，还将此算法与MATLAB自带的DCT变换和IDCT变换函数在相同的条件下进行了运行时间比较，比较结果表明，此算法的运行时间相对于MATLABH带的函数减少了将近50%,具有较高的效率。虽然这种的算法具有很高的效率，但是它依然有它的局限性。这种算法的最大的局限性在于，如果二维矩阵的维数很大，那么这种算法将会难以实现。因为这种算法需要对输入数据进行排列，在进行了一维DCT后还需要进行多次的线性运算，这种排列与线性运算的组合并不具有很明显的规律性。在算法的实现过程中，对于16X16大小的矩阵，需要列出256个式子：但是对于32X32的矩阵，需要的式子达到了1024个。所以对于大型的矩阵，这种方法是不适合的，比较适合于16X16以下的矩阵。五、参考文献C.Ma,“AfastrecursivetwodimensioncosinetransformIntelligentRobotsandComputerVisio:SeventhinaSeries,DavidP.Casasent,Ed.,inProc.SPIE1002,pp.541-548,1988.M.Vetterli,“Fast2-Ddiscretecosinetransform/inProc.ICASSP'85,Mar1985.RDuhamelandC.Guillemot,"'Polynomialtransformcomputationof2-DDC!/'inProc.ICASSP790,pp.1515-1518,Apr1990.N.I.ChoandS.U.Lee,“Afast4X4DCTalgorithmfortherecursive2-DDC]"IEEETrans.Acoust,Speech,SignalProcessing,submittedforpublication.RA.KamangarandK.R.Rao,"Fastalgorithmforthe2-Ddiscretecosinetansform/IEEETrans,Comput,vol.C-31,pp.899-906,Spet.1982.B.G.Lee,“Anewalgorithmtocomputethediscretecosinetransform:IEEETrans.Acoust.,Speech,SignalProcessing,vol.ASSP-32,pp.1243-1245,Dec.1987.[7]杨帆.数字图像处理与分析.北京：北京航空航天大学出版社,2007⑻全红艳，曹桂涛.数字图像处理原理与实现方法.北京：机械工业出版社，2014六、附件（二）函数（二）函数fidcQ源码(一)函数fdct2源码function[Y]=fdct2(X)%快速4X4DCT变换函数%%By黄跃辉2016/11/30%%%检测输入量[a,b]=size(X);ifaM=411b4Y=zeros(4,4);return;end%%A=zeros(4,4);B=zeros(4,4);A(L1)=X(1,1)+X(L4);A(2,l)=X(2/2)+X(23)；A(X1)=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(23)；A(3,2)=X(3,3)-X(3,2);A(4,2)=X(4,4)-X(4,l);A(13)=X(l,2)+X(13)；A(23)=X(24)+X(44)；A(3,3)=X(3,4)+X(3,1);A(43)=X(43)+X(4,2)；A(l,4)=X(l/2)-X(13)；A(2,4)=X(2/4)-X(2,l);A(3,4)=X(3,1)-X(3,4);A(4,4)=X(43)-X(4,2)；%%%获取变换矩阵T=zeros(4);fori=0:3forj=03T(i+l,j+l)=cos(pi*(j+0.5)*i/4);endend%对A中的每一列做一维DCT,得到Bfori=l:4B(:,i)=T*A(:J);endclearAT%%%由B元素的线性组合得到结果Y=zeros(4,4);Y(l,l)=(B(l,l)+B(13))*0.25;Y(2,l)=(B(2,l)+B(2,3))*sqrt(2)/4;Y(3,l)=(B(Xl)+B(33))*sqrt(2)/4;Y(44)=(B(4,l)+B(43))*sqrt(2)/4;TOC\o"1-5"\h\zY(2,2)=0.5*(B(3,2)+B(3,4)+B(l,2))*0.5;Y(l,2)=(B(2,2)+B(4,4))*sqrt(2)/4;Y(4,4)=0.5*(B(L2)-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(33)=0.5*(B(l/l)-B(13))*0.5;Y(23)=0.5*(B(2/l)-B(23)+B(4/l)-B(43))*0.5;Y(13)=(B(3/l)-B(33))*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(l,4)-B(3,4))*0.5;Y(l/4)=(B(4,2)-B(2,4))*sqrt(2)/4;Y(4,2)=0.5”(B(l,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;endfunction[Y]=fidct2(X)function[Y]=fidct2(X)%快速4X4DCT反变换%%By黄跃辉2016/11/30%%%检测输入最[a,b]=size(X);ifa**=411b"=4Y=zeros(4,4);return;end%%做归一化处理X(U)=X(14)*4;X(24)=X(24)*4/sqrt(2);X(3,l)=X(3,l)*4/sqrt(2);X(4,l)=X(44)*4/sqrt(2);X(2,2)=X(2,2)*2;X(l,2)=X(l,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(13)=X(13)*4/sqrt(2);X(43)=X(4