JPEG 编解码中的 dct 变换

/* 包含头文件 */
#include "dct.h"

#if 1 /* 快速的整数运算版本 */
/* 内部常量定义 */
#define DCTSIZE  8

/* 内部全局变量定义 */
static long enfactor[64] =
{
     65536,  47248,  50159,  55733,  65536,  83411, 121094, 237535,
     47248,  34064,  36162,  40181,  47248,  60136,  87304, 171253,
     50159,  36162,  38390,  42656,  50159,  63840,  92681, 181802,
     55733,  40181,  42656,  47397,  55733,  70935, 102982, 202007,
     65536,  47248,  50159,  55733,  65536,  83411, 121094, 237535,
     83411,  60136,  63840,  70935,  83411, 106162, 154124, 302325,
    121094,  87304,  92681, 102982, 121094, 154124, 223753, 438909,
    237535, 171253, 181802, 202007, 237535, 302325, 438909, 860951,
};

static long defactor[64] =
{
    65536,  90901,  85626,  77062,  65536,  51491,  35467,  18081,
    90901, 126083, 118767, 106888,  90901,  71420,  49195,  25079,
    85626, 118767, 111876, 100686,  85626,  67276,  46340,  23624,
    77062, 106888, 100686,  90615,  77062,  60547,  41705,  21261,
    65536,  90901,  85626,  77062,  65536,  51491,  35467,  18081,
    51491,  71420,  67276,  60547,  51491,  40456,  27866,  14206,
    35467,  49195,  46340,  41705,  35467,  27866,  19195,   9785,
    18081,  25079,  23624,  21261,  18081,  14206,   9785,   4988,
};

/*
  46341  / 65536 = 0.707106781f
  25080  / 65536 = 0.382683433f
  35468  / 65536 = 0.541196100f
  85627  / 65536 = 1.306562965f
  92682  / 65536 = 1.414213562f
  121095 / 65536 = 1.847759065f
  70936  / 65536 = 1.082392200f
  171254 / 65536 = 2.613125930f
 */

/* 函数实现 */
void fdct2d8x8(long *data)
{
    long tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
    long tmp10, tmp11, tmp12, tmp13;
    long z1, z2, z3, z4, z5, z11, z13;
    long *dataptr;
    int  ctr, i;

    for (i=0; i<64; i++) data[i] <<= 1;

    /* Pass 1: process rows. */
    dataptr = data;
    for (ctr=0; ctr<DCTSIZE; ctr++)
    {
        tmp0 = dataptr[0] + dataptr[7];
        tmp7 = dataptr[0] - dataptr[7];
        tmp1 = dataptr[1] + dataptr[6];
        tmp6 = dataptr[1] - dataptr[6];
        tmp2 = dataptr[2] + dataptr[5];
        tmp5 = dataptr[2] - dataptr[5];
        tmp3 = dataptr[3] + dataptr[4];
        tmp4 = dataptr[3] - dataptr[4];
   
        /* Even part */
        tmp10 = tmp0 + tmp3;  /* phase 2 */
        tmp13 = tmp0 - tmp3;
        tmp11 = tmp1 + tmp2;
        tmp12 = tmp1 - tmp2;

        dataptr[0] = tmp10 + tmp11;  /* phase 3 */
        dataptr[4] = tmp10 - tmp11;

        z1 = (tmp12 + tmp13) * 46341 >> 16; /* c4 */
        dataptr[2] = tmp13 + z1;  /* phase 5 */
        dataptr[6] = tmp13 - z1;

        /* Odd part */
        tmp10 = tmp4 + tmp5;  /* phase 2 */
        tmp11 = tmp5 + tmp6;
        tmp12 = tmp6 + tmp7;

        /* The rotator is modified from fig 4-8 to avoid extra negations. */
        z5 = (tmp10 - tmp12) * 25080 >> 16;  /* c6 */
        z2 = (tmp10 * 35468 >> 16) + z5;     /* c2-c6 */
        z4 = (tmp12 * 85627 >> 16) + z5;     /* c2+c6 */
        z3 =  tmp11 * 46341 >> 16;           /* c4 */

        z11 = tmp7 + z3;        /* phase 5 */
        z13 = tmp7 - z3;

        dataptr[5] = z13 + z2;  /* phase 6 */
        dataptr[3] = z13 - z2;
        dataptr[1] = z11 + z4;
        dataptr[7] = z11 - z4;

        dataptr += DCTSIZE;     /* advance pointer to next row */
    }

    /* Pass 2: process columns. */
    dataptr = data;
    for (ctr=0; ctr<DCTSIZE; ctr++)
    {
        tmp0 = dataptr[DCTSIZE * 0] + dataptr[DCTSIZE * 7];
        tmp7 = dataptr[DCTSIZE * 0] - dataptr[DCTSIZE * 7];
        tmp1 = dataptr[DCTSIZE * 1] + dataptr[DCTSIZE * 6];
        tmp6 = dataptr[DCTSIZE * 1] - dataptr[DCTSIZE * 6];
        tmp2 = dataptr[DCTSIZE * 2] + dataptr[DCTSIZE * 5];
        tmp5 = dataptr[DCTSIZE * 2] - dataptr[DCTSIZE * 5];
        tmp3 = dataptr[DCTSIZE * 3] + dataptr[DCTSIZE * 4];
        tmp4 = dataptr[DCTSIZE * 3] - dataptr[DCTSIZE * 4];

        /* Even part */
        tmp10 = tmp0 + tmp3;  /* phase 2 */
        tmp13 = tmp0 - tmp3;
        tmp11 = tmp1 + tmp2;
        tmp12 = tmp1 - tmp2;

        dataptr[DCTSIZE * 0] = tmp10 + tmp11;  /* phase 3 */
        dataptr[DCTSIZE * 4] = tmp10 - tmp11;

        z1 = (tmp12 + tmp13) * 46341 >> 16; /* c4 */
        dataptr[DCTSIZE * 2] = tmp13 + z1;  /* phase 5 */
        dataptr[DCTSIZE * 6] = tmp13 - z1;

        /* Odd part */
        tmp10 = tmp4 + tmp5;  /* phase 2 */
        tmp11 = tmp5 + tmp6;
        tmp12 = tmp6 + tmp7;

        /* The rotator is modified from fig 4-8 to avoid extra negations. */
        z5 = (tmp10 - tmp12) * 25080 >> 16;  /* c6 */
        z2 = (tmp10 * 35468 >> 16) + z5;     /* c2-c6 */
        z4 = (tmp12 * 85627 >> 16) + z5;     /* c2+c6 */
        z3 =  tmp11 * 46341 >> 16;           /* c4 */

        z11 = tmp7 + z3;  /* phase 5 */
        z13 = tmp7 - z3;

        dataptr[DCTSIZE * 5] = z13 + z2; /* phase 6 */
        dataptr[DCTSIZE * 3] = z13 - z2;
        dataptr[DCTSIZE * 1] = z11 + z4;
        dataptr[DCTSIZE * 7] = z11 - z4;

        dataptr++;  /* advance pointer to next column */
    }

    for (i=0; i<64; i++)
    {
        data[i]  *= enfactor[i];
        data[i] >>= 20;
    }
}

