C语言实现DCT变换编码1

第.4:“Wherearithmeticprecisionisnotspecified,suchasthecalculationoftheIDCT,theprecisionshallbesufficientsothatsignificanterrorsdonotoccurinthefinalintegervalues.”逆DCT变换的过程这里不再详述，需要实现这个的可以去参考这个标准。在实际应用中一般通过两次1-DIDCT变换来完成2-DIDCT变换，这种方法通常被称为行－列法。一般来说，后者在结构上的对称性更好，并且可以重复使用硬件资源，所以在我们的芯片设计选用一种行-列法来进行IDCT单元的结构研究。二维IDCT可以分解成二次一维IDCT运算，如以下公式。在结构上，上两式所定义的运算使用了相同的运算“核”(如以下公式所示),它们具有相似性。因此利用三角函数的各种关系，可以得到“核”的快速算法。其中，为了便于理解，可将快速算法表示成蝶形图，如下图。将对模块进行一维IDCT变换的结果存储起来，转置输出，再进行一次IDCT变换，即为相应的二维IDCT变换。图4－8折叠结构的二维IDCT单元一行数据(一行有8个像素数据)在该单元中的处理流程是：１—>２—>３—>４—>５—>６—>７—>８。DCT变换探究1前言

此文适合于那些对DCT或对Haar小波的Mallat算法有一定了解的人。

由于我还是高一新丁，文学底子很薄弱，对于一些技术方面的知识，我是有口说不出，无法用文字表达出来，因此这里提供的知识只是我所知道的1/4左右，还有3/4我不知该如何表达，特别是第三节“深入研究DCT”，我个人认为简直是浅入！

如果你只是菜鸟，不但想看懂此文，而且还要看懂其他的类似文章，那么我教你一个最快的学习方法：

设X={10,20}

分解的方法：低频=10+20=30，高频=10-20=-10，

即Y={30,-10}

合并的方法：X(0)=(低频+高频)/2=(30+(-10))/2=10,X(1)=X(0)-高频=10-(-10)=20

即X={10,20}

只要搞清楚低频和高频是怎么来的和如何合并的即可。2DCT简介

DCT全名为DiscreteCosineTransform，中文名为离散余弦变换。在众人皆知的JPEG编码中，就是使用了DCT来压缩图像的。为什么DCT可以压缩图像？我想这个问题有很多人都想知道，但其实这是错误的说法！因为DCT在图像压缩中仅仅起到扶助的作用，给它n个数据，经变换后仍然会得出n个数据，DCT只不过消除了这n个数据的冗余性和相关性。

即，用很少的数据就能大致还原出这n个数据，其他的一些DCT系数只起到修正的作用，可有可无。

DCT有一个缺点，就是计算量很大！因为如果按照DCT的标准变换公式(二维)来实现8x8点阵的变换需要将近上万次计算！后来提出了一种优化方法，即将二维的DCT分解为两个一维的DCT，这样一来计算量就可以减少为原来的1/4。但是计算量依然巨大，不具有使用价值，后来在1988年有人提出了一种快速算法叫AAN,它也是将二维的DCT分解成一维的形式，但是二维计算量已减少到只有600来次了，JPG和MPEG编码中的DCT就是使用AAN算法实现的。

DCT还有一个缺点，就是不能无损变换，因为DCT系数都是一些无理数，目前为止，依然无法解决。3深入研究首先让我们来看看AAN算法的第一阶级变换代码：

ForI=0To3

J=7-I

Y(I)=X(I)+X(J)

Y(J)=X(I)-X(J)

NextI

设X={10,20,30,40,50,60,70,80}

那么Y={90,90,90,90,-10,-30,-50,-70}

可以看出，这一阶级的低频部分(相加得出的数据)全部相等，而高频部分则呈线性或者是有规律的。DCT之所以能以较少的数据大致还原图像，就是因为通过预测高频部分而达到的。那么为何高频部分可以预测呢？请仔细看上面的代码，可以看出DCT是由外到内来进行处理的，由于像素及像素间有一定的关联性，所以靠的越近的像素之间的差就应该越小，越远就因该越大，但也并不是说所有的数据都具有这种规律，因此DCT预测出来的高频数据就会和原高频数据不大相同，它们之间的差便是第二节提出的修正数据。第二阶级变换则是在第一阶级变换的基础上再次分解出低、高频，和预测高频，得出修正值。第三阶级……。最后，再将DCT系数按照重要程度由大到小，由左到右，重排列即可。

