离散傅里叶变换的实现与在图像中的应用

c++实现二维离散傅里叶变换：

基于快速傅里叶蝶形算法

#include <stdio.h>
#include <iostream>
#include <complex>
#include <bitset>
#include <fstream>
#include <algorithm>
#include <iterator>
#include <math.h>
#include "cv.h"
#include "highgui.h"
#include "CImg.h"
#define pi 3.1415926535
using std::iostream;
using std::bitset;
using std::complex;
using namespace std;
using namespace cimg_library;

int rev(int k,int n)
{
     bitset<32> tmp(k);
     bitset<32> dist;
     for(int m = 0;m<n;m++)
     {
        if(tmp.test(m))
        {
            dist.set(n-1-m);
        }
     }
     int revk = (int)dist.to_ulong();
     return revk;
}
//重载形式
complex<double>* FFT(const complex<double> *srcimg,int n)
{
    double flag = log10((double)n)/log10(2.0);
    int N = n;
    if(flag - (int)flag != 0){   //判断是否为2的指数次
        cout <<"the length of srcimg is wrong"<< endl;
        /*填充成2的指数项*/
        cout <<"need padding"<<endl;
        N = pow(2,(double)((int)flag+1));
        flag = log10((double)N)/log10(2.0);
    }
    /*改变成奇偶顺序*/
    complex<double> *arr= new complex<double>[N];
    int sub;
    for(int k =0;k<N;k++)
    {
        sub =rev(k,(int)flag); 
        if(sub <= n-1){
            arr[k] = *(srcimg + sub);
        }else{
            complex<double> t = complex<double>(0,0);
            arr[k] = t;
        }
    }
//    cout<<"------------after padding and retrival-----------------"<<endl;
//    for(size_t y=0;y<N;y++)
//    {
//        cout << arr[y].real()<<"  "<<arr[y].imag()<<endl;
//    }

    /*基于迭代的蝶形快速傅立叶变换，自底向上*/
     for(int s =1;s <= flag;s++)
    {
        int m = pow(2,(double)s);
        complex<double> wm = complex<double>(cos(2*pi/m),sin(2*pi/m));//wm1：从W21开始,周期变换
        for(int p = 0;p<N;p+=m)
        {
            complex<double> w(1,0);
            for(int j = 0;j<=m/2-1;j++)
            {
                complex<double> t = w * arr[p+j+m/2];
                complex<double> u = arr[p+j];
                arr[p+j] = u+t;
                arr[p+j+m/2] = u-t;
                w = w*wm;
            }
        }
    }
    return arr;
}
/***********一维快速傅立叶变换********************
*srcimg : 原始一维序列                          *
*n     ：一维序列的长度
*************************************************/
complex<double>* FFT(const double *srcimg,int n)
{
    double flag = log10((double)n)/log10(2.0);
    int N = n;
    if(flag - (int)flag != 0){   //判断是否为2的指数次
//        cout <<"the length of srcimg is wrong"<< endl;
        /*填充成2的指数项*/
//        cout <<"need padding"<<endl;
        N = pow(2,(double)((int)flag+1));
        flag = log10((double)N)/log10(2.0);
    }
    /*改变成奇偶顺序*/
    complex<double> *arr= new complex<double>[N];
    int sub;
    for(int k =0;k<N;k++)
    {
        sub =rev(k,(int)flag); 
        if(sub <= n-1){
            arr[k] = complex<double>(*(srcimg + sub),0);
        }else{
            complex<double> t = complex<double>(0,0);
            arr[k] = t;
        }
    }
//    cout<<"------------after padding and retrival-----------------"<<endl;
//    for(size_t y=0;y<N;y++)
//    {
//        cout << arr[y].real()<<"  "<<arr[y].imag()<<endl;
//    }

    /*基于迭代的蝶形快速傅立叶变换，自底向上*/
     for(int s =1;s <= flag;s++)
    {
        int m = pow(2,(double)s);
        complex<double> wm = complex<double>(cos(2*pi/m),sin(2*pi/m));//wm1：从W21开始,周期变换
        for(int p = 0;p<N;p+=m)
        {
            complex<double> w(1,0);
            for(int j = 0;j<=m/2-1;j++)
            {
                complex<double> t = w * arr[p+j+m/2];
                complex<double> u = arr[p+j];
                arr[p+j] = u+t;
                arr[p+j+m/2] = u-t;
                w = w*wm;
            }
        }
    }
    return arr;
}
int countPadding(int n);
/*****************一维快速傅立叶逆变换********************/
/*fftimg:原始一维傅立叶序列
  n     : 序列长度
*******************************************************/
complex<double>* IFFT(const complex<double> *fftimg,int n)
{
    n = countPadding(n);
    double flag = log10((double)n)/log10(2.0);
    int N = n;
    if(flag - (int)flag != 0){   //判断是否为2的指数次
        cout <<"the length of srcimg is wrong"<< endl;
        /*填充成2的指数项*/
        cout <<"need padding"<<endl;
        N = pow(2,(double)((int)flag+1));
        flag = log10((double)N)/log10(2.0);
    }
    /*改变成奇偶顺序*/
    complex<double> * spit = new complex<double>[N];
    int sub=0;
    for(int k =0;k<N;k++)
    {
        sub =rev(k,(int)flag); 
        if(sub < n){
            spit[k] = complex<double>(*(fftimg + sub));
        }else{
            spit[k] = complex<double>(0,0);
        }
    }

    for(int s =1;s <= flag;s++)
    {
        int m = pow(2,(double)s);
        complex<double> wm = complex<double>(cos(-2*pi/m),sin(-2*pi/m));//wm1：从W2(-1)开始
        for(int p = 0;p<N;p+=m)
        {
            complex<double> w(1,0);
            for(int j = 0;j<=m/2-1;j++)
            {
                complex<double> t = w * spit[p+j+m/2];
                complex<double> u = spit[p+j];
                spit[p+j] = u+t;
                spit[p+j+m/2] = u-t;
                w = w*wm;
            }
        }
    }

    for(size_t p =0;p<n;p++)
    {
        spit[p] = spit[p]/complex<double>(N,0);
    }
    return spit;
}
/*******使用共轭的快速傅立叶逆变换**************/
complex<double>* IFFT2(const complex<double> *fftimg,int n)
{
    n = countPadding(n);
    complex<double> *gfftimg = new complex<double>[n];
    for(size_t i = 0;i<n;i++){
        gfftimg[i] = complex<double>(fftimg[i].real(),-fftimg[i].imag());
    }
    complex<double> *ifft = FFT(gfftimg,n); 
    for(size_t j = 0;j<n;j++)
    {
        ifft[j] = ifft[j]/complex<double>(n,0);
    }
    delete gfftimg;
    return ifft;
}