void idct2d8x8(long *data)
{
    long  tmp0,  tmp1,  tmp2,  tmp3;
    long  tmp4,  tmp5,  tmp6,  tmp7;
    long  tmp10, tmp11, tmp12, tmp13;
    long  z5, z10, z11, z12, z13;
    long *dataptr;
    int   ctr, i;

    for (i=0; i<64; i++)
    {
        data[i]  *= defactor[i];
        data[i] >>= 13;
    }

    /* Pass 1: process rows. */
    dataptr = data;
    for (ctr=0; ctr<DCTSIZE; ctr++)
    {
        if ( dataptr[1] + dataptr[2] + dataptr[3] + dataptr[4] +
             dataptr[5] + dataptr[6] + dataptr[7] == 0 )
        {
            dataptr[1] = dataptr[0];
            dataptr[2] = dataptr[0];
            dataptr[3] = dataptr[0];
            dataptr[4] = dataptr[0];
            dataptr[5] = dataptr[0];
            dataptr[6] = dataptr[0];
            dataptr[7] = dataptr[0];

            dataptr += DCTSIZE;
            continue;
        }

        /* Even part */
        tmp0 = dataptr[0];
        tmp1 = dataptr[2];
        tmp2 = dataptr[4];
        tmp3 = dataptr[6];

        tmp10 = tmp0 + tmp2;    /* phase 3 */
        tmp11 = tmp0 - tmp2;

        tmp13   = tmp1 + tmp3;   /* phases 5-3 */
        tmp12   = tmp1 - tmp3;   /* 2 * c4 */
        tmp12  *= 92682;
        tmp12 >>= 16;
        tmp12  -= tmp13;

        tmp0 = tmp10 + tmp13;   /* phase 2 */
        tmp3 = tmp10 - tmp13;
        tmp1 = tmp11 + tmp12;
        tmp2 = tmp11 - tmp12;
   
        /* Odd part */
        tmp4 = dataptr[1];
        tmp5 = dataptr[3];
        tmp6 = dataptr[5];
        tmp7 = dataptr[7];

        z13 = tmp6 + tmp5;    /* phase 6 */
        z10 = tmp6 - tmp5;
        z11 = tmp4 + tmp7;
        z12 = tmp4 - tmp7;

        tmp7    = z11 + z13;   /* phase 5 */
        tmp11   = z11 - z13;   /* 2 * c4 */
        tmp11  *= 92682;
        tmp11 >>= 16;

        z5 = (z10 + z12) * 121095 >> 16;  /*  2 * c2 */
        tmp10 =  (z12 *  70936 >> 16) - z5; /*  2 * (c2 - c6) */
        tmp12 = -(z10 * 171254 >> 16) + z5; /* -2 * (c2 + c6) */

        tmp6 = tmp12 - tmp7;    /* phase 2 */
        tmp5 = tmp11 - tmp6;
        tmp4 = tmp10 + tmp5;

        dataptr[0] = tmp0 + tmp7;
        dataptr[7] = tmp0 - tmp7;
        dataptr[1] = tmp1 + tmp6;
        dataptr[6] = tmp1 - tmp6;
        dataptr[2] = tmp2 + tmp5;
        dataptr[5] = tmp2 - tmp5;
        dataptr[4] = tmp3 + tmp4;
        dataptr[3] = tmp3 - tmp4;

        dataptr += DCTSIZE;
    }

    /* Pass 2: process columns. */
    dataptr = data;
    for (ctr=0; ctr<DCTSIZE; ctr++)
    {
        /* Even part */
        tmp0 = dataptr[DCTSIZE * 0];
        tmp1 = dataptr[DCTSIZE * 2];
        tmp2 = dataptr[DCTSIZE * 4];
        tmp3 = dataptr[DCTSIZE * 6];

        tmp10 = tmp0 + tmp2;    /* phase 3 */
        tmp11 = tmp0 - tmp2;

        tmp13   = tmp1 + tmp3;   /* phases 5-3 */
        tmp12   = tmp1 - tmp3;   /* 2 * c4 */
        tmp12  *= 92682;
        tmp12 >>= 16;
        tmp12  -= tmp13;

        tmp0 = tmp10 + tmp13;   /* phase 2 */
        tmp3 = tmp10 - tmp13;
        tmp1 = tmp11 + tmp12;
        tmp2 = tmp11 - tmp12;
   
        /* Odd part */
        tmp4 = dataptr[DCTSIZE * 1];
        tmp5 = dataptr[DCTSIZE * 3];
        tmp6 = dataptr[DCTSIZE * 5];
        tmp7 = dataptr[DCTSIZE * 7];

        z13 = tmp6 + tmp5;    /* phase 6 */
        z10 = tmp6 - tmp5;
        z11 = tmp4 + tmp7;
        z12 = tmp4 - tmp7;

        tmp7    = z11 + z13;   /* phase 5 */
        tmp11   = z11 - z13;   /* 2 * c4 */
        tmp11  *= 92682;
        tmp11 >>= 16;

        z5 = (z10 + z12) * 121095 >> 16;  /*  2 * c2 */
        tmp10 =  (z12 *  70936 >> 16) - z5; /*  2 * (c2 - c6) */
        tmp12 = -(z10 * 171254 >> 16) + z5; /* -2 * (c2 + c6) */

        tmp6 = tmp12 - tmp7;    /* phase 2 */
        tmp5 = tmp11 - tmp6;
        tmp4 = tmp10 + tmp5;

        dataptr[DCTSIZE * 0] = tmp0 + tmp7;
        dataptr[DCTSIZE * 7] = tmp0 - tmp7;
        dataptr[DCTSIZE * 1] = tmp1 + tmp6;
        dataptr[DCTSIZE * 6] = tmp1 - tmp6;
        dataptr[DCTSIZE * 2] = tmp2 + tmp5;
        dataptr[DCTSIZE * 5] = tmp2 - tmp5;
        dataptr[DCTSIZE * 4] = tmp3 + tmp4;
        dataptr[DCTSIZE * 3] = tmp3 - tmp4;

        dataptr++; /* advance pointers to next column */
    }

    for (i=0; i<64; i++) data[i] >>= 6;
}
#endif

