基于MFC和opencv的FFT

在网上折腾了一阵子，终于把这个程序写好了，程序是基于MFC的，图像显示的部分和获取图像的像素点是用到了opencv的一些函数，不过FFT算法没有用opencv的（呵呵，老师不让），网上的二维的FFT程序一般都是把图像分别进行行变换后进行列变换的，在编程过程中遇到了一些问题，是这样的，FFT算法算完后得到的复数矩阵怎么imshow?问题就出现在这，我原来的程序因为归一化到0-255时，程序运行特别慢（用了个CArray，找出array里的最大值和最小值，然后（每一个复数矩阵求模后-最小值）/(最大值-最小值)，不满才怪呵呵，得出FFT结果是全黑的）。参考了别人的归一化

    /*-----------------------------------------------------------------------------
     *  计算功率谱
     *    和归一化
     *
     *-----------------------------------------------------------------------------*/
    // 行
    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;

            // 判断是否超过255
            if (dTemp > 255)
            {
                // 对于超过的，直接设置为255
                dTemp = 255;
            }

            image.at<uchar>(i,j)=(uchar)dTemp;
            // 更新源图像

        }
    }

主要代码：

/*
 * =====================================================================================
 *
 *       Filename:  fft_dlgDlg.cpp
 *      Environment:opencv2.4.4 vs2010
 *            
 *    Description:  基于MFC的FFT程序
 *
 *
 *
 *        Version:  1.0
 *        Created:  2013/10/19 19:24:06
 *         Author:  yuliyang
I*
 *             Mail:  wzyuliyang911@gmail.com
 *             Blog:  http://www.cnblogs.com/yuliyang
 *
 * =====================================================================================
 */


// fft_dlgDlg.cpp : 实现文件
//

#include "stdafx.h"
#include "fft_dlg.h"
#include "fft_dlgDlg.h"
#include "afxdialogex.h"

#include <complex>
#include <math.h>
#include <stdio.h>
#include <Windows.h>
#include "opencvhighgui.h"
#include "opencv2/core/core.hpp"  
#include "opencv2/highgui/highgui.hpp"  
#include "opencv2/imgproc/imgproc.hpp"  
using namespace cv;
using namespace std;
#define PI 3.1415936;

#ifdef _DEBUG
#define new DEBUG_NEW
#endif


// 用于应用程序“关于”菜单项的 CAboutDlg 对话框

class CAboutDlg : public CDialogEx
{
public:
    CAboutDlg();

// 对话框数据
    enum { IDD = IDD_ABOUTBOX };

    protected:
    virtual void DoDataExchange(CDataExchange* pDX);    // DDX/DDV 支持

// 实现
protected:
    DECLARE_MESSAGE_MAP()
};

CAboutDlg::CAboutDlg() : CDialogEx(CAboutDlg::IDD)
{
}

void CAboutDlg::DoDataExchange(CDataExchange* pDX)
{
    CDialogEx::DoDataExchange(pDX);
}

BEGIN_MESSAGE_MAP(CAboutDlg, CDialogEx)
END_MESSAGE_MAP()


// Cfft_dlgDlg 对话框




Cfft_dlgDlg::Cfft_dlgDlg(CWnd* pParent /*=NULL*/)
    : CDialogEx(Cfft_dlgDlg::IDD, pParent)
{
    m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
}

void Cfft_dlgDlg::DoDataExchange(CDataExchange* pDX)
{
    CDialogEx::DoDataExchange(pDX);
}

BEGIN_MESSAGE_MAP(Cfft_dlgDlg, CDialogEx)
    ON_WM_SYSCOMMAND()
    ON_WM_PAINT()
    ON_WM_QUERYDRAGICON()
    ON_BN_CLICKED(IDC_OPEN_FILE, &Cfft_dlgDlg::OnBnClickedOpenFile)
    
END_MESSAGE_MAP()


