离散余弦变换C语言实现(DCT)
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
void initDCTParam(int deg)
{
// deg 为DCT变换数据长度的幂
if(bHasInit)
{
return; //不用再计算查找表
}
int total, halftotal, i, group, endstart, factor;
total = 1 << deg;
if(C != NULL) delete []C;
C = (double *)new double[total];
halftotal = total >> 1;
for(i=0; i < halftotal; i++)
C[total-i-1]=(double)(2*i+1);
for(group=0; group < deg-1; group++)
{
endstart=1 << (deg-1-group);
int len = endstart >> 1;
factor=1 << (group+1);
for(int j = 0;j < len; j++)
C[endstart-j-1] = factor*C[total-j-1];
}
for(i=1; i < total; i++)
C[i] = 2.0*cos(C[i]*PI/(total << 1)); ///C[0]空着,没使用
bHasInit=true;
}
DCT变换过程可分为两部分:前向追底和后向回根
前向追底:
void dct_forward(double *f,int deg)
{
// f中存储DCT数据
int i_deg, i_halfwing, total, wing, wings, winglen, halfwing;
double temp1,temp2;
total = 1 << deg;
for(i_deg = 0; i_deg < deg; i_deg++)
{
wings = 1 << i_deg;
winglen = total >> i_deg;
halfwing = winglen >> 1;
for(wing = 0; wing < wings; wing++)
{
for(i_halfwing = 0; i_halfwing < halfwing; i_halfwing++)
{
temp1 = f[wing*winglen+i_halfwing];
temp2 = f[(wing+1)*winglen-1-i_halfwing];
if(wing%2)
swap(temp1,temp2); // 交换temp1与temp2的值
f[wing*winglen+i_halfwing] = temp1+temp2;
f[(wing+1)*winglen-1-i_halfwing] =
(temp1-temp2)*C[winglen-1-i_halfwing];
}
}
}
}
后向回根: void dct_backward(double *f,int deg)
{
// f中存储DCT数据
int total,i_deg,wing,wings,halfwing,winglen,i_halfwing,temp1,temp2;
total = 1 << deg;
for(i_deg = deg-1; i_deg >= 0; i_deg--)
{
wings = 1 << i_deg;
winglen = 1 << (deg-i_deg);
halfwing = winglen >> 1;
for(wing = 0; wing < wings; wing++)
{
for(i_halfwing = 0; i_halfwing < halfwing; i_halfwing++)
{
//f[i_halfwing+wing*winglen] = f[i_halfwing+wing*winglen];
if(i_halfwing == 0)
{
f[halfwing+wing*winglen+i_halfwing] =
0.5*f[halfwing+wing*winglen+i_halfwing];
}
else
{
temp1=bitrev(i_halfwing,deg-i_deg-1); // bitrev为位反序
temp2=bitrev(i_halfwing-1,deg-i_deg-1); // 第一参数为要变换的数
// 第二参数为二进制长度
f[halfwing+wing*winglen+temp1] =
f[halfwing+wing*winglen+temp1]-f[halfwing+wing*winglen+temp2];
}
}
}
}
}
位反序函数如下:int bitrev(int bi,int deg)
{
int j = 1, temp = 0, degnum, halfnum;
degnum = deg;
//if(deg<0) return 0;
if(deg == 0) return bi;
halfnum = 1 << (deg-1);
while(halfnum)
{
if(halfnum&bi)
temp += j;
j<<=1;
halfnum >>= 1;
}
return temp;
}
完整实现一维DCT变换:void fdct_1D_no_param(double *f,int deg)
{
initDCTParam(deg);
dct_forward(f,deg);
dct_backward(f,deg);
fbitrev(f,deg); // 实现位反序排列
f[0] = 1/(sqrt(2.0))*f[0];
}
void fdct_1D(double *f,int deg)
{
fdct_1D_no_param(f,deg);
int total = 1 << deg;
double param = sqrt(2.0/total);
for(int i = 0; i < total; i++)
f[i] = param*f[i];
}
利用一维DCT变换来实现二维DCT变换:void fdct_2D(double *f,int deg_row,int deg_col)
{
int rows,cols,i_row,i_col;
double two_div_sqrtcolrow;
rows=1 << deg_row;
cols=1 << deg_col;
init2D_Param(rows,cols);
two_div_sqrtcolrow = 2.0/(sqrt(double(rows*cols)));
for(i_row = 0; i_row < rows; i_row++)
{
fdct_1D_no_param(f+i_row*cols,deg_col);
}
for(i_col = 0; i_col < cols; i_col++)
{
for(i_row = 0; i_row < rows; i_row++)
{
temp_2D[i_row] = f[i_row*cols+i_col];
}
fdct_1D_no_param(temp_2D, deg_row);
for(i_row = 0; i_row < rows; i_row++)
{
f[i_row*cols+i_col] = temp_2D[i_row]*two_div_sqrtcolrow;
}
}
bHasInit = false;
}
IDCT快速变换
初始化查找表: void initIDCTParam(int deg)
{
if(bHasInit)
return; //不用再计算查找表
int total, halftotal, i, group, endstart, factor;
total = 1 << deg;
// if(C!=NULL) delete []C;
// C=(double *)new double[total];
// 由于正变换已经为C申请了空间,所以逆变换就需再申请空间了!
halftotal = total >> 1;
for(i = 0; i < halftotal; i++)
C[total-i-1] = (double)(2*i+1);
for(group = 0; group < deg-1; group++)
{
endstart = 1 << (deg-1-group);
int len = endstart>>1;
factor = 1 << (group+
1);
for(int j = 0; j < len; j++)
C[endstart-j-1] = factor*C[total-j-1];
}
for(i = 1; i < total; i++)
C[i] = 1.0/(2.0*cos(C[i]*PI/(total << 1))); // C[0]空着没用
bHasInit=true;
}
IDCT变换过程也可分为两部分:前向追底和后向回根
前向追底 void idct_forward(double *F,int deg)
{
int total,i_deg,wing,wings,halfwing,winglen,i_halfwing,temp1,temp2;
total = 1 << deg;
for(i_deg = 0; i_deg < deg; i_deg++)
{
wings = 1 << i_deg;
winglen = 1 << (deg-i_deg);
halfwing = winglen >> 1;
for(wing = 0; wing < wings; wing++)
{
for(i_halfwing = halfwing-1; i_halfwing >= 0; i_halfwing--)
{
if(i_halfwing == 0)
{
F[halfwing+wing*winglen+i_halfwing] =
2.0*F[halfwing+wing*winglen+i_halfwing];
}
else
{
temp1 = bitrev(i_halfwing,deg-i_deg-1);
temp2 = bitrev(i_halfwing-1,deg-i_deg-1);
F[halfwing+wing*winglen+temp1] = F[halfwing+wing*winglen+temp1]
+F[halfwing+wing*winglen+temp2];
}
}
}
}
}
后向回根void idct_backward(double *F, int deg)
{
int i_deg,i_halfwing,total,wing,wings,winglen,halfwing;
double temp1, temp2;
total = 1 << deg;
for(i_deg = deg-1; i_deg >= 0; i_deg--)
{
wings = 1 << i_deg;
winglen = total >> i_deg;
halfwing = winglen >> 1;
for(wing = 0; wing < wings; wing++)
{
for(i_halfwing = 0; i_halfwing < halfwing; i_halfwing++)
{
temp1 = F[wing*winglen+i_halfwing];
temp2 = F[(wing+1)*winglen-1-i_halfwing]*C[winglen-1-i_halfwing];
if(wing % 2)
{
F[wing*winglen+i_halfwing] = (temp1-temp2)*0.5;
F[(wing+1)*winglen-1-i_halfwing] = (temp1+temp2)*0.5;
}
else
{
F[wing*winglen+i_halfwing] = (temp1+temp2)*0.5;
F[(wing+1)*winglen-1-i_halfwing] = (temp1-temp2)*0.5;
}
}
}
}
}
完整实现一维IDCT变换:void fidct_1D_no_param(double *F, int deg)
{
initIDCTParam(deg);
F[0] = F[0]*sqrt(2.0);
fbitrev(F, deg);
idct_forward
(F, deg);
idct_backward(F, deg);
}
void fdct_1D(double *f, int deg)
{
fdct_1D_no_param(f, deg);
int total = 1 << deg;
double param = sqrt(2.0/total);
for(int i = 0; i < total; i++)
f[i] = param*f[i];
}
利用一维IDCT变换来实现二维IDCT变换:void fidct_2D(double *F, int deg_row, int deg_col)
{
int rows,cols,i_row,i_col;
double sqrtcolrow_div_two;
rows = 1 << deg_row;
cols = 1 << deg_col;
init2D_Param(rows,cols);
sqrtcolrow_div_two = (sqrt(double(rows*cols)))/2.0;
for(i_row = 0; i_row < rows; i_row++)
{
fidct_1D_no_param(F+i_row*cols,deg_col);
}
for(i_col = 0; i_col < cols; i_col++)
{
for(i_row = 0; i_row < rows; i_row++)
{
temp_2D[i_row] = F[i_row*cols+i_col];
}
fidct_1D_no_param(temp_2D, deg_row);
for(i_row = 0; i_row < rows; i_row++)
{
F[i_row*cols+i_col] = temp_2D[i_row]*sqrtcolrow_div_two;
}
}
bHasInit=false;
}
多线程的考量由于DCT变换要花费一定的时间,特别是在数据矩阵尺寸比较大的时候。此时,如果没有增加一个线程来执行DCT变换,操作界面可能因程序忙于DCT变换的计算而失去响应,所以,增加一个用来进行DCT变换的线程是十分必要的。
首先定义一个结构 typedef struct
{
int row;
int col;
double *data;
//double *data2;
//double *data3; // 在计算彩色图象的数据矩阵时,彩色图象用RGB三个分量
bool m_bfinished;
DWORD m_start,m_end; //以毫秒计,用来计算DCT变换所用的时间;
}RUNINFO;
DCT变换进程函数:UINT ThreadProcfastDct(LPVOID pParam)
{
RUNINFO *pinfo = (RUNINFO*)pParam;
pinfo->m_start = ::GetTickCount();
fdct_2D((double *)pinfo->data, GetTwoIndex(pinfo->row), GetTwoIndex(pinfo->col));
pinfo->m_end = ::GetTickCount();
pinfo->m_bfinished = true;
return 1;
}
IDCT变换进程函数:UINT ThreadProcfastIDct(LPVOID pParam)
{
RUNINFO *pinfo = (RUNINFO*)pParam;
pinfo->m_start = ::GetTickCount();
fidct_2D((double *)pinfo->data, GetTwoIndex(pinfo->row), GetTwoIndex(pinfo->col));
pinfo->m_end = ::GetTickCount();
pinfo->m_bfinished = true;
return 1;
}