#if 0 /* 快速的浮点数运算版本 */

/* 内部常量定义 */
#define DCTSIZE  8

/* 内部全局变量定义 */
static float aandctfactor[DCTSIZE] = {
    1.0f, 1.387039845f, 1.306562965f, 1.175875602f,
    1.0f, 0.785694958f, 0.541196100f, 0.275899379f,
};

static float dctfactor[64] = {0};

/* 内部函数实现 */
static initdctfactor()
{
    int i, j;
    for (i=0; i<DCTSIZE; i++)
        for (j=0; j<DCTSIZE; j++)
            dctfactor[i * 8 + j] = aandctfactor[i] * aandctfactor[j];
}

/* 函数实现 */
void fdct2d8x8(float *data)
{
    float tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
    float tmp10, tmp11, tmp12, tmp13;
    float z1, z2, z3, z4, z5, z11, z13;
    float *dataptr;
    int   ctr;

    /* Pass 1: process rows. */
    dataptr = data;
    for (ctr=0; ctr<DCTSIZE; ctr++)
    {
        tmp0 = dataptr[0] + dataptr[7];
        tmp7 = dataptr[0] - dataptr[7];
        tmp1 = dataptr[1] + dataptr[6];
        tmp6 = dataptr[1] - dataptr[6];
        tmp2 = dataptr[2] + dataptr[5];
        tmp5 = dataptr[2] - dataptr[5];
        tmp3 = dataptr[3] + dataptr[4];
        tmp4 = dataptr[3] - dataptr[4];
   
        /* Even part */
        tmp10 = tmp0 + tmp3;  /* phase 2 */
        tmp13 = tmp0 - tmp3;
        tmp11 = tmp1 + tmp2;
        tmp12 = tmp1 - tmp2;

        dataptr[0] = tmp10 + tmp11;  /* phase 3 */
        dataptr[4] = tmp10 - tmp11;

        z1 = (tmp12 + tmp13) * 0.707106781f; /* c4 */
        dataptr[2] = tmp13 + z1;  /* phase 5 */
        dataptr[6] = tmp13 - z1;

        /* Odd part */
        tmp10 = tmp4 + tmp5;  /* phase 2 */
        tmp11 = tmp5 + tmp6;
        tmp12 = tmp6 + tmp7;

        /* The rotator is modified from fig 4-8 to avoid extra negations. */
        z5 = (tmp10 - tmp12) * 0.382683433f; /* c6 */
        z2 = 0.541196100f * tmp10 + z5;      /* c2-c6 */
        z4 = 1.306562965f * tmp12 + z5;      /* c2+c6 */
        z3 = tmp11 * 0.707106781f;           /* c4 */

        z11 = tmp7 + z3;        /* phase 5 */
        z13 = tmp7 - z3;

        dataptr[5] = z13 + z2;  /* phase 6 */
        dataptr[3] = z13 - z2;
        dataptr[1] = z11 + z4;
        dataptr[7] = z11 - z4;

        dataptr += DCTSIZE;     /* advance pointer to next row */
    }

    /* Pass 2: process columns. */
    dataptr = data;
    for (ctr=0; ctr<DCTSIZE; ctr++)
    {
        tmp0 = dataptr[DCTSIZE * 0] + dataptr[DCTSIZE * 7];
        tmp7 = dataptr[DCTSIZE * 0] - dataptr[DCTSIZE * 7];
        tmp1 = dataptr[DCTSIZE * 1] + dataptr[DCTSIZE * 6];
        tmp6 = dataptr[DCTSIZE * 1] - dataptr[DCTSIZE * 6];
        tmp2 = dataptr[DCTSIZE * 2] + dataptr[DCTSIZE * 5];
        tmp5 = dataptr[DCTSIZE * 2] - dataptr[DCTSIZE * 5];
        tmp3 = dataptr[DCTSIZE * 3] + dataptr[DCTSIZE * 4];
        tmp4 = dataptr[DCTSIZE * 3] - dataptr[DCTSIZE * 4];

        /* Even part */
        tmp10 = tmp0 + tmp3;  /* phase 2 */
        tmp13 = tmp0 - tmp3;
        tmp11 = tmp1 + tmp2;
        tmp12 = tmp1 - tmp2;

        dataptr[DCTSIZE * 0] = tmp10 + tmp11;  /* phase 3 */
        dataptr[DCTSIZE * 4] = tmp10 - tmp11;

        z1 = (tmp12 + tmp13) * 0.707106781f; /* c4 */
        dataptr[DCTSIZE * 2] = tmp13 + z1;  /* phase 5 */
        dataptr[DCTSIZE * 6] = tmp13 - z1;

        /* Odd part */
        tmp10 = tmp4 + tmp5;  /* phase 2 */
        tmp11 = tmp5 + tmp6;
        tmp12 = tmp6 + tmp7;

        /* The rotator is modified from fig 4-8 to avoid extra negations. */
        z5 = (tmp10 - tmp12) * 0.382683433f; /* c6 */
        z2 = 0.541196100f * tmp10 + z5;      /* c2-c6 */
        z4 = 1.306562965f * tmp12 + z5;      /* c2+c6 */
        z3 = tmp11 * 0.707106781f;           /* c4 */

        z11 = tmp7 + z3;  /* phase 5 */
        z13 = tmp7 - z3;

        dataptr[DCTSIZE * 5] = z13 + z2; /* phase 6 */
        dataptr[DCTSIZE * 3] = z13 - z2;
        dataptr[DCTSIZE * 1] = z11 + z4;
        dataptr[DCTSIZE * 7] = z11 - z4;

        dataptr++;  /* advance pointer to next column */
    }
}

