离散余弦变换(FDCT和IDCT)

最新推荐文章于 2023-07-01 20:02:59 发布

qiuchangyong

最新推荐文章于 2023-07-01 20:02:59 发布

阅读量1.3w

点赞数

分类专栏：数学与图形学文章标签： c output 语言

数学与图形学专栏收录该内容

24 篇文章 6 订阅

订阅专栏

原文地址：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);
}

qiuchangyong

关注

0
点赞
踩
11

收藏

觉得还不错? 一键收藏
0
评论
离散余弦变换(FDCT和IDCT)

原文地址：http://hi.baidu.com/alphame/blog/item/2a4e52b7aa50a1f130add1a1.html1.为什么要进行变换？图像数据一般有较强的相关性，若所选用的正交矢量空间的基矢量与图像本身的主要特征相近，在该正交矢量空间中描述图像数据则会变得更简单。经过正交变换，会把原来分散在原空间的图像数据在新的坐标空间中
复制链接

扫一扫