离散余弦正逆变换

最近网上有朋友因为要进行图像处理，而问及离散余弦正逆变换的问题，我在网上找了一些资料，网下也找了一些参考书，发现居然有很多错误，尤其是逆变换(IDCT)的公式都是错的！怪不得网友做完DCT后，再IDCT，却得不到原来的出使数据了(四舍五入的误差除外)，实在误人甚深。写此博文，以正其谬。

说明：

下面实现的代码，严格遵循离散余弦正逆变换的公式，仅供理解离散余弦逆变换的原理而写。并未考虑代码优化，也未考虑快速算法，快速算法见快速离散余弦变换代码实现(FDCT)和快速逆离散余弦变换代码实现(FIDCT)。

离散余弦变换(DCT)公式：

离散余弦逆变换(IDCT)公式：

注：严格的离散余弦变换公式中，其中一个N可以写为M，且M和N可以不相等。这里由于为图像处理服务，所以都写成了N。

代码实现：

#include <iostream>

#include <math.h>

using namespace std;

#define NUM 8

#define PI 3.1415926

short round(double a)

{

if (a >= 0)

{

return (short)(a + 0.5);

}

else

{

return (short)(a - 0.5);

}

}

// DCT - Discrete Cosine Transform

void DCT(short data[NUM][NUM])

{

short output[NUM][NUM];

double ALPHA, BETA;

short u = 0;

short v = 0;

short i = 0;

short j = 0;

for(u = 0; u < NUM; u++)

{

for(v = 0; v < NUM; v++)

{

if(u == 0)

{

ALPHA = sqrt(1.0 / NUM);

}

else

{

ALPHA = sqrt(2.0 / NUM);

}

if(v == 0)

{

BETA = sqrt(1.0 / NUM);

}

else

{

BETA = sqrt(2.0 / NUM);

}

double tmp = 0.0;

for(i = 0; i < NUM; i++)

{

for(j = 0; j < NUM; j++)

{

tmp += data[i][j] * cos((2*i+1)*u*PI/(2.0 * NUM)) * cos((2*j+1)*v*PI/(2.0 * NUM));

}

}

output[u][v] = round(ALPHA * BETA * tmp);

}

}

memset(data, 0, NUM * NUM * sizeof(short));

memcpy(data, output, NUM * NUM * sizeof(short));

}

// Inverse DCT

void IDCT(short data[NUM][NUM])

{

short output[NUM][NUM];

double ALPHA, BETA;

short u = 0;

short v = 0;

short i = 0;

short j = 0;

for(i = 0; i < NUM; i++)

{

for(short j = 0; j < NUM; j++)

{

double tmp = 0.0;

for(short u = 0; u < NUM; u++)

{

for(v = 0; v < NUM; v++)

{

if(u == 0)

{

ALPHA = sqrt(1.0 / NUM);

}

else

{

ALPHA = sqrt(2.0 / NUM);

}

if(v == 0)

{

BETA = sqrt(1.0 / NUM);

}

else

{

BETA = sqrt(2.0 / NUM);

}

tmp += ALPHA * BETA * data[u][v] * cos((2*i+1)*u*PI/(2.0 * NUM)) * cos((2*j+1)*v*PI/(2.0 * NUM));

}

}

output[i][j] = round(tmp);

}

}

memset(data, 0, NUM * NUM * sizeof(short));

memcpy(data, output, NUM * NUM * sizeof(short));

}

int main(void)

{

short i = 0;

short j = 0;

short u = 0;

short v = 0;

// 8 x 8 的图像数据

short input[NUM][NUM] =

{

{89, 101, 114, 125, 126, 115, 105, 96},

{97, 115, 131, 147, 149, 135, 123, 113},

{114, 134, 159, 178, 175, 164, 149, 137},

{121, 143, 177, 196, 201, 189, 165, 150},

{119, 141, 175, 201, 207, 186, 162, 144},

{107, 130, 165, 189, 192, 171, 144, 125},

{97, 119, 149, 171, 172, 145, 117, 96},

{88, 107, 136, 156, 155, 129, 97, 75}

};

DCT(input);

cout << "The result of DCT:" << endl;

for(u = 0; u < NUM; u++)

{

for(v = 0; v < NUM; v++)

{

cout << input[u][v] << '/t';

}

cout << endl;

}

IDCT(input);

cout << "The result of IDCT:" << endl;

for(i = 0; i < NUM; i++)

{

for(short j = 0; j < NUM; j++)

{

cout << input[i][j] << '/t';

}

cout << endl;

}

return 0;

}

运算结果如下：

The result of DCT:
1125 -32 -185 -72 -1 -2 2
-22 -16 45 -3 -2 0 -2 -2
-165 3217 2 1-1-3 0
-7 -4 0 2 2-1 -1 2
-2 0 0 3 0 0 2 1
3 1 1 -1-2 0 2 0
0 0 2 -1-1 2 1 -1
0 3 1 -1 2 1 -2 0

The result of IDCT:
89 101 114 125 126 115 105 96
97 115 131 147 149 135 123 113
114 134 159 178 175 164 149 137
121 143 177 196 201 189 165 150
119 141 175 201 207 186 162 144
107 130 165 189 192 171 144 125
96 119 150 171 172 145 116 96
88 107 136 156 155 129 97 75