离散余弦变换(FDCT和IDCT)

转载 2012年03月23日 11:53:33

原文地址:http://hi.baidu.com/alphame/blog/item/2a4e52b7aa50a1f130add1a1.html

1.为什么要进行变换?

       图像数据一般有较强的相关性,若所选用的正交矢量空间的基矢量与图像本身的主要特征相近,在该正交矢量空间中描述图像数据则会变得更简单。
        经过正交变换,会把原来分散在原空间的图像数据在新的坐标空间中得到集中。对于大多数图像,大量变换系数很小,只要删除接近于零的系数,并且对较小的系数进行粗量化,而保留包含图像主要信息的系数,以此进行压缩编码。
       在重建图像进行解码时,所损失的将是一些不重要的信息,几乎不会引起图像的失真。

2.子图像尺寸选择
         在变换编码中,首先要将图像数据分割成子图像,然后对子图像数据块实施某种变换,如DCT变换,那么子图像尺寸取多少好呢?根据实践证明子图像尺寸取4×4、8×8、16×16适合作图像的压缩,这是因为:
        <1> 如果子图像尺寸取得太小,虽然计算速度快,实现简单,但压缩能力有一定的限制。
        <2> 如果子图像尺寸取得太大,虽然去相关效果变好,因为象DFT、DCT等正弦型变换均具有渐近最佳性,但也渐趋饱和。若尺寸太大,由于图像本身的相关性很小,反而使其压缩效果不明显,而且增加了计算的复杂性。

3. 数学公式:

FDCT:
                                              7       7                      2*x+1                   2*y+1
F(u,v) = alpha(u)*alpha(v)* sum sum f(x,y) * cos (------- *u*PI)* cos (--------- *v*PI)
                                            x=0   y=0                    16                          16