// Cfft_dlgDlg 消息处理程序

BOOL Cfft_dlgDlg::OnInitDialog()
{
    CDialogEx::OnInitDialog();

    // 将“关于...”菜单项添加到系统菜单中。

    // IDM_ABOUTBOX 必须在系统命令范围内。
    ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX);
    ASSERT(IDM_ABOUTBOX < 0xF000);

    CMenu* pSysMenu = GetSystemMenu(FALSE);
    if (pSysMenu != NULL)
    {
        BOOL bNameValid;
        CString strAboutMenu;
        bNameValid = strAboutMenu.LoadString(IDS_ABOUTBOX);
        ASSERT(bNameValid);
        if (!strAboutMenu.IsEmpty())
        {
            pSysMenu->AppendMenu(MF_SEPARATOR);
            pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);
        }
    }

    // 设置此对话框的图标。当应用程序主窗口不是对话框时，框架将自动
    //  执行此操作
    SetIcon(m_hIcon, TRUE);            // 设置大图标
    SetIcon(m_hIcon, FALSE);        // 设置小图标

    // TODO: 在此添加额外的初始化代码

    return TRUE;  // 除非将焦点设置到控件，否则返回 TRUE
}

void Cfft_dlgDlg::OnSysCommand(UINT nID, LPARAM lParam)
{
    if ((nID & 0xFFF0) == IDM_ABOUTBOX)
    {
        CAboutDlg dlgAbout;
        dlgAbout.DoModal();
    }
    else
    {
        CDialogEx::OnSysCommand(nID, lParam);
    }
}

// 如果向对话框添加最小化按钮，则需要下面的代码
//  来绘制该图标。对于使用文档/视图模型的 MFC 应用程序，
//  这将由框架自动完成。

void Cfft_dlgDlg::OnPaint()
{
    if (IsIconic())
    {
        CPaintDC dc(this); // 用于绘制的设备上下文

        SendMessage(WM_ICONERASEBKGND, reinterpret_cast<WPARAM>(dc.GetSafeHdc()), 0);

        // 使图标在工作区矩形中居中
        int cxIcon = GetSystemMetrics(SM_CXICON);
        int cyIcon = GetSystemMetrics(SM_CYICON);
        CRect rect;
        GetClientRect(&rect);
        int x = (rect.Width() - cxIcon + 1) / 2;
        int y = (rect.Height() - cyIcon + 1) / 2;

        // 绘制图标
        dc.DrawIcon(x, y, m_hIcon);
    }
    else
    {
        CDialogEx::OnPaint();
    }
}

//当用户拖动最小化窗口时系统调用此函数取得光标
//显示。
HCURSOR Cfft_dlgDlg::OnQueryDragIcon()
{
    return static_cast<HCURSOR>(m_hIcon);
}



