离散傅里叶的MFC算法详解（DFT，FFT，DCT）

算法详解：离散时间傅里叶变换DFT 

转载▼源程序来自某书

---------------------
BOOL CDibImage::Fourier(LPSTR lpDIBBits, LONG lWidth, LONG lHeight)
{        
    unsigned char*    lpSrc;    // 指向源图像的指针
    double    dTemp;    
    LONG    i;                // 循环变量
    LONG    j;    
    LONG    w;                // 进行傅里叶变换的宽度（2的整数次方）
    LONG    h;                // 进行傅里叶变换的高度（2的整数次方）    
    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++;//幂
    }   
    complex<double> *TD = new complex<double>[w * h];
    complex<double> *FD = new complex<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;           
            // 给时域赋值
            TD[j + w * i] = complex<double>(*(lpSrc), 0);
        }
    }   
    for(i = 0; i < h; i++)
    {
        // 对y方向进行快速傅里叶变换
        FFT(&TD[w * i], &FD[w * i], wp);
    }   
    // 保存变换结果
    for(i = 0; i < h; i++)
    {
        for(j = 0; j < w; j++)
        {
            TD[i + h * j] = FD[j + w * i];//y变x，完成y方向一维运算，导入横方向做第二维运算。
        }
    }   
    for(i = 0; i < w; i++)
    {
        // 对x方向进行快速傅里叶变换
        FFT(&TD[i * h], &FD[i * h], hp);
    }       
    for(i = 0; i < h; i++)            // 行
    {        
        for(j = 0; j < w; j++)        // 列
        {
            // 计算频谱，求模
            dTemp = sqrt(FD[j * h + i].real() * FD[j * h + i].real() + FD[j * h + i].imag() * FD[j * h + i].imag()) / 100;
            if (dTemp > 255)
            {
                dTemp = 255;
            }           
            // 指向DIB第(i<h/2 ? i+h/2 : i-h/2)行，第(j<w/2 ? j+w/2 : j-w/2)个
            // 象素的指针，此处不直接取i和j，是为了将变换后的原点移到中心
            // lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight-1-i) + j;
            lpSrc = (unsigned char*)lpDIBBits + lLineBytes *
                (lHeight - 1 - (i<h/2 ? i+h/2 : i-h/2)) + (j<w/2 ? j+w/2 : j-w/2);//平移(h/2,w/2)
            // 更新源图像
            * (lpSrc) = (BYTE)(dTemp);
        }
    }   
    delete TD;
    delete FD;
    return TRUE;
}

算法详解：基2时间抽取FFT    

转载

这是我修改出来的第一个算法，花费了不少时间。主要是卡在c语言和变量传递上。
---------------------------------------------------------
VOID CDibImage::FFTTIME(complex<double> * TD, complex<double> * FD, int r)
{    
    LONG    count;                // 付立叶变换总点数    
    int        i,j,k;                // 循环变量
    int        N,p;                //N为当前每组点数
    double    w;                    // 弧度    
    complex<double> *W,*X1,*X2,*X,*z;
    
    count = 1 << r;                // 计算付立叶变换总点数，r为总幂
    
    // 分配运算所需存储器
    W  = new complex<double>[count/2];                    //count为接收点数(256，512等)
    X1 = new complex<double>[count];
    X2 = new complex<double>[count];
    
    // 计算加权系数(旋转因子W),w为弧度
    for(i = 0; i < count/2; i++)
    {
        w = -i * PI * 2 / count;
        W[i] = complex<double> (cos(w), sin(w));            //cos实部
    }   
    // 将时域点写入X1,X1为输入，X2为输出
//    memcpy(X1, TD, sizeof(complex<double>) * count);
    
// 重新排序,从最低位到最高位检查正序数的每一位i，然后把该位的值赋给倒序数的相应位
    //在倒序数中为第r-i-1位，r为序数的位数。二进制正序右边为低位
        for(j = 0; j < count; j++)
    {
        p = 0;
        for(i = 0; i < r; i++)    //二进制的i位，0...r-1
        {    
            if (j&(1<<i))         //j按位与位数，若非0则进入循环
            {   
                p+=1<<(r-i-1);    //如果正序数中位i的值为1，则倒序数中r-i-1的相应位置1。p+是累加其他所有位的1
            }
        }
           X1[j]=TD[p];        //对位       
    }
    // 快速付立叶变换(基2时间抽取)，k对调r-k-1
    for(k =0 ; k < r; k++)                                  //k为当前幂，k=0为第1级
    {          
            N = 1 <<k+1;    //N为当前幂每组点数，频率抽取第一级最多,k=0,N=2^r，512-256-128....2
            for(i = 0; i < N / 2; i++)          //i为当前幂每组半点，所以后面两个半组一起算
            {
              for(j = 0; j < count/N; j++)//j为当前幂的组数，k=0时第1级有1组(j=0)，k=1时第2级有2组
                {
                p = j * N;                    //p为当前幂组外间隔，即算到第j组时跨过的总点数
                X2[i + p] = X1[i + p] + X1[i + p + N / 2]* W[i * (1<<r-k-1)];  //每组上半，N/2为间组内隔，256-128....1                
                X2[i + p + N / 2] = X1[i + p] - X1[i + p + N / 2] * W[i * (1<<r-k-1)];//每组下半，[]内为W上标,为每组半点*组数
                }
            }       
        X  = X1;
        X1 = X2;
        X2 = X;
    }       
    for(j = 0; j < count; j++)
    {   
        FD[j]=X1[j];              //对位
    }   
    delete W;
    delete X1;
    delete X2;
}