/*************二维快速傅立叶变换**************************
*srcimg:用一维表示的二维原始序列
*width ：宽度
height：高度
********************************************************/
complex<double>* twoDFFT(const double *srcimg,int width,int height)
{
    int w = countPadding(width);
    int h = countPadding(height);
    int pixes = w * h;
    complex<double> *hdirection = new complex<double>[w];
    complex<double> *vdirection = new complex<double>[h];
    complex<double> *fourier = new complex<double>[pixes];
    /********移位居中************/
    /*二维水平方向*/
    for(size_t i = 0;i<h;i++){
        for(size_t j = 0;j<w;j++){
            if(i>=height || j >=width){
                hdirection[j] = complex<double>(0,0);
            }else{
                hdirection[j] = complex<double>(srcimg[i*width + j],0);
            }
    //        cout << hdirection[j] << " ";
        }
    //    cout <<""<<endl;
        complex<double> *hfourier = FFT(hdirection,w);
        for(size_t m = 0;m<w;m++){
            fourier[i*w+m] = hfourier[m];
        }
        delete hfourier;
    }
    /*二维垂直方向*/
    for(size_t ii = 0;ii<w;ii++){
        for(size_t jj = 0;jj<h;jj++){
            vdirection[jj] = fourier[jj*w + ii];
        }
        complex<double> *vfourier = FFT(vdirection,h);
        for(size_t mm = 0;mm < h;mm++){
            fourier[mm*w +ii] = vfourier[mm];
        }
        delete vfourier;
    }
    delete hdirection;
    delete vdirection;
    return fourier;

}
/**************二维快速傅立叶逆变换*************************
*fourier : 一维表示的二维傅立叶变换序列
*宽
*height:高
********************************************/

