FFT、IFFT和DCT、IDCT和WALSH、IWALSH

#include "stdafx.h"
#include "math.h"
#include "complex"
#include "iostream"
using namespace std;

 

#define PI 3.1415926

 

void FFT(complex<double>* TD, complex<double>* FD, int r)
{      
       long count;                            // 傅立叶变换点数
       int i, j, k;                            // 循环变量
       int offset, 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++)
       {
              offset = 1 << (r-k);
              for(j = 0; j < 1 << k; j++)
              {
                     p = j * offset;
                     for(i = 0; i < offset / 2; i++)
                     {
                           
                            X2[i + p] = X1[i + p] + X1[i + p + offset / 2];
                            X2[i + p + offset / 2] = (X1[i + p] - X1[i + p + offset / 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;
}

 

void IFFT(complex<double>* FD, complex<double>* TD, int r)
{      
       long count;                                   // 傅立叶变换点数      
       int i;                                                 // 循环变量      
       complex<double> *X;       
       
       count = 1 << r;                                   // 计算傅立叶变换点数
       X = new complex<double>[count];       // 分配运算所需存储器
       memcpy(X, FD, sizeof(complex<double>) * count);       // 将频域点写入X
      
       // 求共轭
       for(i = 0; i < count; i++)
       {
              X[i] = complex<double> (X[i].real(), -X[i].imag());
       }
      
       FFT(X, TD, r);                                   // 调用快速傅立叶变换
      
       // 求时域点的共轭
       for(i = 0; i < count; i++)
       {
              TD[i] = complex<double> (TD[i].real() / count, -TD[i].imag() / count);
       }
      
       delete []X;
}

 

void DCT(double* TD, double* FD, int r)
{      
       long count;                     // 离散余弦变换点数      
       int i;                            // 循环变量      
       double factor;
       complex<double> *X;
             
       count = 1 << r;                     // 计算离散余弦变换点数      

       X = new complex<double>[count * 2];
       memset(X, 0, sizeof(complex<double>) * count * 2);       // 赋初值为
      
       // 将时域点写入数组X
       for(i=0;i<count;i++)
       {
              X[i] = complex<double> (TD[i], 0);
       }
             
       FFT(X, X, r + 1);                            // 调用快速傅立叶变换             
       factor = 1 / sqrt((double)count);              // 调整系数             
       FD[0] = X[0].real() * factor;       // 求F[0]      
       factor *= sqrt(2.0);
      
       // 求F[u]      
       for(i = 1; i < count; i++)
       {
              FD[i] = (X[i].real() * cos(i * PI / (count * 2)) + X[i].imag() *
                     sin(i * PI / (count * 2))) * factor;
       }
      
       delete []X;
}

void IDCT(double* FD, double* TD, int r)
{
       long count;                     // 离散余弦反变换点数
       int i;                            // 循环变量
       double factor, factor0;
       complex<double> *X;
             
       count = 1 << r;                     // 计算离散余弦变换点数

       X = new complex<double>[count * 2];
       memset(X, 0, sizeof(complex<double>) * count * 2);       // 赋初值为
      
       // 将频域变换后点写入数组X
       for(i = 0;i < count;i++)
       {
              X[i] = complex<double> (FD[i] * cos(i * PI / (count * 2)), FD[i] *
                     sin(i * PI / (count * 2)));
       }
      
       IFFT(X, X, r + 1);              // 调用快速傅立叶反变换
      
       // 调整系数
       factor = sqrt(2.0 / count);
       factor0 = (sqrt(1.0 / count) - factor) * FD[0];
      
       for(i = 0; i < count; i++)
       {
              TD[i] = factor0 + X[i].real() * factor * 2 * count;
       }
      
       delete []X;
}

 

void WALSH(double* TD, double* FD, int r)
{      
       long count;                            // 沃尔什-哈达玛变换点数      
       int i, j, k;                            // 循环变量
       int offset, p;
       double *X1, *X2, *X;
             
       count = 1 << r;                            // 计算快速沃尔什变换点数      
       X1 = new double[count];              // 分配运算所需的数组
       X2 = new double[count];              // 分配运算所需的数组
             
       memcpy(X1, TD, sizeof(double) * count);       // 将时域点写入数组X1
      
       // 蝶形运算
       for(k = 0; k < r; k++)
       {
              offset = 1 << (r - k);
              for(j = 0; j < 1 << k; j++)
              {
                     p = j * offset;
                     for(i = 0; i < offset / 2; i++)
                     {                           
                            X2[i + p] = X1[i + p] + X1[i + p + offset / 2];
                            X2[i + p + offset / 2] = X1[i + p] - X1[i + p + offset / 2];
                     }
              }
             
              // 互换X1和X2 
              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] / count;
       }
      
       delete []X1;
       delete []X2;
}

 

void IWALSH(double* FD, double* TD, int r)
{      
       long count;                            // 变换点数      
       int i;                                   // 循环变量
             
       count = 1 << r;                            // 计算变换点数      
       WALSH(FD, TD, r);                            // 调用快速沃尔什－哈达玛变换进行反变换
             
       for(i = 0; i < count; i++)       // 调整系数
       {
              TD[i] *= count;
       }
}