void idct2d8x8(float *data)
{
    float  tmp0,  tmp1,  tmp2,  tmp3;
    float  tmp4,  tmp5,  tmp6,  tmp7;
    float  tmp10, tmp11, tmp12, tmp13;
    float  z5, z10, z11, z12, z13;
    float *dataptr;
    int    ctr;

    /* Pass 1: process rows. */
    dataptr = data;
    for (ctr=0; ctr<DCTSIZE; ctr++)
    {
        /* Even part */
        tmp0 = dataptr[0];
        tmp1 = dataptr[2];
        tmp2 = dataptr[4];
        tmp3 = dataptr[6];

        tmp10 = tmp0 + tmp2;    /* phase 3 */
        tmp11 = tmp0 - tmp2;

        tmp13  = tmp1 + tmp3;   /* phases 5-3 */
        tmp12  = tmp1 - tmp3;   /* 2 * c4 */
        tmp12 *= 1.414213562f;
        tmp12 -= tmp13;

        tmp0 = tmp10 + tmp13;   /* phase 2 */
        tmp3 = tmp10 - tmp13;
        tmp1 = tmp11 + tmp12;
        tmp2 = tmp11 - tmp12;
   
        /* Odd part */
        tmp4 = dataptr[1];
        tmp5 = dataptr[3];
        tmp6 = dataptr[5];
        tmp7 = dataptr[7];

        z13 = tmp6 + tmp5;    /* phase 6 */
        z10 = tmp6 - tmp5;
        z11 = tmp4 + tmp7;
        z12 = tmp4 - tmp7;

        tmp7   = z11 + z13;   /* phase 5 */
        tmp11  = z11 - z13;   /* 2 * c4 */
        tmp11 *= 1.414213562f;

        z5 = (z10 + z12) * 1.847759065f;  /*  2 * c2 */
        tmp10 =  1.082392200f * z12 - z5; /*  2 * (c2 - c6) */
        tmp12 = -2.613125930f * z10 + z5; /* -2 * (c2 + c6) */

        tmp6 = tmp12 - tmp7;    /* phase 2 */
        tmp5 = tmp11 - tmp6;
        tmp4 = tmp10 + tmp5;

        dataptr[0] = tmp0 + tmp7;
        dataptr[7] = tmp0 - tmp7;
        dataptr[1] = tmp1 + tmp6;
        dataptr[6] = tmp1 - tmp6;
        dataptr[2] = tmp2 + tmp5;
        dataptr[5] = tmp2 - tmp5;
        dataptr[4] = tmp3 + tmp4;
        dataptr[3] = tmp3 - tmp4;

        dataptr += DCTSIZE;
    }

    /* Pass 2: process columns. */
    dataptr = data;
    for (ctr=0; ctr<DCTSIZE; ctr++)
    {
        /* Even part */
        tmp0 = dataptr[DCTSIZE * 0];
        tmp1 = dataptr[DCTSIZE * 2];
        tmp2 = dataptr[DCTSIZE * 4];
        tmp3 = dataptr[DCTSIZE * 6];

        tmp10 = tmp0 + tmp2;    /* phase 3 */
        tmp11 = tmp0 - tmp2;

        tmp13  = tmp1 + tmp3;   /* phases 5-3 */
        tmp12  = tmp1 - tmp3;   /* 2 * c4 */
        tmp12 *= 1.414213562f;
        tmp12 -= tmp13;

        tmp0 = tmp10 + tmp13;   /* phase 2 */
        tmp3 = tmp10 - tmp13;
        tmp1 = tmp11 + tmp12;
        tmp2 = tmp11 - tmp12;
   
        /* Odd part */
        tmp4 = dataptr[DCTSIZE * 1];
        tmp5 = dataptr[DCTSIZE * 3];
        tmp6 = dataptr[DCTSIZE * 5];
        tmp7 = dataptr[DCTSIZE * 7];

        z13 = tmp6 + tmp5;    /* phase 6 */
        z10 = tmp6 - tmp5;
        z11 = tmp4 + tmp7;
        z12 = tmp4 - tmp7;

        tmp7   = z11 + z13;   /* phase 5 */
        tmp11  = z11 - z13;   /* 2 * c4 */
        tmp11 *= 1.414213562f;

        z5 = (z10 + z12) * 1.847759065f;  /*  2 * c2 */
        tmp10 =  1.082392200f * z12 - z5; /*  2 * (c2 - c6) */
        tmp12 = -2.613125930f * z10 + z5; /* -2 * (c2 + c6) */

        tmp6 = tmp12 - tmp7;    /* phase 2 */
        tmp5 = tmp11 - tmp6;
        tmp4 = tmp10 + tmp5;

        dataptr[DCTSIZE * 0] = tmp0 + tmp7;
        dataptr[DCTSIZE * 7] = tmp0 - tmp7;
        dataptr[DCTSIZE * 1] = tmp1 + tmp6;
        dataptr[DCTSIZE * 6] = tmp1 - tmp6;
        dataptr[DCTSIZE * 2] = tmp2 + tmp5;
        dataptr[DCTSIZE * 5] = tmp2 - tmp5;
        dataptr[DCTSIZE * 4] = tmp3 + tmp4;
        dataptr[DCTSIZE * 3] = tmp3 - tmp4;

        dataptr++; /* advance pointers to next column */
    }
}

