快速离散余弦正逆变换算法
网上常见的快速离散余弦变换的代码如下：
#include
using namespace std;
#define DCTSIZE 8
void FDCT(double* lpBuff)
{
double tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
double tmp10, tmp11, tmp12, tmp13;
double z1, z2, z3, z4, z5, z11, z13;
double* dataptr;
int ctr;
/* 第一部分，对行进行计算 */
dataptr = lpBuff;
for (ctr = DCTSIZE-1; ctr >= 0; 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];
/* 对偶数项进行运算 */
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.707106781); /* c4 */
dataptr[2] = tmp13 + z1; /* phase 5 */
dataptr[6] = tmp13 - z1;
/* 对奇数项进行计算 */
tmp10 = tmp4 + tmp5; /* phase 2 */
tmp11 = tmp5 + tmp6;
tmp12 = tmp6 + tmp7;
z5 = (tmp10 - tmp12) * ( 0.382683433); /* c6 */
z2 = (0.541196100) * tmp10 + z5; /* c2-c6 */
z4 = (1.306562965) * tmp12 + z5; /* c2+c6 */
z3 = tmp11 * (0.707106781); /* 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; /* 将指针指向下一行 */
}
/* 第二部分，对列进行计算 */
dataptr = lpBuff;
for (ctr = DCTSIZE-1; ctr >= 0; 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];
/* 对偶数项进行运算 */
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.707106781); /* c4 */
dataptr[DCTSIZE*2] = tmp13 + z1; /* phase 5 */
dataptr[DCTSIZE*6] = tmp13 - z1;
/* 对奇数项进行计算 */
tmp10 = tmp4 + tmp5; /* phase 2 */
tmp11 = tmp5 + tmp6;
tmp12 = tmp6 + tmp7;
z5 = (tmp10 - tmp12) * (0.382683433); /* c6 */
z2 = (0.541196100) * tmp10 + z5; /* c2-c6 */
z4 = (1.306562965) * tmp12 + z5; /* c2+c6 */
z3 = tmp11 * (0.707106781); /* c4 */
z11 = tmp7 + z3; /* phase 5 */
z13 = tmp7 - z3;
dataptr[DCTSIZE*5] = z13 + z2; /* p
hase 6 */
dataptr[DCTSIZE*3] = z13 - z2;
dataptr[DCTSIZE*1] = z11 + z4;
dataptr[DCTSIZE*7] = z11 - z4;
++dataptr; /* 将指针指向下一列 */
}
}
int main(void)
{
const short NUM = 8;
short i = 0;
short j = 0;
short u = 0;
short v = 0;
double input[NUM][NUM] =
{
{89.00, 101.00, 114.00, 125.00, 126.00, 115.00, 105.00, 96.00},
{97.00, 115.00, 131.00, 147.00, 149.00, 135.00, 123.00, 113.00},
{114.00, 134.00, 159.00, 178.00, 175.00, 164.00, 149.00, 137.00},
{121.00, 143.00, 177.00, 196.00, 201.00, 189.00, 165.00, 150.00},
{119.00, 141.00, 175.00, 201.00, 207.00, 186.00, 162.00, 144.00},
{107.00, 130.00, 165.00, 189.00, 192.00, 171.00, 144.00, 125.00},
{97.00, 119.00, 149.00, 171.00, 172.00, 145.00, 117.00, 96.00},
{88.00, 107.00, 136.00, 156.00, 155.00, 129.00, 97.00, 75.00}
};
double *data = new double[NUM * NUM];
for(i = 0; i < NUM; i++)
{
for(j = 0; j < NUM; j++)
{
data[i * NUM + j] = input[i][j];
}
}
FDCT(data);
// print the result of FDCT
for(u = 0; u < NUM; u++)
{
for(v = 0; v < NUM; v++)
{
printf("%7.2f ", data[u*NUM + v]);
}
cout << endl;
}
return 0;
}
运算得到的结果：
9000.00 -353.39 -1933.67 -61.84 16.00 -4.74 -10.33 3.96
-241.34 -243.40 650.60 -42.79 -24.64 -1.12 -9.34 -5.43
-1725.67 464.18 230.84 25.81 13.49 -6.50 -15.02 -0.25
-66.65 -56.94 2.18 17.28 15.93 -5.71 -5.37 4.21
-18.00 2.53 -5.07 28.89 -2.00 -1.57 9.07 2.16
20.23 8.66 11.55 -5.36 -11.20 2.24 5.68 0.07
-0.33 -1.85 9.02 -6.61 -3.49 6.67 3.16 -1.46
-0.24 8.50 3.08 -2.01 3.92 1.92 -2.39 -0.12
这个结果和离散余弦正逆变换中最原始的或者说是真正的离散余弦变换的结果是不同的。原因请看下面程序中的红色字体部分。
=================================================================
// FDCT.h
void fdct(double *block);
// FDCT.cpp
#define USE_ACCURATE_ROUNDING
#define RIGHT_SHIFT(x, shft) ((x) >> (shft))
#ifdef USE_ACCURATE_ROUNDING
#define ONE ((int) 1)
#define DESCALE(x, n) RIGHT_SHIFT((x) + (ONE << ((n) - 1)), n)
#else
#define DESCALE(x, n) RIGHT_SHIFT(x, n)
#endif
#define CONST_BITS 13
#define PASS1_BITS 2
#define FIX_0_298631336 ((int) 2446) /* FIX(0.298631336) */
#define FIX_0_390180644 ((int) 3196) /* FIX(0.390180644) */
#define FIX_0_541196100 ((int) 4433) /* FIX(0.541196100) */
#define FIX_0_765366865 ((int) 6270) /* FIX(0.765366865) */
#define FIX_0_899976223 ((int) 7373) /* FIX(0.899976223) */
#define FIX_1_175875602 ((int) 9633) /* FIX(1.175875602) */
#define FIX_1_501321110 ((int) 12299) /* FIX(1.501321110) */
#define FIX_1_847759065 ((int) 15137) /* FIX(1.847759065) */
#def
ine FIX_1_961570560 ((int) 16069) /* FIX(1.961570560) */
#define FIX_2_053119869 ((int) 16819) /* FIX(2.053119869) */
#define FIX_2_562915447 ((int) 20995) /* FIX(2.562915447) */
#define FIX_3_072711026 ((int) 25172) /* FIX(3.072711026) */
void fdct(double *block)
{
int tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7;
int tmp10, tmp11, tmp12, tmp13;
int z1, z2, z3, z4, z5;
double *blkptr;
int *dataptr;
int data[64];
int i;
/* Pass 1: process rows. */
/* Note results are scaled up by sqrt(8) compared to a true DCT; */
/* furthermore, we scale the results by 2**PASS1_BITS. */
dataptr = data;
blkptr = block;
for (i = 0; i < 8; i++) {
tmp0 = (int)(blkptr[0] + blkptr[7]);
tmp7 = (int)(blkptr[0] - blkptr[7]);
tmp1 = (int)(blkptr[1] + blkptr[6]);
tmp6 = (int)(blkptr[1] - blkptr[6]);
tmp2 = (int)(blkptr[2] + blkptr[5]);
tmp5 = (int)(blkptr[2] - blkptr[5]);
tmp3 = (int)(blkptr[3] + blkptr[4]);
tmp4 = (int)(blkptr[3] - blkptr[4]);
tmp10 = tmp0 + tmp3;
tmp13 = tmp0 - tmp3;
tmp11 = tmp1 + tmp2;
tmp12 = tmp1 - tmp2;
dataptr[0] = (tmp10 + tmp11) << PASS1_BITS;
dataptr[4] = (tmp10 - tmp11) << PASS1_BITS;
z1 = (tmp12 + tmp13) * FIX_0_541196100;
dataptr[2] = DESCALE(z1 + tmp13 * FIX_0_765366865, CONST_BITS - PASS1_BITS);
dataptr[6] = DESCALE(z1 + tmp12 * (-FIX_1_847759065), CONST_BITS - PASS1_BITS);
z1 = tmp4 + tmp7;
z2 = tmp5 + tmp6;
z3 = tmp4 + tmp6;
z4 = tmp5 + tmp7;
z5 = (z3 + z4) * FIX_1_175875602; /* sqrt(2) * c3 */
tmp4 *= FIX_0_298631336; /* sqrt(2) * (-c1+c3+c5-c7) */
tmp5 *= FIX_2_053119869; /* sqrt(2) * ( c1+c3-c5+c7) */
tmp6 *= FIX_3_072711026; /* sqrt(2) * ( c1+c3+c5-c7) */
tmp7 *= FIX_1_501321110; /* sqrt(2) * ( c1+c3-c5-c7) */
z1 *= -FIX_0_899976223; /* sqrt(2) * (c7-c3) */
z2 *= -FIX_2_562915447; /* sqrt(2) * (-c1-c3) */
z3 *= -FIX_1_961570560; /* sqrt(2) * (-c3-c5) */
z4 *= -FIX_0_390180644; /* sqrt(2) * (c5-c3) */
z3 += z5;
z4 += z5;
dataptr[7] = DESCALE(tmp4 + z1 + z3, CONST_BITS - PASS1_BITS);
dataptr[5] = DESCALE(tmp5 + z2 + z4, CONST_BITS - PASS1_BITS);
dataptr[3] = DESCALE(tmp6 + z2 + z3, CONST_BITS - PASS1_BITS);
dataptr[1] = DESCALE(tmp7 + z1 + z4, CONST_BITS - PASS1_BITS);
dataptr += 8; /* advance pointer to next row */
blkptr += 8;
}
/* Pass 2: process columns.
* We remove the PASS1_BITS scaling, but leave the results scaled up
* by an overall factor of 8.
*/
dataptr = data;
for (i = 0; i < 8; i++) {
tmp0 = dataptr[0] + dataptr[56];
tmp7 = dataptr[0] - dataptr[56];
tmp1 = dataptr[8] + dataptr[48];
tmp6 = dataptr[8] - dataptr[48];
tmp2 = dataptr[16] + dataptr[40];
tmp5 = dataptr[16] - dataptr[40];
tmp3 = dataptr[24] + dataptr[32];
tmp4 = dataptr[24] - dataptr[32];
tmp10 = tmp0 + tmp3;
tmp13 = tmp0 - tmp3;
tmp11 = tmp1 + tmp2;
tmp12 = tmp1 - tmp2;
dataptr[0] = DESCALE(tmp10 + tmp11, PASS1
_BITS);
dataptr[32] = DESCALE(tmp10 - tmp11, PASS1_BITS);
z1 = (tmp12 + tmp13) * FIX_0_541196100;
dataptr[16] = DESCALE(z1 + tmp13 * FIX_0_765366865, CONST_BITS + PASS1_BITS);
dataptr[48] = DESCALE(z1 + tmp12 * (-FIX_1_847759065), CONST_BITS + PASS1_BITS);
z1 = tmp4 + tmp7;
z2 = tmp5 + tmp6;
z3 = tmp4 + tmp6;
z4 = tmp5 + tmp7;
z5 = (z3 + z4) * FIX_1_175875602; /* sqrt(2) * c3 */
tmp4 *= FIX_0_298631336; /* sqrt(2) * (-c1+c3+c5-c7) */
tmp5 *= FIX_2_053119869; /* sqrt(2) * ( c1+c3-c5+c7) */
tmp6 *= FIX_3_072711026; /* sqrt(2) * ( c1+c3+c5-c7) */
tmp7 *= FIX_1_501321110; /* sqrt(2) * ( c1+c3-c5-c7) */
z1 *= -FIX_0_899976223; /* sqrt(2) * (c7-c3) */
z2 *= -FIX_2_562915447; /* sqrt(2) * (-c1-c3) */
z3 *= -FIX_1_961570560; /* sqrt(2) * (-c3-c5) */
z4 *= -FIX_0_390180644; /* sqrt(2) * (c5-c3) */
z3 += z5;
z4 += z5;
dataptr[56] = DESCALE(tmp4 + z1 + z3, CONST_BITS + PASS1_BITS);
dataptr[40] = DESCALE(tmp5 + z2 + z4, CONST_BITS + PASS1_BITS);
dataptr[24] = DESCALE(tmp6 + z2 + z3, CONST_BITS + PASS1_BITS);
dataptr[8] = DESCALE(tmp7 + z1 + z4, CONST_BITS + PASS1_BITS);
dataptr++; /* advance pointer to next column */
}
/* descale */
for (i = 0; i < 64; i++)
block[i] = (double)DESCALE(data[i], 3);
}
// 测试程序：calculate.cpp
#include "FDCT.h"
#include
using namespace std;
int main(void)
{
const short NUM = 8;
short i = 0;
short j = 0;
short u = 0;
short v = 0;
double input[NUM][NUM] =
{
{89.00, 101.00, 114.00, 125.00, 126.00, 115.00, 105.00, 96.00},
{97.00, 115.00, 131.00, 147.00, 149.00, 135.00, 123.00, 113.00},
{114.00, 134.00, 159.00, 178.00, 175.00, 164.00, 149.00, 137.00},
{121.00, 143.00, 177.00, 196.00, 201.00, 189.00, 165.00, 150.00},
{119.00, 141.00, 175.00, 201.00, 207.00, 186.00, 162.00, 144.00},
{107.00, 130.00, 165.00, 189.00, 192.00, 171.00, 144.00, 125.00},
{97.00, 119.00, 149.00, 171.00, 172.00, 145.00, 117.00, 96.00},
{88.00, 107.00, 136.00, 156.00, 155.00, 129.00, 97.00, 75.00}
};
double *data = new double[NUM * NUM];
for(i = 0; i < NUM; i++)
{
for(j = 0; j < NUM; j++)
{
data[i * NUM + j] = input[i][j];
}
}
fdct(data);
// print the result of FDCT
for(u = 0; u < NUM; u++)
{
for(v = 0; v < NUM; v++)
{
printf("%7.2f ", data[u*NUM + v]);
//cout << input[u][v] << '\t';
}
cout << endl;
}
return 0;
}
运算结果：
1125.00 -32.00 -185.00 -7.00 2.00 -1.00 -2.00 2.00
-22.00 -16.00 45.00 -3.00 -2.00 0.00 -1.00 -2.00
-165.00 32.00 17.00 2.00 1.00 -1.00 -3.00 0.00
-7.00 -4.00 0.00 2.00 2.00 -1.00 -1.00 2.00
-2.00 0.00 0.00 3.00 0.00 0.00 2.00 1.00
3.00 1.00 1.00 -1.00 -2.00 1.00 2.00 0.00
0.00 0.00 2.00 -1.00 -1.00 2.0
0 1.00 -1.00
0.00 3.00 1.00 -1.00 2.00 1.00 -2.00 0.00
上述运算结果和离散余弦正逆变换一文中的结果保持完全一致。
尽管两种快速算法的结果有很大的不同，但都反应了信号的频谱特性，因此在实际应用中都可以认为是对的。
本文来自CSDN博客，转载请标明出处：https://www.wendangku.net/doc/783461438.html,/pathuang68/archive/2009/05/19/4199987.aspx
逆：
(https://www.wendangku.net/doc/783461438.html,/pathuang68/archive/2009/05/19/4202267.aspx)
快速逆离散余弦变换代码实现(FIDCT) 收藏
下面的二维快速逆离散余弦变换是在Thomas G. Lane 的代码的基础上，修改而成。可直接用于图像或视频处理。测试代码是玄机逸士写的。// fidct.hvoid fidct(short *const block);void fidct_init();// fidct.cpp/***************************************************************************** * * XVID MPEG-4 VIDEO CODEC * - Inverse DCT - * * These routines are from Independent JPEG Group's free JPEG software * Copyright (C) 1991-1998, Thomas G. Lane (see the file README.IJG) ****************************************************************************//* Copyright (C) 1996, MPEG Software Simulation Group. All Rights Reserved. */#include "fidct.h"#define W1 2841 /* 2048*sqrt(2)*cos(1*pi/16) */#define W2 2676 /* 2048*sqrt(2)*cos(2*pi/16) */#define W3 2408 /* 2048*sqrt(2)*cos(3*pi/16) */#define W5 1609 /* 2048*sqrt(2)*cos(5*pi/16) */#define W6 1108 /* 2048*sqrt(2)*cos(6*pi/16) */#define W7 565 /* 2048*sqrt(2)*cos(7*pi/16) *//* private data */static short iclip[1024]; /* clipping table */static short *iclp;/* two dimensional inverse discrete cosine transform */// fidct_init() MUST BE CALLED BEOFRE THE FIRST CALL TO THIS FUNCTION!void fidct(short *const block){ static short *blk; static long i; static long X0, X1, X2, X3, X4, X5, X6, X7, X8; for (i = 0; i < 8; i++) /* idct rows */ { blk = block + (i << 3); if (!((X1 = blk[4] << 11) | (X2 = blk[6]) | (X3 = blk[2]) | (X4 = blk[1]) | (X5 = blk[7]) | (X6 = blk[5]) | (X7 = blk[3]))) { blk[0] = blk[1] = blk[2] = blk[3] = blk[4] = blk[5] = blk[6] = blk[7] = blk[0] << 3; continue; } X0 = (blk[0] << 11) + 128; /* for proper rounding in the fourth stage */ /* first stage */ X8 = W7 * (X4 + X5); X4 = X8 + (W1 - W7) * X4; X5 = X8 - (W1 + W7) * X5; X8 = W3 * (X6 + X7); X6 = X8 - (W3 - W5) * X6; X7 = X8 - (W3 + W5) * X7; /* second stage */ X8 = X0 + X1; X0 -= X1; X1 = W6 * (X3 + X2); X2 = X1 - (W2 + W6) * X2; X3 = X1 + (W2 - W6) * X3; X1 = X4 + X6; X4 -= X6; X6 = X5 + X7; X5 -= X7; /* third stage */ X7 = X8 + X3; X8 -= X3; X3 = X0 + X2; X0 -= X2; X2 = (181 * (X4 + X5) + 128) >> 8; X4 = (181 * (X4 - X5) + 128) >> 8; /* fourth stage */ blk[0] = (short) ((X7 + X1) >> 8); blk[1] = (short) ((X3 + X2) >> 8); blk[2] = (short) ((X0 + X4) >> 8); blk[3] = (short) ((X8 + X6) >> 8); blk[4] = (short) ((X8 - X6) >> 8); blk[5] = (s
hort) ((X0 - X4) >> 8); blk[6] = (short) ((X3 - X2) >> 8); blk[7] = (short) ((X7 - X1) >> 8); } /* end for ( i = 0; i < 8; ++i ) IDCT-rows */ for (i = 0; i < 8; i++) /* idct columns */ { blk = block + i; /* shortcut */ if (! ((X1 = (blk[8 * 4] << 8)) | (X2 = blk[8 * 6]) | (X3 = blk[8 * 2]) | (X4 = blk[8 *1]) | (X5 = blk[8 * 7]) | (X6 = blk[8 * 5]) | (X7 = blk[8 * 3]))) { blk[8 * 0] = blk[8 * 1] = blk[8 * 2] = blk[8 * 3] = blk[8 * 4] = blk[8 * 5] = blk[8 * 6] = blk[8 * 7] = iclp[(blk[8 * 0] + 32) >> 6]; continue; } X0 = (blk[8 * 0] << 8) + 8192; /* first stage */ X8 = W7 * (X4 + X5) + 4; X4 = (X8 + (W1 - W7) * X4) >> 3; X5 = (X8 - (W1 + W7) * X5) >> 3; X8 = W3 * (X6 + X7) + 4; X6 = (X8 - (W3 - W5) * X6) >> 3; X7 = (X8 - (W3 + W5) * X7) >> 3; /* second stage */ X8 = X0 + X1; X0 -= X1; X1 = W6 * (X3 + X2) + 4; X2 = (X1 - (W2 + W6) * X2) >> 3; X3 = (X1 + (W2 - W6) * X3) >> 3; X1 = X4 + X6; X4 -= X6; X6 = X5 + X7; X5 -= X7; /* third stage */ X7 = X8 + X3; X8 -= X3; X3 = X0 + X2; X0 -= X2; X2 = (181 * (X4 + X5) + 128) >> 8; X4 = (181 * (X4 - X5) + 128) >> 8; /* fourth stage */ blk[8 * 0] = iclp[(X7 + X1) >> 14]; blk[8 * 1] = iclp[(X3 + X2) >> 14]; blk[8 * 2] = iclp[(X0 + X4) >> 14]; blk[8 * 3] = iclp[(X8 + X6) >> 14]; blk[8 * 4] = iclp[(X8 - X6) >> 14]; blk[8 * 5] = iclp[(X0 - X4) >> 14]; blk[8 * 6] = iclp[(X3 - X2) >> 14]; blk[8 * 7] = iclp[(X7 - X1) >> 14]; }} void fidct_init(){ int i; iclp = iclip + 512; for (i = -512; i < 512; i++) iclp[i] = (i < -256) ? -256 : ((i > 255) ? 255 : i);}// 测试代码：calculate.cpp#include
T)
本文来自CSDN博客，转载请标明出处：https://www.wendangku.net/doc/783461438.html,/pathuang68/archive/2009/05/19/4202267.aspx