void Cfft_dlgDlg::OnBnClickedOpenFile()
{
    // TODO: 在此添加控件通知处理程序代码
    
    /*-----------------------------------------------------------------------------
     *  选择文件名
     *-----------------------------------------------------------------------------*/
    
    CString strFileName,strszFilter,strtitle,strext;
    strszFilter="位图文件(*.bmp)|*.bmp|位图文件(*.jpg)|*.jpg|全部文件(*.*)|*.*||";
    CFileDialog bmpdlg(TRUE,NULL,NULL,OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT,strszFilter,NULL);
    if(IDOK == bmpdlg.DoModal())
    {
        strFileName = bmpdlg.GetPathName();
        strtitle=bmpdlg.GetFileTitle();
        strext=bmpdlg.GetFileExt();

    }
    if (strFileName.IsEmpty())
    {
        //MessageBox((LPCTSTR)strFileName);
        MessageBox("请选择一副图像");

        return ;
    }

    /*-----------------------------------------------------------------------------
     *  opencv读入图像
     *
     *
     *-----------------------------------------------------------------------------*/
    Mat image=imread(strFileName.GetBuffer(0),0);

    //unsigned char*    lpSrc;
    // 中间变量
    double    dTemp;
    // 循环变量
    LONG    i;
    LONG    j;
    // 进行付立叶变换的宽度和高度（2的整数次方）
    LONG    w;
    LONG    h;
    int        wp;
    int        hp;
    // 赋初值
    w = 1;
    h = 1;
    wp = 0;
    hp = 0;

    /*-----------------------------------------------------------------------------
     *  填充到2的幂次方
     *
     *
     *-----------------------------------------------------------------------------*/
    // 计算进行付立叶变换的宽度和高度（2的整数次方）
    while(w * 2 <= image.cols)
    {
        w *= 2;
        wp++;
    }

    while(h * 2 <= image.rows)
    {
        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++)
        {
            // 给时域赋值
            TD[j + w * i] = complex<double>(image.at<uchar>(i,j), 0); /* opencv函数读取图像像素到复数矩阵里 */
        }
    }
    /*-----------------------------------------------------------------------------
     *  把二维的FFT换成分别对行方向和列方向进行一维的FFT
     *
     *
     *-----------------------------------------------------------------------------*/
    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];
        }
    }

    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;

            // 判断是否超过255
            if (dTemp > 255)
            {
                // 对于超过的，直接设置为255
                dTemp = 255;
            }

            image.at<uchar>(i,j)=(uchar)dTemp;
            // 更新源图像

        }
    }

    /*-----------------------------------------------------------------------------
     *  释放内存
     *
     *
     *-----------------------------------------------------------------------------*/
    // 删除临时变量
    delete TD;
    delete FD;
    
    /*-----------------------------------------------------------------------------
     *  进行图像的中心化
     *
     *
     *-----------------------------------------------------------------------------*/
    int cy = image.rows/2;                      /* 中心点位置 cx ,cy */
    int cx = image.cols/2;
    uchar tmp13,tmp24;

    //imshow("未中心化的功率谱",image);

    IplImage  center_image=IplImage(image);
    for( j = 0; j < cy; j++ ){  
        for( i = 0; i < cx; i++ ){  
            //中心化，将整体份成四块进行对角交换  
            tmp13 = CV_IMAGE_ELEM( &center_image, uchar, j, i);  
            CV_IMAGE_ELEM( &center_image, uchar, j, i) = CV_IMAGE_ELEM(  
                &center_image, uchar, j+cy, i+cx);  
            CV_IMAGE_ELEM( &center_image, uchar, j+cy, i+cx) = tmp13;  

            tmp24 = CV_IMAGE_ELEM( &center_image, uchar, j, i+cx);  
            CV_IMAGE_ELEM( &center_image, uchar, j, i+cx) =  
                CV_IMAGE_ELEM( &center_image, uchar, j+cy, i);  
            CV_IMAGE_ELEM( &center_image, uchar, j+cy, i) = tmp24;  
        }  
    }  


    /*-----------------------------------------------------------------------------
     *  用于保存FFT图像
     *
     *
     *-----------------------------------------------------------------------------*/
    Mat img(&center_image,0);
    if (BST_CHECKED == ::IsDlgButtonChecked(m_hWnd,IDC_CHECK_SAVE))
    {

        strtitle="FFT_"+strtitle+"."+strext;
        
        //MessageBox(strtitle.GetBuffer(0));
        imwrite(strtitle.GetBuffer(0),img);

    }
    imshow("中心化后的功率谱",img);

    waitKey();
    // 返回
    //return TRUE;
    


}

/* 
 * ===  FUNCTION  ======================================================================
 *         Name:  FFT
 *  Description:  计算一维FFT
 * =====================================================================================
 */
void Cfft_dlgDlg::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 =-2*i*PI;
        angle =angle/(double)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;
}

运行效果如下：

提供程序一份：

http://pan.baidu.com/s/1ehvwy