/* how to use */
static void dcttest(void)
{
    float data[64];
    int   i;

    /* FDCT */
    fdct2d8x8(data);
    for (i=0; i<64; i++) data[i] /= dctfactor[i];
    for (i=0; i<64; i++) data[i] /= 8.0f;

    /* IDCT */
    for (i=0; i<64; i++) data[i] *= dctfactor[i];
    idct2d8x8(data);
    for (i=0; i<64; i++) data[i] /= 8.0f;
}

#endif

#if 0 /* 矩阵变换版本 */
/* 内部常量定义 */
#define M_PI  3.14159265358979323846f

/* 内部全局变量定义 */
static float fdctmatrix[8][8] = {0};  /* fdct 变换矩阵 */
static float idctmatrix[8][8] = {0};  /* idct 变换矩阵 */

/* 内部函数定义 */
static float c(int u)
{
    if (u == 0) return (float)sqrt(1.0f / 8.0f);
    else        return (float)sqrt(2.0f / 8.0f);
}

/* 初始化 dct 变换矩阵 */
void initdctmatrix(void)
{
    static int inited = 0;
    int    u, x;

    /* 避免重复初始化 */
    if (inited) return;

    /* init fdct matrix */
    for (u=0; u<8; u++)
    {
        for (x=0; x<8; x++)
        {
            fdctmatrix[u][x] = (float)(c(u) * cos((2.0f * x + 1.0f) * u * M_PI / 16.0f));
        }
    }

    /* init idct matrix */
    for (u=0; u<8; u++)
    {
        for (x=0; x<8; x++)
        {
            idctmatrix[x][u] = (float)(c(u) * cos((2.0f * x + 1.0f) * u * M_PI / 16.0f));
        }
    }

    inited = 1;
}

/* 1d 8-points forward dct */
void fdct1d8(float *u, float *x)
{
    int i;
    int j;

    for (i=0; i<8; i++)
    {
        u[i] = 0.0f;
        for (j=0; j<8; j++)
        {
            u[i] += fdctmatrix[i][j] * x[j];
        }
    }
}

/* 1d 8-points invert dct */
void idct1d8(float *x, float *u)
{
    int i;
    int j;

    for (i=0; i<8; i++)
    {
        x[i] = 0.0f;
        for (j=0; j<8; j++)
        {
            x[i] += idctmatrix[i][j] * u[j];
        }
    }
}

/* 2d 8x8 forward dct */
void fdct2d8x8(float *data)
{
    float temp[64];
    int   i, j;

    /* 初始化变换矩阵 */
    initdctmatrix();

    /* 逐行进行 1d fdct */
    for (i=0; i<8; i++)
    {
        fdct1d8(temp + 8 * i, data + 8 * i);
    }

    /* 转置矩阵 */
    for (i=0; i<8; i++)
        for (j=0; j<8; j++)
            *(data + 8 * i + j) = *(temp + 8 * j + i);

    /* 逐行进行 1d fdct */
    for (i=0; i<8; i++)
    {
        fdct1d8(temp + 8 * i, data + 8 * i);
    }

    /* 转置矩阵 */
    for (i=0; i<8; i++)
        for (j=0; j<8; j++)
            *(data + 8 * i + j) = *(temp + 8 * j + i);
}

