📄 dct.c
字号:
/************************************************************************ *DCT快速变换 ************************************************************************/#include "macros.h"#include "sim.h"#include <math.h>#ifndef PI# ifdef M_PI# define PI M_PI# else# define PI 3.14159265358979323846# endif#endifint zigzag[8][8] = { {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},};/********************************************************************** * * 函数名: Dct * Description: 8×8编码块的快速DCT变换,并利用zigzag扫描完成系数赋值 * * 输入: 包含64个像素的一维数组 * 返回值: 64个变换系数,一维数组形式 **********************************************************************/int Dct( int *block, int *coeff){ int j1, i, j, k; float b[8]; float b1[8]; float d[8][8]; float f0=(float).7071068; float f1=(float).4903926; float f2=(float).4619398; float f3=(float).4157348; float f4=(float).3535534; float f5=(float).2777851; float f6=(float).1913417; float f7=(float).0975452; //变换源数据赋值 for (i = 0, k = 0; i < 8; i++, k += 8) { for (j = 0; j < 8; j++) { b[j] = (float)block[k+j]; } /* 水平方向(行)变换 */ for (j = 0; j < 4; j++) { j1 = 7 - j; b1[j] = b[j] + b[j1]; b1[j1] = b[j] - b[j1]; }
//按照递减的顺利进行计算 b[0] = b1[0] + b1[3]; b[1] = b1[1] + b1[2]; b[2] = b1[1] - b1[2]; b[3] = b1[0] - b1[3]; b[4] = b1[4]; b[5] = (b1[6] - b1[5]) * f0; b[6] = (b1[6] + b1[5]) * f0; b[7] = b1[7];
//偶数项的提升处理 d[i][0] = (b[0] + b[1]) * f4; d[i][4] = (b[0] - b[1]) * f4; d[i][2] = b[2] * f6 + b[3] * f2; d[i][6] = b[3] * f6 - b[2] * f2; b1[4] = b[4] + b[5]; b1[7] = b[7] + b[6]; b1[5] = b[4] - b[5]; b1[6] = b[7] - b[6];
//奇数项的提升处理 d[i][1] = b1[4] * f7 + b1[7] * f1; d[i][5] = b1[5] * f3 + b1[6] * f5; d[i][7] = b1[7] * f7 - b1[4] * f1; d[i][3] = b1[6] * f3 - b1[5] * f5; } /* 竖直方向(列)变换 */ for (i = 0; i < 8; i++) { for (j = 0; j < 4; j++) { j1 = 7 - j; b1[j] = d[j][i] + d[j1][i]; b1[j1] = d[j][i] - d[j1][i]; }
//与行变换进行相同的处理 b[0] = b1[0] + b1[3]; b[1] = b1[1] + b1[2]; b[2] = b1[1] - b1[2]; b[3] = b1[0] - b1[3]; b[4] = b1[4]; b[5] = (b1[6] - b1[5]) * f0; b[6] = (b1[6] + b1[5]) * f0; b[7] = b1[7];
//偶数项的提升 d[0][i] = (b[0] + b[1]) * f4; d[4][i] = (b[0] - b[1]) * f4; d[2][i] = b[2] * f6 + b[3] * f2; d[6][i] = b[3] * f6 - b[2] * f2; b1[4] = b[4] + b[5]; b1[7] = b[7] + b[6]; b1[5] = b[4] - b[5]; b1[6] = b[7] - b[6];
//奇数项的提升 d[1][i] = b1[4] * f7 + b1[7] * f1; d[5][i] = b1[5] * f3 + b1[6] * f5; d[7][i] = b1[7] * f7 - b1[4] * f1; d[3][i] = b1[6] * f3 - b1[5] * f5; } /* Zigzag扫描并输出变换系数 */ for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) { *(coeff + zigzag[i][j]) = (int)(d[i][j]); } } return 0;}#ifdef FASTIDCT/********************************************************************** * * Name: idct * Description: Descans zigzag-scanned coefficients and does * inverse dct on 64 coefficients * single precision floats * * Input: 64 coefficients, block for 64 pixels * Returns: 0 * Side effects: * * Date: 930128 Author: Robert.Danielsen@nta.no * **********************************************************************/int idct(int *coeff,int *block){ int j1, i, j; double b[8], b1[8], d[8][8]; double f0=.7071068; double f1=.4903926; double f2=.4619398; double f3=.4157348; double f4=.3535534; double f5=.2777851; double f6=.1913417; double f7=.0975452; double e, f, g, h; /* Horizontal */ /* Descan coefficients first */ for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) { b[j] = *( coeff + zigzag[i][j]); } e = b[1] * f7 - b[7] * f1; h = b[7] * f7 + b[1] * f1; f = b[5] * f3 - b[3] * f5; g = b[3] * f3 + b[5] * f5; b1[0] = (b[0] + b[4]) * f4; b1[1] = (b[0] - b[4]) * f4; b1[2] = b[2] * f6 - b[6] * f2; b1[3] = b[6] * f6 + b[2] * f2; b[4] = e + f; b1[5] = e - f; b1[6] = h - g; b[7] = h + g; b[5] = (b1[6] - b1[5]) * f0; b[6] = (b1[6] + b1[5]) * f0; b[0] = b1[0] + b1[3]; b[1] = b1[1] + b1[2]; b[2] = b1[1] - b1[2]; b[3] = b1[0] - b1[3]; for (j = 0; j < 4; j++) { j1 = 7 - j; d[i][j] = b[j] + b[j1]; d[i][j1] = b[j] - b[j1]; } } /* Vertical */ for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) { b[j] = d[j][i]; } e = b[1] * f7 - b[7] * f1; h = b[7] * f7 + b[1] * f1; f = b[5] * f3 - b[3] * f5; g = b[3] * f3 + b[5] * f5; b1[0] = (b[0] + b[4]) * f4; b1[1] = (b[0] - b[4]) * f4; b1[2] = b[2] * f6 - b[6] * f2; b1[3] = b[6] * f6 + b[2] * f2; b[4] = e + f; b1[5] = e - f; b1[6] = h - g; b[7] = h + g; b[5] = (b1[6] - b1[5]) * f0; b[6] = (b1[6] + b1[5]) * f0; b[0] = b1[0] + b1[3]; b[1] = b1[1] + b1[2]; b[2] = b1[1] - b1[2]; b[3] = b1[0] - b1[3]; for (j = 0; j < 4; j++) { j1 = 7 - j; d[j][i] = b[j] + b[j1]; d[j1][i] = b[j] - b[j1]; } } for (i = 0; i < 8; i++) { for (j = 0; j < 8; j++) { *(block + i * 8 + j) = mnint(d[i][j]); } } return 0;}#else/* Perform IEEE 1180 reference (64-bit floating point, separable 8x1 * direct matrix multiply) Inverse Discrete Cosine Transform*//* Here we use math.h to generate constants. Compiler results may vary a little *//* private data *//* cosine transform matrix for 8x1 IDCT */static double c[8][8];/* initialize DCT coefficient matrix */void init_idctref(){ int freq, time; double scale; for (freq=0; freq < 8; freq++) { scale = (freq == 0) ? sqrt(0.125) : 0.5; for (time=0; time<8; time++) c[freq][time] = scale*cos((PI/8.0)*freq*(time + 0.5)); }}/* perform IDCT matrix multiply for 8x8 coefficient block */void idctref(int *coeff, int *block){ int i, j, k, v; double partial_product; double tmp[64]; int tmp2[64]; extern int zigzag[8][8]; for (i=0; i<8; i++) for (j=0; j<8; j++) tmp2[j+i*8] = *(coeff + zigzag[i][j]); for (i=0; i<8; i++) for (j=0; j<8; j++) { partial_product = 0.0; for (k=0; k<8; k++) partial_product+= c[k][j]*tmp2[8*i+k]; tmp[8*i+j] = partial_product; } /* Transpose operation is integrated into address mapping by switching loop order of i and j */ for (j=0; j<8; j++) for (i=0; i<8; i++) { partial_product = 0.0; for (k=0; k<8; k++) partial_product+= c[k][i]*tmp[8*k+j]; v = (int)floor(partial_product+0.5); block[8*i+j] = (v<-256) ? -256 : ((v>255) ? 255 : v); }}#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -