数字图像典型算法- 离散余弦变换

 离散余弦变换
     算法描述：
    离散余弦变换（DCT for Discrete Cosine Transform）是与傅里叶变换相关的一种变换，它类似于离散傅里叶变换（DFT for Discrete Fourier Transform），但是只使用实数。
    离散余弦变换相当于一个长度大概是它两倍的离散傅里叶变换，这个离散傅里叶变换是对一个实偶函数进行的（因为一个实偶函数的傅里叶变换仍然是一个实偶函数），在有些变形里面需要将输入或者输出的位置移动半个单位（DCT有8种标准类型，其中4种是常见的）。

    程序实现：

函数名称：FFT()
参数:
complex<double> * TD - 指向时域数组的指针
complex<double> * FD - 指向频域数组的指针
r         －2的幂数，即迭代次数
返回值:无。
说明:该函数用来实现快速付立叶变换。
VOID WINAPI FFT(complex<double> * TD, complex<double> * FD, int r)
{
 // 付立叶变换点数
 LONG count;
 // 循环变量
 int  i,j,k;
  // 中间变量
 int  bfsize,p;
  // 角度
 double angle;
 
 complex<double> *W,*X1,*X2,*X;
  // 计算付立叶变换点数
 count = 1 << r;
 
 // 分配运算所需存储器
 W  = new complex<double>[count / 2];
 X1 = new complex<double>[count];
 X2 = new complex<double>[count];
 
 // 计算加权系数
 for(i = 0; i < count / 2; i++)
 {
  angle = -i * PI * 2 / count;
  W[i] = complex<double> (cos(angle), sin(angle));
 }
 
 // 将时域点写入X1
 memcpy(X1, TD, sizeof(complex<double>) * count);
 
 // 采用蝶形算法进行快速付立叶变换
 for(k = 0; k < r; k++)
 {
  for(j = 0; j < 1 << k; j++)
  {
   bfsize = 1 << (r-k);
   for(i = 0; i < bfsize / 2; i++)
   {
    p = j * bfsize;
    X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
    X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];
   }
  }
  X  = X1;
  X1 = X2;
  X2 = X;
 }
 
 // 重新排序
 for(j = 0; j < count; j++)
 {
  p = 0;
  for(i = 0; i < r; i++)
  {
   if (j&(1<<i))
   {
    p+=1<<(r-i-1);
   }
  }
  FD[j]=X1[p];
 }
 
 // 释放内存
 delete W;
 delete X1;
 delete X2;
}

函数名称：DCT()
参数:
double * f    - 指向时域值的指针
double * F    - 指向频域值的指针
r      －2的幂数
返回值:无。
说明:该函数用来实现快速离散余弦变换，利用2N点的快速付立叶变换来实现离散余弦变换。
VOID WINAPI DCT(double *f, double *F, int r)
{
 // 离散余弦变换点数
 LONG count;
 // 循环变量
 int  i;
  // 中间变量
 double dTemp;
 
 complex<double> *X;
  // 计算离散余弦变换点数
 count = 1<<r;
 
 // 分配内存
 X = new complex<double>[count*2];
  // 赋初值为0
 memset(X, 0, sizeof(complex<double>) * count * 2);
  // 将时域点写入数组X
 for(i=0;i<count;i++)
 {
  X[i] = complex<double> (f[i], 0);
 }
 
 // 调用快速付立叶变换
 FFT(X,X,r+1);
 // 调整系数
 dTemp = 1/sqrt(count);
  // 求F[0]
 F[0] = X[0].real() * dTemp;
 dTemp *= sqrt(2);
 // 求F[u] 
 for(i = 1; i < count; i++)
 {
  F[i]=(X[i].real() * cos(i*PI/(count*2)) + X[i].imag() * sin(i*PI/(count*2))) * dTemp;
 }
 
 // 释放内存
 delete X;
}

函数名称：DIBDct()
参数:
LPSTR lpDIBBits    - 指向源DIB图像指针
LONG  lWidth       - 源图像宽度（象素数）
LONG  lHeight      - 源图像高度（象素数）
返回值:BOOL               - 成功返回TRUE，否则返回FALSE。
说明:该函数用来对图像进行离散余弦变换。
BOOL WINAPI DIBDct(LPSTR lpDIBBits, LONG lWidth, LONG lHeight)
{
  // 指向源图像的指针
 unsigned char* lpSrc;
  // 循环变量
 LONG i;
 LONG j;
  // 进行付立叶变换的宽度和高度（2的整数次方）
 LONG w;
 LONG h;
  // 中间变量
 double dTemp;
 int  wp;
 int  hp;
 
 // 图像每行的字节数
 LONG lLineBytes;
  // 计算图像每行的字节数
 lLineBytes = WIDTHBYTES(lWidth * 8);
 
 // 赋初值
 w = 1;
 h = 1;
 wp = 0;
 hp = 0;
 
 // 计算进行离散余弦变换的宽度和高度（2的整数次方）
 while(w * 2 <= lWidth)
 {
  w *= 2;
  wp++;
 }
 
 while(h * 2 <= lHeight)
 {
  h *= 2;
  hp++;
 }
  // 分配内存
 double *f = new double[w * h];
 double *F = new double[w * h];
 
 // 行
 for(i = 0; i < h; i++)
 {
  // 列
  for(j = 0; j < w; j++)
  {
   // 指向DIB第i行，第j个象素的指针
   lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
   
   // 给时域赋值
   f[j + i * w] = *(lpSrc);
  }
 }
 
 for(i = 0; i < h; i++)
 {
  // 对y方向进行离散余弦变换
  DCT(&f[w * i], &F[w * i], wp);
 }
 
 // 保存计算结果
 for(i = 0; i < h; i++)
 {
  for(j = 0; j < w; j++)
  {
   f[j * h + i] = F[j + w * i];
  }
 }
 
 for(j = 0; j < w; j++)
 {
  // 对x方向进行离散余弦变换
  DCT(&f[j * h], &F[j * h], hp);
 }
  // 行
 for(i = 0; i < h; i++)
 {
  // 列
  for(j = 0; j < w; j++)
  {
   // 计算频谱
   dTemp = fabs(F[j*h+i]);
   
   // 判断是否超过255
   if (dTemp > 255)
   {
    // 对于超过的，直接设置为255
    dTemp = 255;
   }
   
   // 指向DIB第y行，第x个象素的指针
   lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
   
   // 更新源图像
   * (lpSrc) = (BYTE)(dTemp);
  }
 }
 
 // 释放内存
 delete f;
 delete F;

 // 返回
 return TRUE;
}