u,v = 0,1,...,7

                  { 1/sqrt(8) (u==0)
alpha(u) = {
                 { 1/2            (u!=0)

IDCT:
               7      7                                                  2*x+1                    2*y+1
f(x,y) = sum sum alpha(u)*alpha(v)*F(u,v)*cos (------- *u*PI)* cos (------ *v*PI)
              u=0 v=0                                                 16                         16

x,y=0,1...7

    上述公式是二维的, 二维离散余弦变换可以分解成一维DCT来实现,也就是说可以先行,后列或者先列后行。

                              7                   2*x+1
F(u) = alpha(u)* sum f(x) * cos (------- *u*PI)
                            x=0                   16

下面是C语言的实现,输入都是8x8的整数矩阵,每个元素在-128到127之间。FDCT0和IDCT0是用来参考的,DCT1和DCT2是两种实现方式,DCT1规则一些,DCT2使用了较少的乘法运算。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

static const double PI=3.14159265358979323;
const int C1_16 = 4017; // cos( pi/16) x4096
const int C2_16 = 3784; // cos(2pi/16) x4096
const int C3_16 = 3406; // cos(3pi/16) x4096
const int C4_16 = 2896; // cos(4pi/16) x4096
const int C5_16 = 2276; // cos(5pi/16) x4096
const int C6_16 = 1567; // cos(6pi/16) x4096
const int C7_16 = 799; // cos(7pi/16) x4096

void display(int *DataIn){
int i ;
for (i=0;i<64;i++){
    printf("%d ",DataIn[i]);
    if(i%8 == 7 ) printf("\n");
}
}

void FDCT0(int *BlockIn, int *BlockOut){
int u,v;
int i,x,y;
double tmp;
for(i=0;i<64;i++){
    u= i % 8;
    v= i / 8;
    tmp = 0;
    for(y=0;y<8;y++){
      for(x=0;x<8;x++){
        tmp += BlockIn[y*8+x]*cos(((2*x+1)*u*PI)/16)*cos(((2*y+1)*v*PI)/16);
      }
    }
    if (v==0) tmp *= sqrt(2)/2;
    if (u==0) tmp *= sqrt(2)/2;
    BlockOut[i] = tmp/4;
}
display(BlockOut);
}

void FDCT1(int *data){
    int i;
    int S07,S16,S25,S34,S0734,S1625;
    int D07,D16,D25,D34,D0734,D1625;
    for(i=0;i<8;i++){
      S07=(data[0]+data[7]); S16=(data[1]+data[6]);
      S25=(data[2]+data[5]); S34=(data[3]+data[4]);
      D07=(data[0]-data[7]); D16=(data[1]-data[6]);
      D25=(data[2]-data[5]); D34=(data[3]-data[4]);

      S0734=S07+S34; S1625=S16+S25;
      D0734=S07-S34; D1625=S16-S25;

      *data++=(C4_16*(S0734+S1625))>>11;
      *data++=(C1_16*D07+C3_16*D16+C5_16*D25+C7_16*D34)>>11;
      *data++=(C2_16*D0734+C6_16*D1625)>>11;
      *data++=(C3_16*D07-C7_16*D16-C1_16*D25-C5_16*D34)>>11;
      *data++=(C4_16*(S0734-S1625))>>11;
      *data++=(C5_16*D07-C1_16*D16+C7_16*D25+C3_16*D34)>>11;
      *data++=(C6_16*D0734-C2_16*D1625)>>11;
      *data++=(C7_16*D07-C5_16*D16+C3_16*D25-C1_16*D34)>>11;
    }
    data -= 64;
    for(i=0;i<8;i++){
      S07=data[0]+data[56];
      S16=data[8]+data[48];
      S25=data[16]+data[40];
      S34=data[24]+data[32];
      D07=data[0]-data[56];
      D16=data[8]-data[48];
      D25=data[16]-data[40];
      D34=data[24]-data[32];

      S0734=S07+S34; S1625=S16+S25;
      D0734=S07-S34; D1625=S16-S25;

      data[0]=(C4_16*(S0734+S1625))>>15;
      data[8]=(C1_16*D07+C3_16*D16+C5_16*D25+C7_16*D34)>>15;
      data[16]=(C2_16*D0734+C6_16*D1625)>>15;
      data[24]=(C3_16*D07-C7_16*D16-C1_16*D25-C5_16*D34)>>15;
      data[32]=(C4_16*(S0734-S1625))>>15;
      data[40]=(C5_16*D07-C1_16*D16+C7_16*D25+C3_16*D34)>>15;
      data[48]=(C6_16*D0734-C2_16*D1625)>>15;
      data[56]=(C7_16*D07-C5_16*D16+C3_16*D25-C1_16*D34)>>15;
      data++;
    }
    data -=8;
    display(data);
}
void FDCT2(int *data){
static const int c1=1004 ; /* cos(pi/16)<<10 */
static const int s1=200 ; /* sin(pi/16)<<10 */
static const int c3=851 ; /* cos(3pi/16)<<10 */
static const int s3=569 ; /* sin(3pi/16)<<10 */
static const int r2c6=554; /* sqrt(2)*cos(6pi/16)<<10 */
static const int r2s6=1337;
static const int r2=181;   /* sqrt(2)<<7 */
int x0,x1,x2,x3,x4,x5,x6,x7,x8;
int i;
for(i=0;i<8;i++){
    x0=data[0], x1=data[1], x2=data[2], x3=data[3],
    x4=data[4], x5=data[5], x6=data[6], x7=data[7];
    /* Stage 1 */
    x8=x7+x0; x0-=x7; x7=x1+x6; x1-=x6;
    x6=x2+x5; x2-=x5; x5=x3+x4; x3-=x4;
    /* Stage 2 */
    x4=x8+x5; x8-=x5; x5=x7+x6; x7-=x6;
    x6=c1*(x1+x2); x2=(-s1-c1)*x2+x6; x1=(s1-c1)*x1+x6;
    x6=c3*(x0+x3); x3=(-s3-c3)*x3+x6; x0=(s3-c3)*x0+x6;
    /* Stage 3 */
    x6=x4+x5; x4-=x5;
    x5=r2c6*(x7+x8); x7=(-r2s6-r2c6)*x7+x5; x8=(r2s6-r2c6)*x8+x5;
    x5=x0+x2;x0-=x2; x2=x3+x1; x3-=x1;
    /* Stage 4, round, and output */
    *data++=x6;
    *data++=(x2+x5+512)>>10;
    *data++=(x8+512)>>10;
    *data++=(x3*r2+65536)>>17;
    *data++=x4;
    *data++=(x0*r2+65536)>>17;
    *data++=(x7+512)>>10;
    *data++=(x2-x5+512)>>10;
}
data -=64;
for(i=0;i<8;i++){
    x0=data[0], x1=data[8], x2=data[16], x3=data[24],
    x4=data[32], x5=data[40], x6=data[48], x7=data[56];
    /* Stage 1 */
    x8=x7+x0; x0-=x7; x7=x1+x6; x1-=x6;
    x6=x2+x5; x2-=x5; x5=x3+x4; x3-=x4;
    /* Stage 2 */
    x4=x8+x5; x8-=x5; x5=x7+x6; x7-=x6;
    x6=c1*(x1+x2); x2=(-s1-c1)*x2+x6; x1=(s1-c1)*x1+x6;
    x6=c3*(x0+x3); x3=(-s3-c3)*x3+x6; x0=(s3-c3)*x0+x6;
    /* Stage 3 */
    x6=x4+x5; x4-=x5;
    x5=r2c6*(x7+x8); x7=(-r2s6-r2c6)*x7+x5; x8=(r2s6-r2c6)*x8+x5;
    x5=x0+x2;x0-=x2; x2=x3+x1; x3-=x1;
    /* Stage 4, round, and output */
    data[0]=x6>>3;
    data[8]=(x2+x5+512)>>13;
    data[16]=(x8+512)>>13;
    data[24]=(x3*r2+65536)>>20;
    data[32]=x4>>3;
    data[40]=(x0*r2+65536)>>20;
    data[48]=(x7+512)>>13;
    data[56]=(x2-x5+512)>>13;
    data++;
}
data -=8;
display(data);
}
void IDCT0(int *BlockIn, int*BlockOut){
int x,y;
int i,u,v;
double tmp, tmp1;
for(i=0;i<64;i++){
    x= i % 8;
    y= i / 8;
    tmp = 0;
    for(u=0;u<8;u++){
      for(v=0;v<8;v++){
        tmp1 = BlockIn[v*8+u]*cos(((2*x+1)*u*PI)/16)*cos(((2*y+1)*v*PI)/16);
        if (u==0) tmp1 *= sqrt(2)/2;
        if (v==0) tmp1 *= sqrt(2)/2;
tmp += tmp1;
      }
    }
    BlockOut[i] = tmp / 4;
}
display(BlockOut);
}
void IDCT1(int *data){
int i;
int s0,s1,s2,s3,s4,s5,s6,s7;
int t0,t1,t2,t3,t4,t5,t6,t7;
for(i=0;i<8;i++) {
    s0 = (data[0] + data[4]) * C4_16;
    s1 = (data[0] - data[4]) * C4_16;
    s3 = (data[2] * C2_16) + (data[6] * C6_16);
    s2 = (data[2] * C6_16) - (data[6] * C2_16);
    s7 = (data[1] * C1_16) + (data[7] * C7_16);
    s4 = (data[1] * C7_16) - (data[7] * C1_16);
    s6 = (data[5] * C5_16) + (data[3] * C3_16);
    s5 = (data[5] * C3_16) - (data[3] * C5_16);

    t0 = s0 + s3; t1 = s1 + s2; t2 = s1 - s2; t3 = s0 - s3;
    t4 = s4 + s5; t5 = s4 - s5; t6 = s7 - s6; t7 = s7 + s6;

    s6 = (t5 + t6) * 181 / 256; // 1/sqrt(2)
    s5 = (t6 - t5) * 181 / 256; // 1/sqrt(2)

    *data++ = (t0 + t7) >> 11; *data++ = (t1 + s6) >> 11;
    *data++ = (t2 + s5) >> 11; *data++ = (t3 + t4) >> 11;
    *data++ = (t3 - t4) >> 11; *data++ = (t2 - s5) >> 11;
    *data++ = (t1 - s6) >> 11; *data++ = (t0 - t7) >> 11;
}
data -= 64;
for(i=0;i<8;i++){
    s0 = (data[ 0] + data[32]) * C4_16;
    s1 = (data[ 0] - data[32]) * C4_16;
    s3 = data[16] * C2_16 + data[48] * C6_16;
    s2 = data[16] * C6_16 - data[48] * C2_16;
    s7 = data[ 8] * C1_16 + data[56] * C7_16;
    s4 = data[ 8] * C7_16 - data[56] * C1_16;
    s6 = data[40] * C5_16 + data[24] * C3_16;
    s5 = data[40] * C3_16 - data[24] * C5_16;

    t0 = s0 + s3; t1 = s1 + s2; t2 = s1 - s2; t3 = s0 - s3;
    t4 = s4 + s5; t5 = s4 - s5; t6 = s7 - s6; t7 = s6 + s7;

    s5 = (t6 - t5) * 181 / 256; // 1/sqrt(2)
    s6 = (t5 + t6) * 181 / 256; // 1/sqrt(2)

    data[ 0] = ((t0 + t7) >> 15); data[56] = ((t0 - t7) >> 15);
    data[ 8] = ((t1 + s6) >> 15); data[48] = ((t1 - s6) >> 15);
    data[16] = ((t2 + s5) >> 15); data[40] = ((t2 - s5) >> 15);
    data[24] = ((t3 + t4) >> 15); data[32] = ((t3 - t4) >> 15);
    data++;
}
data-=8;
display(data);
}
void IDCT2(int *data)
{
static const int c1=251 ; /* cos(pi/16)<<8 */
static const int s1=50 ; /* sin(pi/16)<<8 */
static const int c3=213 ; /* cos(3pi/16)<<8 */
static const int s3=142 ; /* sin(3pi/16)<<8 */
static const int r2c6=277; /* cos(6pi/16)*sqrt(2)<<9 */
static const int r2s6=669;
static const int r2=181; /* sqrt(2)<<7 */
int i;
int x0,x1,x2,x3,x4,x5,x6,x7,x8;
/* Stage 4 */
for(i=0;i<8;i++){
    x0=data[0]<<9, x1=data[1]<<7, x2=data[2],
    x3=data[3]*r2, x4=data[4]<<9, x5=data[5]*r2,
    x6=data[6], x7=data[7]<<7;
    x8=x7+x1; x1 -= x7;
    /* Stage 3 */
    x7=x0+x4; x0-=x4; x4=x1+x5; x1-=x5; x5=x3+x8; x8-=x3;
    x3=r2c6*(x2+x6);x6=x3+(-r2c6-r2s6)*x6;x2=x3+(-r2c6+r2s6)*x2;
    /* Stage 2 */
    x3=x7+x2; x7-=x2; x2=x0+x6; x0-= x6;
    x6=c3*(x4+x5);x5=(x6+(-c3-s3)*x5)>>6;x4=(x6+(-c3+s3)*x4)>>6;
    x6=c1*(x1+x8);x1=(x6+(-c1-s1)*x1)>>6;x8=(x6+(-c1+s1)*x8)>>6;
    /* Stage 1, rounding and output */
    x7+=512; x2+=512;x0+=512;x3+=512;
    *data++=(x3+x4)>>10; *data++=(x2+x8)>>10;
    *data++=(x0+x1)>>10; *data++=(x7+x5)>>10;
    *data++=(x7-x5)>>10; *data++=(x0-x1)>>10;
    *data++=(x2-x8)>>10; *data++=(x3-x4)>>10;
}
data -= 64;
for(i=0;i<8;i++){
    x0=data[0]<<9, x1=data[8]<<7, x2=data[16],
    x3=data[24]*r2, x4=data[32]<<9, x5=data[40]*r2,
    x6=data[48], x7=data[56]<<7;
    x8=x7+x1; x1 -= x7;
    /* Stage 3 */
    x7=x0+x4; x0-=x4; x4=x1+x5; x1-=x5; x5=x3+x8; x8-=x3;
    x3=r2c6*(x2+x6);x6=x3+(-r2c6-r2s6)*x6;x2=x3+(-r2c6+r2s6)*x2;
    /* Stage 2 */
    x3=x7+x2; x7-=x2; x2=x0+x6; x0-= x6;
    x6=c3*(x4+x5);x5=(x6+(-c3-s3)*x5)>>6;x4=(x6+(-c3+s3)*x4)>>6;
    x6=c1*(x1+x8);x1=(x6+(-c1-s1)*x1)>>6;x8=(x6+(-c1+s1)*x8)>>6;
    /* Stage 1, rounding and output */
    x7+=512; x2+=512;x0+=512;x3+=512;
    data[0]=(x3+x4)>>11; data[8]=(x2+x8)>>11;
    data[16]=(x0+x1)>>11; data[24]=(x7+x5)>>11;
    data[32]=(x7-x5)>>11; data[40]=(x0-x1)>>11;
    data[48]=(x2-x8)>>11; data[56]=(x3-x4)>>11;
    data ++;
}
data -=8;
display(data);
}


main(){
int Data1[64];
int Data2[64];
int i ;
for(i=0;i<64;i++){
    Data1[i]=rand()%256 - 128;
}
printf("\n -- DataIn -- \n");
display(Data1);
printf("\n -- FDCT0 -- \n");
FDCT0(Data1,Data2);

memcpy(Data2,Data1,sizeof(int)*64);
printf("\n -- FDCT1 -- \n");
FDCT1(Data2);

memcpy(Data2,Data1,sizeof(int)*64);
printf("\n -- FDCT2 -- \n");
FDCT2(Data2);

memcpy(Data2,Data1,sizeof(int)*64);
printf("\n -- IDCT0 -- \n");
IDCT0(Data1,Data2);

memcpy(Data2,Data1,sizeof(int)*64);
printf("\n -- IDCT1 -- \n");
IDCT1(Data2);

memcpy(Data2,Data1,sizeof(int)*64);
printf("\n -- IDCT2 -- \n");
IDCT2(Data2);
}

相关文章推荐

快速离散余弦变换代码实现(FDCT)

网上常见的快速离散余弦变换的代码如下:#include using namespace std;#define DCTSIZE 8void FDCT(double* lpBuff){  double ...

离散余弦变换(FDCT和IDCT)

1.为什么要进行变换?        图像数据一般有较强的相关性,若所选用的正交矢量空间的基矢量与图像本身的主要特征相近,在该正交矢量空间中描述图像数据则会变得更简单。          经过...
  • htjx99
  • htjx99
  • 2014年09月16日 12:09
  • 980

FDCT和IDCT的计算公式

DCT变换在信号处理领域非常重要,是语音图像压缩的基础。它有一维和二维两种形式。...

IDCT的原理以及代码分析

代码版本: HM16.2IDCT(逆DCT)变换是DCT变换的逆变换。在HEVC的解码端对图像的像素残差进行IDCT变换,再进行反量化,就可以得到最初的残差值, 原像素 = 预测值 + 残差值。这篇文...

Delphi7高级应用开发随书源码

  • 2003年04月30日 00:00
  • 676KB
  • 下载

H.264编码之IDCT变换公式推导

H.264编码中的DCT变换与反变换的具体细节,可能比较难,再加上实际工作中用的不是很多,很多朋友可能看的并不多,也不是特别细,下面将详细介绍下H.264中的IDCT变换的具体细节和推导过程!...

DCT快速变换

DCT快速变换作者:陈祖尚 (blog:darnshong)下载源代码一、引言DCT变换是数字图像处理中重要的变换,很多重要的图像算法、图像应用都是基于DCT变换的,如JPEG图像编码方式。对于大尺寸...
  • dd167
  • dd167
  • 2006年02月15日 17:41
  • 3602

离散余弦变换 DCT

MPEG采用了Ahmed(一个巨牛的数学家) 等人于70年代提出的离散余弦变换(DCT-Discrete Cosine Transform)压缩算法,降低视频信号的空间冗余度。DCT将运动补偿误差或...
  • crfoxzl
  • crfoxzl
  • 2007年11月03日 14:47
  • 15043

Delphi7高级应用开发随书源码

  • 2003年04月30日 00:00
  • 676KB
  • 下载

离散余弦变换(含源码)

图像处理中常用的正交变换除了傅里叶变换外,还有其他一些有用的正交变换,其中离散余弦就是一种。离散余弦变换表示为DCT( Discrete Cosine Transformation),常用于图像处理和...
  • luoweifu
  • luoweifu
  • 2012年11月23日 09:38
  • 21009
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:离散余弦变换(FDCT和IDCT)
举报原因:
原因补充:

(最多只允许输入30个字)