例：X={10,20,30,40,50,60,70,80}

经过FDCT后：Y={127,-64,0,-7,-0,-2,0,-1}

其中127是最最重要的，而-64次之，以此类推。可以发现,-7,-2,-1的能量都很小，说明这三个修正值可以忽略，当忽略后，

得Y={127,-64,0,0,0,0,0,0}

经过IDCT后：X={14,18,27,39,51,63,72,76}

这及原始数据：X={10,20,30,40,50,60,70,80}

是非常接近的，肉眼很难发觉。4为何JPEG2000放弃DCT

在JPEG2000里，放弃了基于块的DCT，而改为了小波变换。为何要放弃DCT呢？我认为最根本的原因还是跟DCT的计算量有关，第二节已经指出，为了减少计算量，我们不得不使用只能处理8X8点阵的AAN快速算法（目前，也只有基于8X8点阵的），对于一幅图像，必须将其分割成无数个8X8大小的“块”，对块进行变换。在低码率下，就会产生方块效应，要解决这个问题，唯有不使用基于区块的AAN快速算法，而是使用直接变换法，但计算量惊人！由于小波变换计算量很少，便于直接处理图像数据，因此就不会产生块效应，但假如用小波也进行基于8X8点阵的块变换，在低码率下，同样也会有块效应！只要是基于块变换的，那么在低码率下就会出现块效应，无论是DCT还是小波。因此，如果忽略DCT直接处理的计算量问题的话，我认为压缩效率会比JPEG2000更好！（具体原因暂不讨论）5DCT的改进

下面的代码是我对DCT变换的改进，它具有以下特性

l无损变换

l计算量少

l原位计算

经改进后，它已不再叫作DCT了，可以认为是一种新的算法，只不过是在DCT的基础上修改而来。

以下是正变换：X为输入端，Y为输出端

设X={10,20,30,40,50,60,70,80}那么Y={45,40,0,0,0,0,0,0}

'第一阶级

ForI=0To3

J=7-I

Y(J)=X(J)-X(I)

Y(I)=X(I)+Fix(Y(J)/2)

NextI

'第二阶级

ForH=0To4Step4

ForI=0To1

J=3-I

X(J+H)=Y(J+H)-Y(I+H)

X(I+H)=Y(I+H)+Fix(X(J+H)/2)

NextI

NextH

'第三阶级

ForI=0To6Step2

Y(I+1)=X(I+1)-X(I)

Y(I)=X(I)+Fix(Y(I+1)/2)

NextI

'预测

Y(3)=Y(3)-Y(2)

Y(6)=Y(6)-Y(7)

Y(7)=Y(7)-Y(4)

'重要性排序及AAN一样，皆为{0,4,2,6,1,5,7,3}，此略为何能无损？为何能原位？和具体实现原理暂时略，以后我会补上6参考

[1]丁贵广，计文平，郭宝龙VisualC++6.0数字图像编码p44,p57,p170快速DCT变换仿效FFT的FDCT方法有及DCT无关的复数运算部分，选用代数分解法可以降低运算量，达到高速运算的目的。代数分解法实现如下：对一维DCT表达式直接展开，寻找各点表达式中共同项，仿FFT蝶形关系，将表达式中的共同项作为下一级节点，依次进行多次，最后得到变换结果。一、DCT部分例子：Definecos(n*pi/16)Cn

F(0，v)=0.5*C(0)*[x(0)+x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)]F(1，v)=0.5*C(0)*[x(0)*C1+x(1)*C3+x(2)*C5+x(3)*C7+x(4)*C9

+x(5)*C11+x(6)*C13+x(7)*C15]

=0.5*{[x(0)-X(7)]C1+[X(1)-X(6)]*C3+[X(2)－x(5)]*C5

+[x(3)-x(4)]*C7]从上面的式子可以看到07,16,25,34可以作为第一次运算的相加节点，将所有节点的表达式列出后，可发现一个规律，得到一蝶形图，按之编程，如下：#define

C10.9808

#define

C20.9239

