基于MFC和opencv的FFT

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

     /*-----------------------------------------------------------------------------
      *  计算功率谱
      *    和归一化
      *
      *-----------------------------------------------------------------------------*/
     // 行
     for(i = ; i < h; i++)
     {
         // 列
         for(j = ; 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()) / ;

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

             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:  [e-mail]
  *             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()), );

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

         // 绘制图标
         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(),);

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

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

     while(h *  <= image.rows)
     {
         h *= ;
         hp++;
     }
     // 分配内存
     complex<double> *TD = new complex<double>[w * h];
     complex<double> *FD = new complex<double>[w * h];
     // 行
     for(i = ; i < h; i++)
     {
         // 列
         for(j = ; j < w; j++)
         {
             // 给时域赋值
             TD[j + w * i] = complex<double>(image.at<uchar>(i,j), ); /* opencv函数读取图像像素到复数矩阵里 */
         }
     }
     /*-----------------------------------------------------------------------------
      *  把二维的FFT换成分别对行方向和列方向进行一维的FFT
      *
      *
      *-----------------------------------------------------------------------------*/
     for(i = ; i < h; i++)
     {
         // 对y方向进行快速付立叶变换
         FFT(&TD[w * i], &FD[w * i], wp);
     }
     // 保存变换结果
     for(i = ; i < h; i++)
     {
         for(j = ; j < w; j++)
         {
             TD[i + h * j] = FD[j + w * i];
         }
     }

     for(i = ; i < w; i++)
     {
         // 对x方向进行快速付立叶变换
         FFT(&TD[i * h], &FD[i * h], hp);
     }

     /*-----------------------------------------------------------------------------
      *  计算功率谱
      *    和归一化
      *
      *-----------------------------------------------------------------------------*/
     // 行
     for(i = ; i < h; i++)
     {
         // 列
         for(j = ; 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()) / ;

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

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

         }
     }

     /*-----------------------------------------------------------------------------
      *  释放内存
      *
      *
      *-----------------------------------------------------------------------------*/
     // 删除临时变量
     delete TD;
     delete FD;

     /*-----------------------------------------------------------------------------
      *  进行图像的中心化
      *
      *
      *-----------------------------------------------------------------------------*/
     int cy = image.rows/;                      /* 中心点位置 cx ,cy */
     int cx = image.cols/;
     uchar tmp13,tmp24;

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

     IplImage  center_image=IplImage(image);
     for( j = ; j < cy; j++ ){
         for( i = ; 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,);
     if (BST_CHECKED == ::IsDlgButtonChecked(m_hWnd,IDC_CHECK_SAVE))
     {

         strtitle="FFT_"+strtitle+"."+strext;

         //MessageBox(strtitle.GetBuffer(0));
         imwrite(strtitle.GetBuffer(),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 =  << r;
     // 分配运算所需存储器
     W  = new complex<double>[count / ];
     X1 = new complex<double>[count];
     X2 = new complex<double>[count];
     // 计算加权系数
     for(i = ; i < count / ; i++)
     {
         angle =-*i*PI;
         angle =angle/(double)count;
         W[i] = complex<double> (cos(angle), sin(angle));
     }
     // 将时域点写入X1
     memcpy(X1, TD, sizeof(complex<double>) * count);

     /*-----------------------------------------------------------------------------
      *
      *                   采用蝶形算法进行快速付立叶变换
      *
      *-----------------------------------------------------------------------------*/

     for(k = ; k < r; k++)
     {
         for(j = ; j <  << k; j++)
         {
             bfsize =  << (r-k);
             for(i = ; i < bfsize / ; i++)
             {
                 p = j * bfsize;
                 X2[i + p] = X1[i + p] + X1[i + p + bfsize / ];
                 X2[i + p + bfsize / ] = (X1[i + p] - X1[i + p + bfsize / ]) * W[i * (<<k)];
             }
         }
         X  = X1;
         X1 = X2;
         X2 = X;
     }
     // 重新排序
     for(j = ; j < count; j++)
     {
         p = ;
         for(i = ; i < r; i++)
         {
             if (j&(<<i))
             {
                 p+=<<(r-i-);
             }
         }
         FD[j]=X1[p];
     }
     // 释放内存
     delete W;
     delete X1;
     delete X2;
 }

运行效果如下：

提供程序一份：

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