基于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 "opencv\highgui.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( ¢er_image, uchar, j, i);
CV_IMAGE_ELEM( ¢er_image, uchar, j, i) = CV_IMAGE_ELEM(
¢er_image, uchar, j+cy, i+cx);
CV_IMAGE_ELEM( ¢er_image, uchar, j+cy, i+cx) = tmp13;

tmp24 = CV_IMAGE_ELEM( ¢er_image, uchar, j, i+cx);
CV_IMAGE_ELEM( ¢er_image, uchar, j, i+cx) =
CV_IMAGE_ELEM( ¢er_image, uchar, j+cy, i);
CV_IMAGE_ELEM( ¢er_image, uchar, j+cy, i) = tmp24;
}
}

/*-----------------------------------------------------------------------------
*  用于保存FFT图像
*
*
*-----------------------------------------------------------------------------*/
Mat img(¢er_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