.DCTb变换

Contents

1.       DCT变换编码C语言

2.     MPEG4中逆DCT变换

3.     DCT变换探究

4.     快速DCT变换

 

 

 

DCT变换编码C语言

#include <memory.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#define PI 3.1415926
#define CLK_TCK CLOCKS_PER_SEC
int N;
void DCT(double *f,double *F)
{
  int n,m,x;
 double *dTemp = new double[N*N];//中间矩阵
 double *coff = new double[N*N];//变换系数
 coff[0] = 1/sqrt(N);
 for( m=1; m<N; m++ )
  coff

­ = sqrt(2)/sqrt(N);
 memset( dTemp, 0, sizeof(double)*N*N );
 memset( F, 0, sizeof(double)*N*N );

 //一维变换
 for(n=0;n<N;n++)
  for(m=0;m<N;m++)
   for(x=0;x<N;x++)
    dTemp[m*N+n] += f[x*N+n] * coff

­ * cos( (2*x+1) * PI * m/(2*N) );

 //第二次一维变换
 for(m=0;m<N;m++)c
  for(n=0;n<N;n++)
   for(x=0;x<N;x++)
    F[m*N+n] += dTemp[m*N+x] * coff[n] * 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 double[N*N];//中间矩阵
 double *coff=new double[N*N];//变换系数
 coff[0]=1/sqrt(N);
 for(m=1;m<N;m++)
  coff

­=sqrt(2)/sqrt(N);
 memset(dTemp,0,sizeof(double)*N*N);
 memset(F,0,sizeof(double)*N*N);

 //一维变换
 for(x=0;x<N;x++)
  for(y=0;y<N;y++)
   for(m=0;m<N;m++)
    dTemp[x*N+y]+=F[x*N+m]*coff

­*cos((2*y+1)*PI*m/(2*N));

 //第二次一维变换
 for(y=0;y<N;y++)
  for(x=0;x<N;x++)
   for(m=0;m<N;m++)
    F[x*N+y]+=dTemp[m*N+y]*coff

­*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 double[N*N];//初始矩阵
 double *F=new double[N*N];//变换后输出矩阵 

memset(F,0,sizeof(double)*N*N);//初始化为0

 for(i=0;i<N*N;i++)
 {
  printf("f[%d][%d]:",i/N,i%N);
  scanf("%lf",&f[i]);
 }

 printf("循环次数:");
 scanf("%d",&L);

//输出初始矩阵
 printf("变换前:\n");
 for(i=1;i<=N*N;i++)
 {
  printf("%f\t",f[i-1]);
  if(i%N==0)
   printf("\n");
 }

 for(i=0;i<L;i++)
  DCT(f,F);//变换
//输出变换后矩阵
 printf("变换后:\n");

 for(i=1;i<=N*N;i++)
 {
  printf("%f\t",F[i-1]);
  if(i%N==0)
   printf("\n");
 } 

for(i=0;i<L;i++)
 iDCT(f,F);
 //输出反变换后矩阵
 printf("反变换后:\n");
 for(i=1;i<=N*N;i++)
 {  printf("%f\t",f[i-1]);
  if(i%N==0)
   printf("\n");
 }
 //printf("\n");

 delete []f;
 delete []F;
 end=clock();
 printf("耗时:%f\n",(double)(end-start)/CLK_TCK);
 return 0; 
}

 

===================================================================

MPEG4中逆DCT变换

一旦DCT系数F[u][v]被恢复，那么就可以用逆DCT变换来获得逆变换值f[y][x]，这些只要被饱和到-256≤f[y][x]≤255。

对短头格式，由于不存在隔行模式，因此全部用帧 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比特，动态范围为[-2ⁿ⁺³:+2ⁿ⁺³-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, December 6, 1990，不过有下列修改：

1)                IEEE规范中的3.2小节的item(1)中最后一句话被替换为：<<Data sets of 1 000 000 (one million) blocks each should be generated for (L=256, H=255), (L=H=5) and (L=384, H=383). >>

2)                IEEE规范中的3.3小节的text被替换为：<<For any pixel location, the peak error shall not exceed 2 in magnitude.  There is no other accuracy requirement  for this test.>>

3)                Let F be the set of 4096 blocks Bi[y][x] (i=0..4095) defined as follows :

a) Bi[0][0] = i - 2048

                                   b) Bi[7][7] = 1 if  Bi[0][0] is even, Bi[7][7] = 0 if  Bi[0][0] is odd