complex<double>* twoDIFFT(const complex<double> *fourier,int width,int height)
{
    width = countPadding(width);
    height = countPadding(height);
    int fpoints = width * height;
    complex<double> *hdirection = new complex<double>[width];
    complex<double> *vdirection = new complex<double>[height];
    complex<double> *ifourier = new complex<double>[fpoints];

    for(size_t ii = 0;ii<height;ii++)
    {
        for(size_t jj = 0;jj<width;jj++){
            hdirection[jj] = fourier[ii*width+jj];
        }
        complex<double> *hifour = IFFT(hdirection,width);//临时变量
        for(size_t mm = 0;mm<width;mm++){
            ifourier[ii*width+mm] = hifour[mm];
        }
        delete hifour;
    }
    for(size_t i = 0;i<width;i++){
        for(size_t j = 0;j<height;j++){
            vdirection[j] = ifourier[j*width+i];
        }
        complex<double> *vifour = IFFT(vdirection,height);
        for(size_t m = 0;m<height;m++){
            ifourier[m*width+i] = vifour[m];
        }
        delete vifour;
    }    
    delete hdirection;
    delete vdirection;
    return ifourier;
}
/******计算填充数*********************/
int countPadding(int n)
{
    double lg = log10((double)n)/log10(2.0);
    if((lg - (int)lg) == 0){
        return n;
    }
    int N = pow(2.0,((int)lg+1));
    return N;
}
/*测试一维填充形式傅立叶变换*/
/*
int main()
{
    ifstream infile("D:\ftest.txt");
    if(!infile.is_open())
    {
        cout<<"file read error"<<endl;
        return -1;
    }
    istream_iterator<double> df(infile);
    istream_iterator<double> eof;
    int N = *df++;
    double *src = new double[N];
    size_t i = 0;
    while(df != eof)
    {
        src[i] = *df++;
        i++;
    }
/*    for(size_t a = 0;a < N;a++)
    {
        cout << src[a]<<endl;
    }
*/    
//    complex<double> *fourier = new complex<double>[N];
//    complex<double> *spit = new complex<double>[N];
/*    complex<double> *fourier = FFT(src,N);
    cout <<"--------fourier sieres----------"<<endl;
    for(size_t x = 0;x <16;x++)
    {
        cout << (*(fourier +x)).real()<<"  "<<(*(fourier+x)).imag() << endl;
    }
    cout <<"----------------inverse fourier -----------------------"<<endl;
    complex<double> *spit = IFFT(fourier,N);
    for(size_t y = 0;y <16;y++)
    {
        cout << spit[y].real()<<"  "<<spit[y].imag() << endl;
    }
    delete src;
//    delete fourier;
//    delete spit;
    return 1;
}*/
/*测试二维零填充傅立叶正负变换*/

/*
int main()
{
    ifstream infile("D:\ftest2.txt");
    if(!infile.is_open())
    {
        cout<<"file read error"<<endl;
        return -1;
    }
    istream_iterator<double> df(infile);
    istream_iterator<double> eof;
    int N = *df++;
    double *src = new double[N];
    size_t i = 0;
    while(df != eof)
    {
        src[i] = *df++;
        i++;
    }
    complex<double> *fourier = twoDFFT(src,3,6);
    cout <<"--------fourier sieres----------"<<endl;
    for(size_t x = 0;x <32;x++)
    {
        cout << (*(fourier +x)).real()<<"  "<<(*(fourier+x)).imag() << endl;
    }
    cout <<"----------------inverse fourier -----------------------"<<endl;
    complex<double> *spit = twoDIFFT(fourier,3,6);
    for(size_t y = 0;y <32;y++)
    {
        cout << spit[y].real()<<"  "<<spit[y].imag() << endl;
    }
    delete src;
//    delete fourier;
//    delete spit;
    return 1;
}*/