/* 2d 8x8 invert dct */
void idct2d8x8(float *data)
{
    float temp[64];
    int   i, j;

    /* 初始化变换矩阵 */
    initdctmatrix();

    /* 逐行进行 1d idct */
    for (i=0; i<8; i++)
    {
        idct1d8(temp + 8 * i, data + 8 * i);
    }

    /* 转置矩阵 */
    for (i=0; i<8; i++)
        for (j=0; j<8; j++)
            *(data + 8 * i + j) = *(temp + 8 * j + i);

    /* 逐行进行 1d idct */
    for (i=0; i<8; i++)
    {
        idct1d8(temp + 8 * i, data + 8 * i);
    }

    /* 转置矩阵 */
    for (i=0; i<8; i++)
        for (j=0; j<8; j++)
            *(data + 8 * i + j) = *(temp + 8 * j + i);
}
#endif

#if 0 /* 数学表达式版本 */
/* 内部常量定义 */
#define M_PI  3.14159265358979323846f

/* 函数实现 */
static float alpha(int n)
{
    if (n == 0) return 1.0f / (float)sqrt(8);
    else        return 1.0f / 2.0f;
}

void fdct2d8x8(float *data)
{
    int u, v;
    int x, y, i;
    float buf[64];
    float temp;

    for (u=0; u<8; u++)
    {
        for (v=0; v<8; v++)
        {
            temp = 0;
            for (x=0; x<8; x++)
            {
                for (y=0; y<8; y++)
                {
                    temp += data[y * 8 + x]
                          * (float)cos((2.0f * x + 1.0f) / 16.0f * u * M_PI)
                          * (float)cos((2.0f * y + 1.0f) / 16.0f * v * M_PI);
                }
            }
            buf[v * 8 + u] = alpha(u) * alpha(v) * temp;
        }
    }

    for (i=0; i<64; i++) data[i] = buf[i];
}

void idct2d8x8(float *data)
{
    int u, v;
    int x, y, i;
    float buf[64];
    float temp;

    for (x=0; x<8; x++)
    {
        for (y=0; y<8; y++)
        {
            temp = 0;
            for (u=0; u<8; u++)
            {
                for (v=0; v<8; v++)
                {
                    temp += alpha(u) * alpha(v) * (data[v * 8 + u]
                          * (float)cos((2.0f * x + 1.0f) / 16.0f * u * M_PI)
                          * (float)cos((2.0f * y + 1.0f) / 16.0f * v * M_PI));
                    
                }
            }
            buf[y * 8 + x] = (char)(temp + 0.5);
        }
    }

    for (i=0; i<64; i++) data[i] = buf[i];
}
#endif

值得一提的是，libjpeg 中的快速 DCT 变换算法并不是完整的算法，变换的结果还需要与一个系数矩阵相乘才行，具体看上面的代码。这个问题也是困扰了我很久，我很早就想采用这个算法，但是从 libjpeg 中抽出来之后，发现变换的结果不对，仔细分析 libjpeg 的代码，才发现原因。

dct不用这么复杂吧？能说说您的算法哪里比较hao?

大家可以从最原始的版本，即表达式版本开始，到矩阵运算版本，再到基于 AA&N 的快速算法，以及从浮点到整型运算的改进，逐步看出算法的优化。

我之所以要在这里总结出来，就是因为我曾经搜遍了互联网，都没有找到一个像样的代码，许多网站的文章都写得不明不白，而要想找到一个真正可用的，清晰可读的代码，几乎就是不可能，稍微好一点的论文，TMD 的还要叫你注册会员，然后交钱，千方百计弄到手来，一看就是个垃圾，没有什么参考价值，这就是当代本科生和硕士生的论文，呵呵，读过大学的都知道是怎么回事儿。

我主要参考的文档如下：

另外就是范式哈夫曼编解码的相关文档了

我这里有你以前写的huffman程式，都不敢看，真是复杂
范式哈夫曼编解码，原理很简单。但和JPEG结合起来，就是想不明白。
你要是能写篇文章讲讲就好了。

最新的REG里面，发现少了不少东西

由于前段时间电脑被盗了，很多东西都丢了，我发布的 RGE v0.2 也是一个保留下来的早期版本。新买的电脑在 XP 下对 dos 的支持都不好，所以 RGE 的开发暂时停止了。

哥们，你的DCT变换效率不行。

我的是采用的 AN&N 的 DCT 算法，并且将浮点运算转换为了定点运算，用到的也就是加法和移位。如果不考虑采用使用处理器专门提供的 SIMD 指令（例如 SSE），这应该是 C 语言所能写出的比较快的算法了。其实这个也不是我的功劳，我也是参考了 libjpeg 才写出来的。x264 里面的代码我没看过，但是就从算法来讲应该都是差不多的，只是在具体优化上的不同了，比如结合特定的处理器架构，采用专门的汇编指令进行优化。

顶斑竹

