C语言实现DCT变换编码.doc

Contents1. DCT变换编码C语言2. MPEG4中逆DCT变换3. DCT变换探究4. 快速DCT变换DCT变换编码C语言#include #include #include #include #define PI 3.1415926#define CLK_TCK CLOCKS_PER_SECint N;void DCT(double *f,double *F)int n,m,x;double *dTemp = new doubleN*N;/中间矩阵double *coff = new doubleN*N;/变换系数coff0 = 1/sqrt(N);for( m=1; mN; m+ )coffm = sqrt(2)/sqrt(N);memset( dTemp, 0, sizeof(double)*N*N );memset( F, 0, sizeof(double)*N*N );/一维变换for(n=0;nN;n+)for(m=0;mN;m+)for(x=0;xN;x+)dTempm*N+n += fx*N+n * coffm * cos( (2*x+1) * PI * m/(2*N) );/第二次一维变换for(m=0;mN;m+)cfor(n=0;nN;n+)for(x=0;xN;x+)Fm*N+n += dTempm*N+x * coffn * cos( (2*x+1) * PI * n/(2*N) );delete dTemp;delete coff;void iDCT(double *f,double *F)int m,y,x;double *dTemp=new doubleN*N;/中间矩阵double *coff=new doubleN*N;/变换系数coff0=1/sqrt(N);for(m=1;mN;m+)coffm=sqrt(2)/sqrt(N);memset(dTemp,0,sizeof(double)*N*N);memset(F,0,sizeof(double)*N*N);/一维变换for(x=0;xN;x+)for(y=0;yN;y+)for(m=0;mN;m+)dTempx*N+y+=Fx*N+m*coffm*cos(2*y+1)*PI*m/(2*N);/第二次一维变换for(y=0;yN;y+)for(x=0;xN;x+)for(m=0;mN;m+)Fx*N+y+=dTempm*N+y*coffm*cos(2*x+1)*PI*m/(2*N);delete dTemp;delete coff;int main()clock_t start,end;start=clock();int i;long L;printf(变换维数:);scanf(%d,&N);double *f=new doubleN*N;/初始矩阵double *F=new doubleN*N;/变换后输出矩阵memset(F,0,sizeof(double)*N*N);/初始化为0for(i=0;iN*N;i+)printf(f%d%d:,i/N,i%N);scanf(%lf,&fi);printf(循环次数:);scanf(%d,&L);/输出初始矩阵printf(变换前:n);for(i=1;i=N*N;i+)printf(%ft,fi-1);if(i%N=0)printf(n);for(i=0;iL;i+)DCT(f,F);/变换/输出变换后矩阵printf(变换后:n);for(i=1;i=N*N;i+)printf(%ft,Fi-1);if(i%N=0)printf(n);for(i=0;iL;i+)iDCT(f,F);/输出反变换后矩阵printf(反变换后:n);for(i=1;i=N*N;i+)printf(%ft,fi-1);if(i%N=0)printf(n);/printf(n);delete f;delete F;end=clock();printf(耗时:%fn,(double)(end-start)/CLK_TCK);return 0;=MPEG4中逆DCT变换一旦DCT系数Fuv被恢复，那么就可以用逆DCT变换来获得逆变换值fyx，这些只要被饱和到-256fyx255。对短头格式，由于不存在隔行模式，因此全部用帧 DCT变换，就是一般的情况。非短头格式时，如果使用隔行模式，并且dct_type等于1，此时使用场DCT变换，场DCT变换的算法同帧DCT变换完全一样，只是输出的时候需要将按场组织的宏块转换为按帧组织的宏块。下面简单介绍一下DCT变换和逆变换的过程。矩阵大小为NxN的二维DCT变换为：u, v, x, y = 0, 1, 2, N-1其中x, y 是原始域中的空间坐标u, v 是变换域中的空间坐标逆DCT变换定义为：如果每个像素为n比特，则DCT变换的输入为n+1比特，逆DCT变换的输出为n+1比特。DCT变换后的DCT系数的长度为n+4比特，动态范围为-2n+3:+2n+3-1。对我们来说这里的n等于8。NxN的逆DCT变换的实现必须符合IEEE的关于8x8的逆DCT变换的实现的标准，即IEEE Standard Specification for the Implementations of 8 by 8 Inverse Discrete Cosine Transform, Std 1180-1990, December6, 1990，不过有下列修改：1) IEEE规范中的3.2小节的item(1)中最后一句话被替换为：2) IEEE规范中的3.3小节的text被替换为：3) Let F be the set of 4096 blocks Biyx (i=0.4095) defined as follows :a) Bi00 = i - 2048b) Bi77 = 1 if Bi00 is even, Bi77 = 0 if Bi00 is oddc) All other coefficients Biyx other than Bi00 and Bi77 are equal to 0For each block Biyx that belongs to set F defined above, an IDCT that claims to be compliant shall output a block fyx that as a peak error of 1 or less compared to the reference saturated mathematical integer-number IDCT f(x,y). In other words, | fyx - f(x,y)| shall be 。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 只要搞清楚低频和高频是怎么来的和如何合并的即可。 2 DCT简介 DCT全名为Discrete Cosine Transform，中文名为离散余弦变换。在众人皆知的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算法的第一阶级变换代码： For I = 0 To 3 J = 7 - I Y(I) = X(I) + X(J) Y(J) = X(I) - X(J) Next I 设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更好！（具体原因暂不讨论） 5 DCT的改进下面的代码是我对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 -第一阶级 For I = 0 To 3 J = 7 - I Y(J) = X(J) - X(I) Y(I) = X(I) + Fix(Y(J) / 2) Next I -第二阶级 For H = 0 To 4 Step 4 For I = 0 To 1 J = 3 - I X(J + H) = Y(J + H) - Y(I + H) X(I + H) = Y(I + H) + Fix(X(J + H) / 2) Next I Next H -第三阶级 For I = 0 To 6 Step 2 Y(I + 1) = X(I + 1) - X(I) Y(I) = X(I) + Fix(Y(I + 1) / 2) Next I -预测 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丁贵广，计文平，郭宝龙 Visual C+6.0数字图像编码 p44,p57,p170快速DCT变换仿效FFT的FDCT方法有与DCT无关的复数运算部分，选用代数分解法可以降低运算量，达到高速运算的目的。代数分解法实现如下：对一维DCT表达式直接展开，寻找各点表达式中共同项，仿FFT蝶形关系，将表达式中的共同项作为下一级节点，依次进行多次，最后得到变换结果。一、DCT部分例子：Define cos(n*pi/16) CnF(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 C1 0.9808#define C2 0.9239#define C3 0.8315#define C4 0.7071#define C5 0.5556#define C6 0.3827#define C7 0.1951/先做行DCTvoid fdctrow(double *blk)double S07,S16,S25,S34,S0734,S1625;double D07,D16,D25,D34,D0734,D1625;S07=blk0+blk7;S16=blk1+blk6;S25=blk2+blk5;S34=blk3+blk4;S0734=S07+S34;S1625=S16+S25;D07=blk0-blk7; D16=blk1-blk6;D25=blk2-blk5;D34=blk3-blk4;D0734=S07-S34;D1625=S16-S25;blk0=0.5*(C4*(S0734+S1625);blk1=0.5*(C1*D07+C3*D16+C5*D25+C7*D34);blk2=0.5*(C2*D0734+C6*D1625);blk3=0.5*(C3*D07-C7*D16-C1*D25-C5*D34);blk4=0.5*(C4*(S0734-S1625);blk5=0.5*(C5*D07-C1*D16+C7*D25+C3*D34);blk6=0.5*(C6*D0734-C2*D1625);blk7=0.5*(C7*D07-C5*D16+C3*D25-C1*D34);/再做列DCTvoid fdctcol(double *blk)double S07,S16,S25,S34,S0734,S1625;double D07,D16,D25,D34,D0734,D1625;S07=blk0*8+blk7*8;S16=blk1*8+blk6*8;S25=blk2*8+blk5*8;S34=blk3*8+blk4*8;S0734=S07+S34;S1625=S16+S25;D07=blk0*8-blk7*8; D16=blk1*8-blk6*8;D25=blk2*8-blk5*8;D34=blk3*8-blk4*8;D0734=S07-S34;D1625=S16-S25;blk0*8=0.5*(C4*(S0734+S1625);blk1*8=0.5*(C1*D07+C3*D16+C5*D25+C7*D34);blk2*8=0.5*(C2*D0734+C6*D1625);blk3*8=0.5*(C3*D07-C7*D16-C1*D25-C5*D34);blk4*8=0.5*(C4*(S0734-S1625);blk5*8=0.5*(C5*D07-C1*D16+C7*D25+C3*D34);blk6*8=0.5*(C6*D0734-C2*D1625);blk7*8=0.5*(C7*D07-C5*D16+C3*D25-C1*D34);void fdct(double *block) int i; for (i=0; i8; i+) fdctrow(block+8*i); for (i=0; i8; i+) fdctcol(block+i);二、IDCT部分图片来源：G. G. Pechanek, C. W. Kurak, C. J. Glossner, C. H. L. Moller, and S. J. Walsh IBM Microelectronics Division, Research Triangle Park, N.C. M.F.A.S.T.: A HIGHLY PARALLEL SINGLE CHIP DSP WITH A 2D IDCT EXAMPLE图中未给出系数，需自行算出，仿上述DCT方法直接展开表达式搜寻规律即可。编程如下：#define C1 0.9808#define C2 0.9239#define C3 0.8315#define C4 0.7071#define C5 0.5556#define C6 0.3827#define C7 0.1951/对行做DCTvoid idctrow(double *blk)double tmp16;/first steptmp0=blk0*C4+blk2*C2;tmp1=blk4*C4+blk6*C6;tmp2=blk0*C4+blk2*C6;tmp3=-blk4*C4-blk6*C2;tmp4=blk0*C4-blk2*C6;tmp5=-blk4*C4+blk6*C2;tmp6=blk0*C4-blk2*C2;tmp7=blk4*C4-blk6*C6;tmp8=blk1*C7-blk3*C5;tmp9=blk5*C3-blk7*C1;tmp10=blk1*C5-blk3*C1;tmp11=blk5*C7+blk7*C3;tmp12=blk1*C3-blk3*C7;tmp13=-blk5*C1-blk7*C5;tmp14=blk1*C1+blk3*C3;tmp15=blk5*C5+blk7*C7;/second steptmp0=0.5*(tmp0+tmp1);tmp1=0.5*(tmp2+tmp3);tmp2=0.5*(tmp4+tmp5);tmp3=0.5*(tmp6+tmp7);tmp4=0.5*(tmp8+tmp9);tmp5=0.5*(tmp10+tmp11);tmp6=0.5*(tmp12+tmp13);tmp7=0.5*(tmp14+tmp15);/third stepblk0=tmp0+tmp7;blk1=tmp1+tmp6;blk2=tmp2+tmp5;blk3=tmp3+tmp4;blk4=tmp3-tmp4;blk5=tmp2-tmp5;blk6=tmp1-tmp6;blk7=tmp0-tmp7;/*blk0=0.5*(Y0C4+Y2C2+Y4C4+Y6C6+Y1C1+Y3C3+Y5C5+Y7C7);blk1=0.5*(Y0C4+Y2C6-Y4C4-Y6C2+Y1C3-Y3C7-Y5C1-Y7C5);blk2=0.5*(Y0C4-Y2C6-Y4C4+Y6C2+Y1C5-Y3C1+Y5C7+Y7C3);blk3=0.5*(Y0C4-Y2C2+Y4C4-Y6C6+Y1C7-Y3C5+Y5C3-Y7C1);blk4=0.5*(Y0C4-Y2C2+Y4C4-Y6C6-Y1C7+Y3C5-Y5C3+Y7C1);blk5=0.5*(Y0C4-Y2C6-Y4C4+Y6C2-Y1C5+Y3C1-Y5C7-Y7C3);blk6=0.5*(Y0C4+Y2C6-Y4C4-Y6C2-Y1C3+Y3C7+Y5C1+Y7C5);blk7=0.5*(Y0C4+Y2C2+Y4C4+Y6C6-Y1C1-Y3C3-Y5C5-Y7C7);*/在对列做DCTvoid idctcol(double *blk)double tmp16;/first steptmp0=blk0*8*C4+blk2*8*C2;tmp1=blk4*8*C4+blk6*8*C6;tmp2=blk0*8*C4+blk2*8*C6;tmp3=-blk4*8*C4-blk6*8*C2;tmp4=blk0*8*C4-blk2*8*C6;tmp5=-blk4*8*C4+blk6*8*C2;tmp6=blk0*8*C4-blk2*8*C2;tmp7=blk4*8*C4-blk6*8*C6;tmp8=blk1*8*C7-blk3*8*C5;tmp9=blk5*8*C3-blk7*8*C1;tmp10=blk1*8*C5-blk3*8*C1;tmp11=blk5*8*C7+blk7*8*C3;tmp12=blk1*8*C3-blk3*8*C7;tmp13=-blk5*8*C1-blk7*8*C5;tmp14=blk1*8*C1+blk3*8*C3;tmp15=blk5*8*C5+blk7*8*C7;/second steptmp0=0.5*(tmp0+tmp1);tmp1=0.5*(tmp2+tmp3);tmp2=0.5*(tmp4+tmp5);tmp3=0.5*(tmp6+tmp7);tmp4=0.5*(tmp8+tmp9);tmp5=0.5*(tmp10+tmp11);tmp6=0.5*(tmp12+tmp13);tmp7=0.5*(tmp14+tmp15);/third stepblk0*8=tmp0+tmp7;blk1*8=tmp1+tmp6;blk2*8=tmp2+tmp5;blk3*8=tmp3+tmp4;blk4*8=tmp3-tmp4;blk5*8=tmp2-tmp5;blk6*8=tmp1-tmp6;blk7*8=tmp0-tmp7;/*blk0*8=0.5*(Y0C4+Y2C2+Y4C4+Y6C6+Y1C1+Y3C3+Y5C5+Y7C7);blk1*8=0.5*(Y0C4+Y2C6-Y4C4-Y6C2+Y1C3-Y3C7-Y5C1-Y7C5);blk2*8=0.5*(Y0C4-Y2C6-Y4C4+Y6C2+Y1C5-Y3C1+Y5C7+Y7C3);blk3*8=0.5*(Y0C4-Y2C2+Y4C4-Y6C6+Y1C7-Y3C5+Y5C3-Y7C1);blk4*8=0.5*(Y0C4-Y2C2+Y4C4-Y6C6-Y1C7+Y3C5-Y5C3+Y7C1);blk5*8=0.5*(Y0C4-Y2C6-Y4C4+Y6C2-Y1C5+Y3C1-Y5C7-Y7C3);blk6*8=0.5*(Y0C4+Y2C6-Y4C4-Y6C2-Y1C3+Y3C7+Y5C1+Y7C5);blk7*8=0.5*(Y0C4+Y2C2+Y4C4+Y6C6-Y1C1-Y3C3-Y5C5-Y7C7);*/void idct(double *block) int i; for (i=0; i8; i+) idctrow(block+8*i); for (i=0; i8; i+) idctcol(block+i);经简单测试，上面两种算法是简单矩阵实现DCT速度的100倍以上。