c) All other coefficients Bi[y][x] other than Bi[0][0] and Bi[7][7] are equal to 0

For each block Bi[y][x] that belongs to set F defined above, an IDCT that claims to be compliant shall output a block f[y][x] that as a peak error of 1 or less compared to the reference saturated mathematical integer-number IDCT fíí(x,y).  In other words, | f[y][x] - fíí(x,y)| shall be <= 1 for all x and y.

NOTE 1  lause  2.3 Std 1180-1990 “Considerations of Specifying IDCT Mismatch Errors” requires the specification of periodic intra-picture coding in order to control the accumulation of mismatch errors. Every macroblock is required to be refreshed before it is coded 132 times as predictive macroblocks. Macroblocks in B-pictures (and skipped macroblocks in P-pictures) are excluded from the counting because they do not lead to the accumulation of mismatch errors. This requirement is the same as indicated in 1180-1990 for visual telephony according to ITU-T Recommendation H.261.

NOTE 2  Whilst the IEEE IDCT standard mentioned above is a necessary condition for the satisfactory implementation of the IDCT function it should be understood that this is not sufficient. In particular attention is drawn to the following sentence from subclause 5.4: “Where arithmetic precision is not specified, such as the calculation of the IDCT, the precision shall be sufficient so that significant errors do not occur in the final integer values.”

逆DCT变换的过程这里不再详述，需要实现这个的可以去参考这个标准。

在实际应用中一般通过两次1-D IDCT变换来完成2-D IDCT变换，这种方法通常被称为行－列法。

一般来说，后者在结构上的对称性更好，并且可以重复使用硬件资源，所以在我们的芯片设计选用一种行-列法来进行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}
只要搞清楚低频和高频是怎么来的和如何合并的即可。

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) 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  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

//先做行DCT
void fdctrow(double *blk)
{ 
double S07,S16,S25,S34,S0734,S1625;
double D07,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);

}

//再做列DCT

void fdctcol(double *blk)
{
double S07,S16,S25,S34,S0734,S1625;
double D07,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);

}

void fdct(double *block)
{
  int i;

  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, 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

//对行做DCT

void idctrow(double *blk)
{ double tmp[16];

 //first step
 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;

 //second step
 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]);

 //third step
 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);
 */

}

//在对列做DCT

void idctcol(double *blk)
{
 double tmp[16];

 //first step
 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;

 //second step
 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]);

 //third step
 blk[0*8]=tmp[0]+tmp[7];
 blk[1*8]=tmp[1]+tmp[6];
 blk[2*8]=tmp[2]+tmp[5];
 blk[3*8]=tmp[3]+tmp[4];
 blk[4*8]=tmp[3]-tmp[4];
 blk[5*8]=tmp[2]-tmp[5];
 blk[6*8]=tmp[1]-tmp[6];
 blk[7*8]=tmp[0]-tmp[7];

 /*
 blk[0*8]=0.5*(Y0C4+Y2C2+Y4C4+Y6C6+Y1C1+Y3C3+Y5C5+Y7C7);
 blk[1*8]=0.5*(Y0C4+Y2C6-Y4C4-Y6C2+Y1C3-Y3C7-Y5C1-Y7C5);
 blk[2*8]=0.5*(Y0C4-Y2C6-Y4C4+Y6C2+Y1C5-Y3C1+Y5C7+Y7C3);
 blk[3*8]=0.5*(Y0C4-Y2C2+Y4C4-Y6C6+Y1C7-Y3C5+Y5C3-Y7C1);

 blk[4*8]=0.5*(Y0C4-Y2C2+Y4C4-Y6C6-Y1C7+Y3C5-Y5C3+Y7C1);
 blk[5*8]=0.5*(Y0C4-Y2C6-Y4C4+Y6C2-Y1C5+Y3C1-Y5C7-Y7C3);
 blk[6*8]=0.5*(Y0C4+Y2C6-Y4C4-Y6C2-Y1C3+Y3C7+Y5C1+Y7C5);
 blk[7*8]=0.5*(Y0C4+Y2C2+Y4C4+Y6C6-Y1C1-Y3C3-Y5C5-Y7C7);
 */

}

void idct(double *block)
{
  int i;

  for (i=0; i<8; i++)
    idctrow(block+8*i);

  for (i=0; i<8; i++)
    idctcol(block+i);
}

经简单测试，上面两种算法是简单矩阵实现DCT速度的100倍以上。