这是我的转换代码。。。。请斑竹指点下。。。
long B2J()
{
  if(type==24)
  {
    JPG a;
    long i,j,c,k;
    long m,n,tp;
    long u,v,x,y;//用于二维DCT变换
    float cv,cu,pi;
    unsigned char count;
    double sum[6]={0};
    unsigned char R,G,B;
    signed char *Y,*U,*V;
    signed char MCU[6][64]={0};
    signed char Q_DCT[6][64]={0};
    float DCT[6][64]={0};
    short sgn;
    signed char Z[64]={0};
    signed char transt[317]={0};
    signed char PreDC[3]={0};
    //int location=0,run=0;
    unsigned char run=0;
    int  r_tp=0;
    long r_len=0;


    bit *buff;
    //unsigned char y,u,v;

    unsigned char ZigZag[64]={0,1,5,6,14,15,27,28,
                         2,4,7,13,16,26,29,42,
                         3,8,12,17,25,30,41,43,
                         9,11,18,24,31,40,44,53,
                         10,19,23,32,39,45,52,54,
                         20,22,33,38,46,51,55,60,
                         21,34,37,47,50,56,59,61,
                         35,36,48,49,57,58,62,63
                        };

    pi=(float)11.25;
    c=h;
    k=w;
    while(c%16!=0) c++;
    while(k%16!=0) k++;

    Y=new signed char[c*k];//获取BMP主要图像信息
    U=new signed char[c*k];
    V=new signed char[c*k];
    memset(Y,0,c*k);
    memset(U,0,c*k);
    memset(V,0,c*k);

    p_temp=new unsigned char[c*k*3>>1];//4:1:1比例最大为一半，记录最后的编码
    memset(p_temp,0,c*k*3>>1);

    buff=new bit[c*k*12];
    for(i=0;i<c*k*12;i++)
        buff[i].a =1;
    //YUV
    for(i=0;i<h;i++)
    {
        for(j=0;j<w;j++)
        {
            R=p_data[i*wf+3*j+2];//Y
            G=p_data[i*wf+3*j+1];//U
            B=p_data[i*wf+3*j];//V
            Y[i*w+j]=0.2990*R+0.5870*G+0.1140*B;
            U[i*w+j]=-0.1687*R-0.3313*G+0.5000*B;
            V[i*w+j]=0.5000*R-0.4187*G-0.0813*B;   
        }
    }
    //4:1:1 sample
    for(i=0;i<c>>4;i++)
        for(j=0;j<k>>4;j++)
        {
            for(n=i*8;n<i*8+8;n++)
                for(m=j*8;m<j*8+8;m++)
                {
                    tp=(n-i*8)*8+(m-j*8);
                    MCU[0][tp]=Y[n*k+m]-128;
                    MCU[1][tp]=Y[n*k+m+8]-128;
                    MCU[2][tp]=Y[(n+8)*k+m]-128;
                    MCU[3][tp]=Y[(n+8)*k+m+8]-128;
                    MCU[4][tp]=U[2*n*k+2*m];
                    MCU[5][tp]=V[2*n*k+2*m];
                }

            /*   
                printf("************************MCU\n*****************");
                for(count=0;count<6;count++)
                {
                    printf("MCU编号%d\n",count);
                    for(y=0;y<64;y++)
                        {
                            if(1) printf("%d ",MCU[count][y]);
                            if((y+1)%8==0) printf("\n");

                    }

                }
                            getch();
                            printf("\n");
            */
          //DCT
                for(v=0;v<8;v++)
                {
                    if(v==0) cv=1/sqrt(2);
                    else     cv=1;

                    for(u=0;u<8;u++)
                    {
                        if(u==0)    cu=1/sqrt(2);
                        else    cu=1;

                        for(count=0;count<6;count++)
                            sum[count]=0;

                        for(y=0;y<8;y++)
                            for(x=0;x<8;x++)
                            {
                                for(count=0;count<6;count++)
                                    sum[count]=sum[count]+MCU[count][y*8+x]*cos((2*x+1)*u*pi)*cos((2*y+1)*v*pi);
                                    
                            }

                            for(count=0;count<6;count++)
                            {
                                DCT[count][v*8+u]=sum[count]*cu*cv/4;
                            }                                
                    }
                }



            //Q
                for(count=0;count<6;count++)
                {
                    for(y=0;y<64;y++)
                        {
                            if(DCT[count][y]>0) sgn=1; else sgn=-1;
                            if(count<4)    Q_DCT[count][y]=(DCT[count][y]+sgn*a.dqt0.QT[y]/2)/a.dqt0.QT[y];
                            if(count>3)    Q_DCT[count][y]=(DCT[count][y]+sgn*a.dqt1.QT[y]/2)/a.dqt1.QT[y];
                            //if(1) printf("%d ",Q_DCT[count][y]);
                            //if((y+1)%8==0) printf("\n");
                            //getch();
                            //printf("\n");
                        }
                    //printf("************\n");
                }
                //ZigZag
                for(count=0;count<6;count++)
                {


                    for(y=0;y<64;y++)
                    {
                        Z[ZigZag[y]]=Q_DCT[count][y];
                    }
                    //Z字编码显示
                    /*
                    printf("原Q块\n");
                    for(y=0;y<64;y++)
                    {
                        printf("%5d",Q_DCT[count][y]);
                        if((y+1)%8==0) printf("\n");
                    }
                    printf("编码块\n");
                    for(y=0;y<64;y++)
                    {
                        printf("%5d",Z[y]);
                        if((y+1)%8==0) printf("\n");
                    }
                    getch();*/
                    memcpy(Q_DCT[count],Z,64);
                }

            //HuffmanEncoder
            for(count=0;count<6;count++)
            {
                run=0;//记录一个8*8像素块的过渡段行程
                //过渡DC
                if(count<4)
                {
                    transt[run+1]=Q_DCT[count][0]-PreDC[0];
                    x=0;
                    while(abs(transt[run+1])>>x!=0) x++;
                    transt[run]=x;
                    PreDC[0]=Q_DCT[count][0];//记录当前DC为下一次差分做准备
                }
                if(count==4)
                {
                    transt[run+1]=Q_DCT[count][0]-PreDC[1];
                    x=0;
                    while(abs(transt[run+1])>>x!=0) x++;
                    transt[run]=x;
                    PreDC[1]=Q_DCT[count][0];//记录当前DC为下一次差分做准备               
                }
                if(count==5)
                {
                    transt[run+1]=Q_DCT[count][0]-PreDC[2];
                    x=0;
                    while(abs(transt[run+1])>>x!=0) x++;
                    transt[run]=x;
                    PreDC[2]=Q_DCT[count][0];//记录当前DC为下一次差分做准备               
                }
                run=run+2;
               
                x=0;
                //过渡AC
                for(y=1;y<64;y++)
                {
                    if(Q_DCT[count][y]==0)
                    {
                        x++;
                        if(x>15)//x=16
                        {
                            r_tp=y+1;
                            while(Q_DCT[count][r_tp]==0 && r_tp<64)    r_tp++;
                            if(r_tp==64)
                            {
                                transt[run]=transt[run+1]=0;//00标记结尾
                                run=run+2;
                                y=64;
                                goto transt_end;//标记结束编码
                            }
                            else
                            {
                                transt[run]=15;
                                transt[run+1]=0;
                                transt[run+2]=0;
                                run=run+3;
                                x=0;//重新计0
                            }
                           
                        }
                    }
                    else
                    {
                        transt[run]=x;//记录0个数
                        x=0;
                        while(abs(Q_DCT[count][y])>>x!=0) x++;
                        transt[run+1]=x;//记录区间
                        transt[run+2]=Q_DCT[count][y];
                        x=0;//记录下一个行程
                        run=run+3;
                    }               

                }
                    if(Q_DCT[count][63]==0)
                    {
                        transt[run]=transt[run+1]=0;//00标记结尾
                        run=run+2;
                    }
transt_end:;
                /*
                //过渡块显示输出
                for(y=0;y<64;y++)
                {
                            if(1) printf("%d ",Q_DCT[count][y]);
                            if((y+1)%8==0) printf("\n");                    
                }
                            if(count<4) printf("亮度块");
                            if(count==4)   printf("色度块");
                            if(count==5)   printf("深度块");
                            printf("****************************\n");

                for(y=0;y<run;y++)
                if(1) printf("%d ",transt[y]);
                printf("\n对应中间编码,共%d个\n",run);
                getch();
               
                */
               
                //过渡段编码
                //p_temp,buff,transt,x,y,r_tp,r_len

                unsigned char b_cmp=0;
                unsigned short *DCtable,*ACtable;
                unsigned char *DCL;//,*ACL;
               
                //unsigned char *DCtable,*ACtable;
                if(count<4)
                {
                    DCtable=a.DCTY.C ;
                    DCL=a.DCTY.L;
                    ACtable=a.ACTY.C ;
                    //ACL=a.dhtac0.sum;
                }
                else
                {
                    DCtable=a.DCTC.C ;
                    DCL=a.DCTC.L;
                    ACtable=a.ACTC.C ;
                    //ACL=a.dhtac0.sum;
                }
                //DC
                //r_len=0;
                x=DCL[transt[0]];
                while(x!=0)
                {
                    x--;
                    buff[r_len++].a=DCtable[transt[0]]>>x;
                }
               
                x=transt[0];
                if(transt[1]<0)  b_cmp=((unsigned char)0x01 <<x)-1-abs(transt[1]);
                else b_cmp=transt[1];
                while(x!=0)
                {
                    x--;
                    buff[r_len++].a=b_cmp>>x;
                }

                //AC
                for(y=2;y<run-2;y=y+3)
                {
                    x=0;
                    r_tp=transt[y]*16+transt[y+1];//RRRRSSSS
                    

                    if(r_tp==1 || (r_tp==0 && count>3) ||(r_tp==2 && count<4)) x=2;
                    else
                    while(ACtable[r_tp]>>x!=0) x++;
                    
                    //x=ACL[r_tp];

                    while(x!=0)
                    {
                        x--;
                        buff[r_len++].a=ACtable[r_tp]>>x;
                    }   
                    
                    x=transt[y+1];
                    if(transt[y+2]<0)  b_cmp=((unsigned char)0x01 <<x)-1-abs(transt[y+2]);
                    else b_cmp=transt[y+2];
                    while(x!=0)
                    {
                        x--;
                        buff[r_len++].a=b_cmp>>x;
                    }
                }
                //结尾
                    x=0;
                    r_tp=0;//RRRRSSSS
                    

                    if(count>3)
                        x=2;
                    else
                        x=4;

                    while(x!=0)
                    {
                        x--;
                        buff[r_len++].a=ACtable[r_tp]>>x;
                    }

            }               
        }
        
        //开始转换字节
        r_tp=0;
        while(r_len%8==0)
        {
            r_len++;
            r_tp=1;
        }
        y=0;
        for(x=0;x<r_len;x=x+8)
        {
            p_temp[y]=((unsigned char)buff[x].a<<7)+((unsigned char)buff[x+1].a<<6)+((unsigned char)buff[x+2].a<<5)+((unsigned char)buff[x+3].a<<4)+((unsigned char)buff[x+4].a<<3)+((unsigned char)buff[x+5].a<<2)+((unsigned char)buff[x+6].a<<1)+(unsigned char)buff[x+7].a;
            
            if(p_temp[y]==0xFF)
            {
                p_temp[y+1]=0x00;
                y=y+2;
            }
            else
                y=y+1;
        }
        if(r_tp==1) //如果有位填充需要补FF
        {
            p_temp[y++]=0xFF;
            p_temp[y++]=0x00;
        }
            printf("\n转换长度%ld",y);
    //delete []p_temp;p_temp=0;//由于不是正常处理图像，所以需要清空指针
    delete []buff;
    return y;
    //exit(0);
  }
  else
          return -1;
}