Fourier Transform

http://www.academictutorials.com/graphics/graphics-fourier-transform.asp


/* For description look into the help() function. */

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
#include <vector>
#include <complex>

using namespace std;
using namespace cv;

#define N 128

Mat fRe = Mat::zeros(1,N, CV_32FC1 );//the function's real part, imaginary part, and amplitude
Mat fIm = Mat::zeros(1,N, CV_32FC1 );
Mat fAmp = Mat::zeros(1,N, CV_32FC1 );
Mat FRe = Mat::zeros(1,N, CV_32FC1 ); //the FT's real part, imaginary part and amplitude
Mat FIm = Mat::zeros(1,N, CV_32FC1 );
Mat FAmp = Mat::zeros(1,N, CV_32FC1 );

//Calculates the DFT
void shiftDFT(int n,bool inverse, InputArray gRe, InputArray gIm, OutputArray GRe, OutputArray GIm)
{
        Mat tempgRe = gRe.getMat();
        Mat tempgIm = gIm.getMat();
        Mat tempGRe = GRe.getMat();
        Mat tempGIm = GIm.getMat();
        float a;
        for(int iw = 0;iw < FRe.rows;iw++)
        {
                for(int jw = 0;jw < FRe.cols;jw++)
                {
                        for(int i = 0;i < fRe.rows;i++)
                        {
                                for(int j = 0;j < fRe.cols;j++)
                                {
                                        a = (float)(-2 * CV_PI * jw * j / N);
                                        if(inverse) a = -a;
                                        double ca = cos(a);
                                        double sa = sin(a);
                                        tempGRe.at<float>(iw,jw) += (float)(tempgRe.at<float>(i,j) * ca - tempgIm.at<float>(i,j)* sa);
                                        tempGIm.at<float>(iw,jw) += (float)(tempgRe.at<float>(i,j) * sa + tempgIm.at<float>(i,j)* ca);
                                }
                        }
                }
        }

        //if(!inverse)
        //{
        // for(int iw = 0;iw < FRe.rows;iw++)
        // {
        // for(int jw = 0;jw < FRe.cols;jw++)
        // {
        // tempGRe.at<float>(iw,jw) /= N;
        // tempGIm.at<float>(iw,jw) /= N;
        // }
        // }
        //}

        if(!inverse)
        {
                tempGRe /= N;
                tempGIm /= N;
        }

        //shift
        Mat copyGRe = tempGRe.clone();
        Mat copyGIm = tempGIm.clone();
        for(int iw = 0;iw < FRe.rows;iw++)
        {
                for(int jw = 0;jw < FRe.cols;jw++)
                {
                        if(jw<N/2)
                        {
                                tempGRe.at<float>(iw,jw) = copyGRe.at<float>(iw,jw+N/2);
                                tempGIm.at<float>(iw,jw) = copyGIm.at<float>(iw,jw+N/2);
                        }
                        else
                        {
                                tempGRe.at<float>(iw,jw) = copyGRe.at<float>(iw,jw-N/2);
                                tempGIm.at<float>(iw,jw) = copyGIm.at<float>(iw,jw-N/2);
                        }
                }
        }
}

void FFT(int n,bool inverse, InputArray gRe, InputArray gIm, OutputArray GRe, OutputArray GIm)
{
        Mat tempgRe = gRe.getMat();
        Mat tempgIm = gIm.getMat();
        Mat tempGRe = GRe.getMat();
        Mat tempGIm = GIm.getMat();

        //Calculate m=log_2(n)
        int m = 0;
        int p = 1;
        while(p < n)
        {
                p *= 2;
                m++;
        }

        //Bit reversal
        tempGRe.at<float>(0,n - 1) = tempgRe.at<float>(0,n - 1);
        tempGIm.at<float>(0,n - 1) = tempgIm.at<float>(0,n - 1);

        int j = 0;
        for(int iw = 0; iw < n - 1; iw++)
        {
                tempGRe.at<float>(0,iw) = tempgRe.at<float>(0,j);
                tempGIm.at<float>(0,iw) = tempgIm.at<float>(0,j);
                int k = n / 2;
                while(k <= j)
                {
                        j -= k;
                        k /= 2;
                }
                j += k;
        }

        //Calculate the FFT
        float ca = -1.0;
        float sa = 0.0;
        int l1 = 1, l2 = 1;

        for(int l = 0; l < m; l++)
        {
                l1 = l2;
                l2 *= 2;
                float u1 = 1.0;
                float u2 = 0.0;
                for(int j = 0; j < l1; j++)
                {
                        for(int i = j; i < n; i += l2)
                        {
                                int i1 = i + l1;
                                float t1 = u1 * tempGRe.at<float>(0,i1) - u2 * tempGIm.at<float>(0,i1);
                                float t2 = u1 * tempGIm.at<float>(0,i1) + u2 * tempGRe.at<float>(0,i1);
                                tempGRe.at<float>(0,i1) = tempGRe.at<float>(0,i) - t1;
                                tempGIm.at<float>(0,i1) = tempGIm.at<float>(0,i) - t2;
                                tempGRe.at<float>(0,i) += t1;
                                tempGIm.at<float>(0,i) += t2;
                        }
                        float z = u1 * ca - u2 * sa;
                        u2 = u1 * sa + u2 * ca;
                        u1 = z;
                }
                sa = sqrt((float)((1.0 - ca) / 2.0));
                if(!inverse) sa =- sa;
                ca = sqrt((float)((1.0 + ca) / 2.0));
        }

         //Divide through n if it isn't the IDFT
        if(!inverse)
        {
                tempGRe /= N;
                tempGIm /= N;
        }
                
        //shift
        Mat copyGRe = tempGRe.clone();
        Mat copyGIm = tempGIm.clone();
        for(int iw = 0;iw < FRe.rows;iw++)
        {
                for(int jw = 0;jw < FRe.cols;jw++)
                {
                        if(jw<N/2)
                        {
                                tempGRe.at<float>(iw,jw) = copyGRe.at<float>(iw,jw+N/2);
                                tempGIm.at<float>(iw,jw) = copyGIm.at<float>(iw,jw+N/2);
                        }
                        else
                        {
                                tempGRe.at<float>(iw,jw) = copyGRe.at<float>(iw,jw-N/2);
                                tempGIm.at<float>(iw,jw) = copyGIm.at<float>(iw,jw-N/2);
                        }
                }
        }
}