#define

C30.8315

#define

C40.7071

#define

C50.5556

#define

C60.3827

#define

C70.1951//先做行DCT

voidfdctrow(double*blk)

{

doubleS07,S16,S25,S34,S0734,S1625;

doubleD07,D16,D25,D34,D0734,D1625;S07=blk[0]+blk[7];

S16=blk[1]+blk[6];

S25=blk[2]+blk[5];

S34=blk[3]+blk[4];

S0734=S07+S34;

S1625=S16+S25;D07=blk[0]-blk[7];

D16=blk[1]-blk[6];

D25=blk[2]-blk[5];

D34=blk[3]-blk[4];

D0734=S07-S34;

D1625=S16-S25;blk[0]=0.5*(C4*(S0734+S1625));

blk[1]=0.5*(C1*D07+C3*D16+C5*D25+C7*D34);

blk[2]=0.5*(C2*D0734+C6*D1625);

blk[3]=0.5*(C3*D07-C7*D16-C1*D25-C5*D34);

blk[4]=0.5*(C4*(S0734-S1625));

blk[5]=0.5*(C5*D07-C1*D16+C7*D25+C3*D34);

blk[6]=0.5*(C6*D0734-C2*D1625);

blk[7]=0.5*(C7*D07-C5*D16+C3*D25-C1*D34);//再做列DCTvoidfdctcol(double*blk)

{

doubleS07,S16,S25,S34,S0734,S1625;

doubleD07,D16,D25,D34,D0734,D1625;S07=blk[0*8]+blk[7*8];

S16=blk[1*8]+blk[6*8];

S25=blk[2*8]+blk[5*8];

S34=blk[3*8]+blk[4*8];

S0734=S07+S34;

S1625=S16+S25;D07=blk[0*8]-blk[7*8];

D16=blk[1*8]-blk[6*8];

D25=blk[2*8]-blk[5*8];

D34=blk[3*8]-blk[4*8];

D0734=S07-S34;

D1625=S16-S25;blk[0*8]=0.5*(C4*(S0734+S1625));

blk[1*8]=0.5*(C1*D07+C3*D16+C5*D25+C7*D34);

blk[2*8]=0.5*(C2*D0734+C6*D1625);

blk[3*8]=0.5*(C3*D07-C7*D16-C1*D25-C5*D34);

blk[4*8]=0.5*(C4*(S0734-S1625));

blk[5*8]=0.5*(C5*D07-C1*D16+C7*D25+C3*D34);

blk[6*8]=0.5*(C6*D0734-C2*D1625);

blk[7*8]=0.5*(C7*D07-C5*D16+C3*D25-C1*D34);voidfdct(double*block)

{

inti;

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

fdctrow(block+8*i);

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

fdctcol(block+i);

}二、IDCT部分

图片来源：G.G.Pechanek,C.W.Kurak,C.J.Glossner,C.H.L.Moller,andS.J.WalshIBMMicroelectronicsDivision,ResearchTrianglePark,N.C."M.F.A.S.T.:AHIGHLYPARALLELSINGLECHIPDSPWITHA2DIDCTEXAMPLE"图中未给出系数，需自行算出，仿上述DCT方法直接展开表达式搜寻规律即可。编程如下：#define

C10.9808

#define

C20.9239

#define

C30.8315

#define

C40.7071

#define

C50.5556

#define

C60.3827

#define

C70.1951//对行做DCTvoididctrow(double*blk)

{

doubletmp[16];

//firststep

tmp[0]=blk[0]*C4+blk[2]*C2;

tmp[1]=blk[4]*C4+blk[6]*C6;

tmp[2]=blk[0]*C4+blk[2]*C6;

tmp[3]=-blk[4]*C4-blk[6]*C2;

tmp[4]=blk[0]*C4-blk[2]*C6;

tmp[5]=-blk[4]*C4+blk[6]*C2;

tmp[6]=blk[0]*C4-blk[2]*C2;

tmp[7]=blk[4]*C4-blk[6]*C6;

tmp[8]=blk[1]*C7-blk[3]*C5;

tmp[9]=blk[5]*C3-blk[7]*C1;

tmp[10]=blk[1]*C5-blk[3]*C1;

tmp[11]=blk[5]*C7+blk[7]*C3;

tmp[12]=blk[1]*C3-blk[3]*C7;

tmp[13]=-blk[5]*C1-blk[7]*C5;

tmp[14]=blk[1]*C1+blk[3]*C3;

tmp[15]=blk[5]*C5+blk[7]*C7;

//secondstep

tmp[0]=0.5*(tmp[0]+tmp[1]);

tmp[1]=0.5*(tmp[2]+tmp[3]);

tmp[2]=0.5*(tmp[4]+tmp[5]);

tmp[3]=0.5*(tmp[6]+tmp[7]);

tmp[4]=0.5*(tmp[8]+tmp[9]);

tmp[5]=0.5*(tmp[10]+tmp[11]);

tmp[6]=0.5*(tmp[12]+tmp[13]);

tmp[7]=0.5*(tmp[14]+tmp[15]);

//thirdstep

blk[0]=tmp[0]+tmp[7];

blk[1]=tmp[1]+tmp[6];

blk[2]=tmp[2]+tmp[5];

blk[3]=tmp[3]+tmp[4];

blk[4]=tmp[3]-tmp[4];

blk[5]=tmp[2]-tmp[5];

blk[6]=tmp[1]-tmp[6];

blk[7]=tmp[0]-tmp[7];

/*

blk[0]=0.5*(Y0C4+Y2C2+Y4C4+Y6C6+Y1C1+Y3C3+Y5C5+Y7C7);

blk[1]=0.5*(Y0C4+Y2C6-Y4C4-Y6C2+Y1C3-Y3C7-Y5C1-Y7C5);

blk[2]=0.5*(Y0C4-Y2C6-Y4C4+Y6C2+Y1C5-Y3C1+Y5C7+Y7C3);

blk[3]=0.5*(Y0C4-Y2C2+Y4C4-Y6C6+Y1C7-Y3C5+Y5C3-Y7C1);

blk[4]=0.5*(Y0C4-Y2C2+Y4C4-Y6C6-Y1C7+Y3C5-Y5C3+Y7C1);

blk[5]=0.5*(Y0C4-Y2C6-Y4C4+Y6C2-Y1C5+Y3C1-Y5C7-Y7C3);

blk[6]=0.5*(Y0C4+Y2C6-Y4C4-Y6C2-Y1C3+Y3C7+Y5C1+Y7C5);

blk[7]=0.5*(Y0C4+Y2C2+Y4C4+Y6C6-Y1C1-Y3C3-Y5C5-Y7C7);

*///在对列做DCTvoididctcol(double*blk)

{

doubletmp[16];

//firststep

tmp[0]=blk[0*8]*C4+blk[2*8]*C2;

tmp[1]=blk[4*8]*C4+blk[6*8]*C6;

tmp[2]=blk[0*8]*C4+blk[2*8]*C6;

tmp[3]=-blk[4*8]*C4-blk[6*8]*C2;

tmp[4]=blk[0*8]*C4-blk[2*8]*C6;

tmp[5]=-blk[4*8]*C4+blk[6*8]*C2;

tmp[6]=blk[0*8]*C4-blk[2*8]*C2;

tmp[7]=blk[4*8]*C4-blk[6*8]*C6;

tmp[8]=blk[1*8]*C7-blk[3*8]*C5;

tmp[9]=blk[5*8]*C3-blk[7*8]*C1;

tmp[10]=blk[1*8]*C5-blk[3*8]*C1;

tmp[11]=blk[5*8]*C7+blk[7*8]*C3;

tmp[12]=blk[1*8]*C3-blk[3*8]*C7;

tmp[13]=-blk[5*8]*C1-blk[7*8]*C5;

tmp[14]=blk[1*8]*C1+blk[3*8]*C3;

tmp[15]=blk[5*8]*C5+blk[7*8]*C7;

//secondstep

tmp[0]=0.5*(tmp[0]+tmp[1]);

tmp[1]=0.5*(tmp[2]+tmp[3]);

tmp[2]=0.5*(tmp[4]+tmp[5]);

tmp[3]=0.5*(tmp[6]+tmp[7]);

tmp[4]=0.5*(tmp[8]+tmp[9]);

tmp[5]=0.5*(tmp[10]+tmp[11]);

tmp[6]=