/*测试一维傅里叶图像变换*//*
int main(int argc,char ** argv)
{
    IplImage *img;
    img = cvLoadImage("F:\Fig0222(a)(face).tif",0);
    int dim = img->imageSize;
    //从图像矩阵复制出单维数组；
    double * src = new double[dim];
    size_t j =0;
    for(int y =0;y<img->height;y++)
    {
        uchar * ptr = (uchar*)(img->imageData + y * img->widthStep);
        for(int x =0;x<img->width;x++){
            src[j] = (double)ptr[x]/256;
            j++;
        }
    }
    //一维数组做傅立叶变换
    complex<double>* fourier = FFT(src,dim);
    //计算填充后傅氏数组大小
    double lg = log10((double)dim)/log10(2.0);
    int n = pow(2,(double)((int)lg+1));
    //傅立叶逆变换
    complex<double> *ifourier = IFFT(fourier,n);
    double ipix = 0;
    size_t po = 0;
    //重填充图像
    for(int y =0;y<img->height;y++)
    {
        uchar * ptr = (uchar*)(img->imageData + y * img->widthStep);
        for(int x =0;x<img->width;x++){
            ipix = ifourier[po].real();
            ptr[x] = (uchar)(ipix * 256);
            po++;
        }
    }
    cvNamedWindow("fourier_t",CV_WINDOW_AUTOSIZE);
    cvShowImage("fourier_t",img);
    cvWaitKey(0);
    cvReleaseImage(&img);
    cvDestroyWindow("fourier_t");

    return 1;
/*    ifstream infile("F:\Fig0222(a)(face).tif",ios::binary);
    vector<double> collection;
    istream_iterator<byte> in_iter(infile);
    istream_iterator<byte> eof;
    double temp = 0;
    unsigned char c = 0;
    while(in_iter != eof)
    {
        c = *in_iter++;
        temp = (double)c/256;
        collection.push_back(temp);
    }
    int amount = collection.size();
    double * src = new double[amount];  
    size_t j = 0;
    for(vector<double>::iterator iter = collection.begin();iter != collection.end();++iter)
    {
        src[j] = *iter;
        j++;
    }*/
/*    delete src;
    delete fourier;
    delete ifourier;
} 
*/

int main(int argc,char **argv)
{
    IplImage *img;
    if((img = cvLoadImage("F:\Fig0222(a)(face).tif",0))!=0){
        int dim = img->imageSize;
    //从图像矩阵复制出单维数组；
        double * src = new double[dim];
        size_t j =0;
        for(int y =0;y<img->height;y++)
        {
            uchar * ptr = (uchar*)(img->imageData + y * img->widthStep);
            for(int x =0;x<img->width;x++){
                src[j] = (double)ptr[x]*pow(-1,(double)(x+y))/256;
                j++;
            }
        }
        int w = img->width;
        int h = img->height;
        int pwid = countPadding(w);
        int phei = countPadding(h);
        complex<double> *twodfourier = twoDFFT(src,w,h);
        double *vals = new double[pwid*phei];
        char * pixval = new char[pwid*phei];
        CvMat freq;
        double max = 0;
        double min = 255;
        for(size_t p = 0;p<pwid*phei;p++){
            vals[p]=log(1+sqrt(pow(twodfourier[p].real(),2.0)+pow(twodfourier[p].imag(),2.0)));//对数级的幅度铺
            if(vals[p] > max){
                max = vals[p];
            }
            if(vals[p] < min){
                min = vals[p];
            }
        }
        cout<<min<< " "<<max<<endl;
        for(size_t s = 0;s<pwid*phei;s++){
            pixval[s] = (char)((vals[s]-min)/(max-min)*255);
        }
/*        for(size_t q = 0;q < pwid*phei;q++){
            cout <<vals[q]<<" "<<endl;
        }*/
        cvInitMatHeader(&freq,pwid,phei,CV_8UC1,pixval);
        IplImage *ipl = cvCreateImage(cvGetSize(&freq),8,1);
        cvGetImage(&freq,ipl);
        complex<double> *twodifourier = twoDIFFT(twodfourier,w,h);
        double ipix = 0;
        size_t po = 0;
        for(int y =0;y<img->height;y++)
        {
            uchar * ptr = (uchar*)(img->imageData + y * img->widthStep);
            for(int x =0;x<img->width;x++){
                ipix = twodifourier[po*pwid+x].real();
                ptr[x] = (uchar)(ipix * 256)*pow(-1,(double)(x+y));
            }
            po++;
        }
        cvNamedWindow("frequence domain",CV_WINDOW_AUTOSIZE);
        cvShowImage("frequence domain",ipl);
        cvNamedWindow("fourier_t",CV_WINDOW_AUTOSIZE);
        cvShowImage("fourier_t",img);
        cvWaitKey(0);
        cvReleaseImage(&ipl);
        cvReleaseImage(&img);
        cvDestroyWindow("fourier_t");
        cvDestroyWindow("frequence domain");
        delete vals;
        delete pixval;
        delete src;
        delete twodfourier;
        delete twodifourier;
        return 1;
    }
    return 0;
}

/*
int main()
{
    CImg<unsigned char> srcimg("F:\Fig0222(a)(face).tif");
    //CImgDisplay main_disp(srcimg,"lina");
    srcimg.display();
    return 1;

}
*/

 输出：

原图：

频谱图：