void calculateAmp(OutputArray ga, InputArray gRe, InputArray gIm)
{
        Mat tempgRe = gRe.getMat();
        Mat tempgIm = gIm.getMat();
        Mat tempga = gIm.getMat();

        for(int i = 0;i < fRe.rows;i++)
        {
                for(int j = 0;j < fRe.cols;j++)
                {
                        tempga.at<float>(i,j) = sqrt(tempgRe.at<float>(i,j) * tempgRe.at<float>(i,j) + tempgIm.at<float>(i,j) * tempgIm.at<float>(i,j));
                }
        }
}

void upplot( OutputArray out, InputArray background ,InputArray f,CvScalar color)
{
        Mat tempbackground = background.getMat();
        Mat tempf = f.getMat();

        for(int i = 0;i < fRe.rows;i++)
        {
                for(int j = 0;j < fRe.cols;j++)
                {
                        line( tempbackground,Point( j*4, (int)(tempbackground.rows/4 - tempf.at<float>(i,j))),Point( j*4, (int)(tempbackground.rows/4 - tempf.at<float>(i,j))),color,2,8);
                }
        }
}

void downplot( OutputArray out, InputArray background ,InputArray f,CvScalar color)
{
        Mat tempbackground = background.getMat();
        Mat tempf = f.getMat();

        for(int i = 0;i < fRe.rows;i++)
        {
                for(int j = 0;j < fRe.cols;j++)
                {
                        line( tempbackground,Point( j*4, (int)(tempbackground.rows*3/4 - tempf.at<float>(i,j))),Point( j*4, (int)(tempbackground.rows*3/4 - tempf.at<float>(i,j))),color,3,8);
                }
        }
}
int main( )
{
         Load 图像
        //Mat src = imread( "lena.png");
        double t;
        Mat background = Mat::zeros( 128*4,128*4, CV_32FC3 );
        //rectangle( background,Point( 0,0),Point(background.cols,background.rows),Scalar( 255, 255, 255 ),-1,8);
        line( background,Point( 0,background.rows/4),Point(background.cols,background.rows/4),Scalar( 128, 128, 128 ),1,8);
        line( background,Point( 0,background.rows*3/4),Point(background.cols,background.rows*3/4),Scalar( 128, 128, 128 ),1,8);
        line( background,Point( background.cols/2,0),Point(background.cols/2,background.rows),Scalar( 128, 128, 128 ),1,8);
        
        for(int i = 0;i < fRe.rows;i++)
        {
                for(int j = 0;j < fRe.cols;j++)
                {
                        fRe.at<float>(i,j) = static_cast<float>(25.0 * sin(j / 2.0)); //generate a sine signal. Put anything else you like here.
                }
        }
        calculateAmp(fAmp, fRe, fIm);

        upplot( background , background ,fRe,Scalar( 0, 0, 255));
        upplot( background , background ,fIm,Scalar( 128, 255, 128));
        upplot( background , background ,fAmp,Scalar( 160, 160, 160));

        t = (double)getTickCount();
        shiftDFT(N,0,fRe, fIm, FRe, FIm);
        t = (double)getTickCount() - t;
        cout<<"shiftDFT in "<<t*1000/getTickFrequency()<<"ms"<<endl;
        calculateAmp(FAmp, FRe, FIm);
        downplot( background , background ,12.0*FRe,Scalar( 0, 0, 255));
        downplot( background , background ,12.0*FIm,Scalar( 128, 255, 128));
        downplot( background , background ,12.0*FAmp,Scalar( 160, 160, 160));

        namedWindow( "shiftDFT", CV_WINDOW_AUTOSIZE );
        imshow( "shiftDFT", background);


        t = (double)getTickCount();
        FFT(N,0,fRe, fIm, FRe, FIm);
        t = (double)getTickCount() - t;
        cout<<"shiftFFT in "<<t*1000/getTickFrequency()<<"ms"<<endl;
        calculateAmp(FAmp, FRe, FIm);
        
        rectangle( background,Point(0,background.rows/3),Point(background.cols,background.rows),Scalar( 0, 0, 0 ),-1,8);
        line( background,Point( 0,background.rows*3/4),Point(background.cols,background.rows*3/4),Scalar( 128, 128, 128 ),1,8);
        line( background,Point( background.cols/2,0),Point(background.cols/2,background.rows),Scalar( 128, 128, 128 ),1,8);
        downplot( background , background ,12.0*FRe,Scalar( 0, 0, 255));
        downplot( background , background ,12.0*FIm,Scalar( 128, 255, 128));
        downplot( background , background ,12.0*FAmp,Scalar( 160, 160, 160));
        
        namedWindow( "shiftFFT", CV_WINDOW_AUTOSIZE );
        imshow( "shiftFFT", background);

        
        waitKey(0);
        return 0;

}