离散余弦反变换的实现

Question

我现在实现了离散余弦变换.程序如下:
BOOL WINAPI DIBDct(CDib *pDib)
{

// 指向源图像的指针
unsigned char* lpSrc;

//图象的宽度和高度
LONG    lWidth;
LONG    lHeight; // 循环变量
LONG i;
LONG j;

// 离散余弦变换的宽度和高度，必须为2的整数次方
LONG lW = 1;
LONG lH = 1;

int wp = 0;
int hp = 0;

// 中间变量
double dTemp; //得到图象的宽度和高度
CSize   SizeDim;
SizeDim = pDib->GetDimensions();
lWidth = SizeDim.cx;
lHeight = SizeDim.cy; // 图像每行的字节数
LONG lLineBytes;

//得到实际的Dib图象存储大小
CSize   SizeRealDim;
SizeRealDim = pDib->GetDibSaveDim(); // 计算图像每行的字节数
lLineBytes = SizeRealDim.cx;   //图像数据的指针
LPBYTE  lpDIBBits =pDib->m_lpImage;

// 保证离散余弦变换的宽度和高度为2的整数次方
while(lW * 2 <= lWidth)
{
lW = lW * 2;
wp++;
}

while(lH * 2 <= lHeight)
{
lH = lH * 2;
hp++;
} // 分配内存
double *dpf = new double[lW * lH];
double *dpF = new double[lW * lH]; // 时域赋值
for(i = 0; i < lH; i++)
{
for(j = 0; j < lW; j++)
{
// 指针指向位图i行j列的象素
lpSrc = lpDIBBits + lLineBytes * ( i) + j;
// lpSrc = lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
// 将象素值赋给时域数组
dpf[j + i * lW] = *(lpSrc);
}
} for(i = 0; i < lH; i++)
// y方向进行离散余弦变换
DCT(&dpf[lW * i], &dpF[lW * i], wp); // 保存计算结果
for(i = 0; i < lH; i++)
for(j = 0; j < lW; j++)
dpf[j * lH + i] = dpF[j + lW * i];


for(j = 0; j < lW; j++)
// x方向进行离散余弦变换
DCT(&dpf[j * lH], &dpF[j * lH], hp);

// 频谱的计算
for(i = 0; i < lH; i++)
{
for(j = 0; j < lW; j++)
{
dTemp = fabs(dpF[j*lH+i]);

// 是否超过255
if (dTemp > 255)
// 如果超过，设置为255
dTemp = 255;

// 指针指向位图i行j列的象素
lpSrc = lpDIBBits + lLineBytes * ( i) + j;
// lpSrc = lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;

// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
}

// 释放内存
delete dpf;
delete dpF; // 返回
return TRUE;
}
并且我有一维快速离散余弦反变换的函数IDCT,要如何编写离散余弦反变换?
我是新手,请指教!

Answer 1

VOID WINAPI FFT_1D(complex<double> * pCTData, complex<double> * pCFData, int nLevel)
{
// 循环控制变量
int i;
int     j;
int     k; double PI = 3.1415926;  // 傅立叶变换点数
int nCount =0 ; // 计算傅立叶变换点数
nCount =(int)pow(2,nLevel) ;

// 某一级的长度
int nBtFlyLen;
nBtFlyLen = 0 ;

// 变换系数的角度 ＝2 * PI * i / nCount
double dAngle;

complex<double> *pCW ;

// 分配内存，存储傅立叶变化需要的系数表
pCW  = new complex<double>[nCount / 2];    // 计算傅立叶变换的系数
for(i = 0; i < nCount / 2; i++)
{
dAngle = -2 * PI * i / nCount;
pCW[i] = complex<double> ( cos(dAngle), sin(dAngle) );
} // 变换需要的工作空间
complex<double> *pCWork1,*pCWork2; 

// 分配工作空间
pCWork1 = new complex<double>[nCount]; pCWork2 = new complex<double>[nCount];
// 临时变量
complex<double> *pCTmp;

// 初始化，写入数据
memcpy(pCWork1, pCTData, sizeof(complex<double>) * nCount); // 临时变量
int nInter; 
nInter = 0; // 蝶形算法进行快速傅立叶变换
for(k = 0; k < nLevel; k++)
{
for(j = 0; j < (int)pow(2,k); j++)
{
//计算长度
nBtFlyLen = (int)pow( 2,(nLevel-k) );

//倒序重排，加权计算
for(i = 0; i < nBtFlyLen/2; i++)
{
nInter = j * nBtFlyLen;
pCWork2[i + nInter] = 
pCWork1[i + nInter] + pCWork1[i + nInter + nBtFlyLen / 2];
pCWork2[i + nInter + nBtFlyLen / 2] =
(pCWork1[i + nInter] - pCWork1[i + nInter + nBtFlyLen / 2]) 
* pCW[(int)(i * pow(2,k))];
}
} // 交换 pCWork1和pCWork2的数据
pCTmp   = pCWork1 ;
pCWork1 = pCWork2 ;
pCWork2 = pCTmp ;
}

// 重新排序
for(j = 0; j < nCount; j++)
{
nInter = 0;
for(i = 0; i < nLevel; i++)
{
if ( j&(1<<i) )
{
nInter += 1<<(nLevel-i-1);
}
}
pCFData[j]=pCWork1[nInter];
}

// 释放内存空间
delete pCW;
delete pCWork1;
delete pCWork2;
pCW = NULL ;
pCWork1 = NULL ;
pCWork2 = NULL ;}