基于MFC和opencv的FFT

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

 1     /*-----------------------------------------------------------------------------
 2      *  计算功率谱
 3      *    和归一化
 4      *
 5      *-----------------------------------------------------------------------------*/
 6     // 行
 7     for(i = 0; i < h; i++)
 8     {
 9         // 列
10         for(j = 0; j < w; j++)
11         {
12             // 计算频谱
13             dTemp = sqrt(FD[j * h + i].real() * FD[j * h + i].real() + 
14                 FD[j * h + i].imag() * FD[j * h + i].imag()) / 100;
15 
16             // 判断是否超过255
17             if (dTemp > 255)
18             {
19                 // 对于超过的，直接设置为255
20                 dTemp = 255;
21             }
22 
23             image.at<uchar>(i,j)=(uchar)dTemp;
24             // 更新源图像
25 
26         }
27     }

主要代码：

  1 /*
  2  * =====================================================================================
  3  *
  4  *       Filename:  fft_dlgDlg.cpp
  5  *      Environment:opencv2.4.4 vs2010
  6  *            
  7  *    Description:  基于MFC的FFT程序
  8  *
  9  *
 10  *
 11  *        Version:  1.0
 12  *        Created:  2013/10/19 19:24:06
 13  *         Author:  yuliyang
 14 I*
 15  *             Mail:  wzyuliyang911@gmail.com
 16  *             Blog:  http://www.cnblogs.com/yuliyang
 17  *
 18  * =====================================================================================
 19  */
 20 
 21 
 22 // fft_dlgDlg.cpp : 实现文件
 23 //
 24 
 25 #include "stdafx.h"
 26 #include "fft_dlg.h"
 27 #include "fft_dlgDlg.h"
 28 #include "afxdialogex.h"
 29 
 30 #include <complex>
 31 #include <math.h>
 32 #include <stdio.h>
 33 #include <Windows.h>
 34 #include "opencv\highgui.h"
 35 #include "opencv2/core/core.hpp"  
 36 #include "opencv2/highgui/highgui.hpp"  
 37 #include "opencv2/imgproc/imgproc.hpp"  
 38 using namespace cv;
 39 using namespace std;
 40 #define PI 3.1415936;
 41 
 42 #ifdef _DEBUG
 43 #define new DEBUG_NEW
 44 #endif
 45 
 46 
 47 // 用于应用程序“关于”菜单项的 CAboutDlg 对话框
 48 
 49 class CAboutDlg : public CDialogEx
 50 {
 51 public:
 52     CAboutDlg();
 53 
 54 // 对话框数据
 55     enum { IDD = IDD_ABOUTBOX };
 56 
 57     protected:
 58     virtual void DoDataExchange(CDataExchange* pDX);    // DDX/DDV 支持
 59 
 60 // 实现
 61 protected:
 62     DECLARE_MESSAGE_MAP()
 63 };
 64 
 65 CAboutDlg::CAboutDlg() : CDialogEx(CAboutDlg::IDD)
 66 {
 67 }
 68 
 69 void CAboutDlg::DoDataExchange(CDataExchange* pDX)
 70 {
 71     CDialogEx::DoDataExchange(pDX);
 72 }
 73 
 74 BEGIN_MESSAGE_MAP(CAboutDlg, CDialogEx)
 75 END_MESSAGE_MAP()
 76 
 77 
 78 // Cfft_dlgDlg 对话框
 79 
 80 
 81 
 82 
 83 Cfft_dlgDlg::Cfft_dlgDlg(CWnd* pParent /*=NULL*/)
 84     : CDialogEx(Cfft_dlgDlg::IDD, pParent)
 85 {
 86     m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);
 87 }
 88 
 89 void Cfft_dlgDlg::DoDataExchange(CDataExchange* pDX)
 90 {
 91     CDialogEx::DoDataExchange(pDX);
 92 }
 93 
 94 BEGIN_MESSAGE_MAP(Cfft_dlgDlg, CDialogEx)
 95     ON_WM_SYSCOMMAND()
 96     ON_WM_PAINT()
 97     ON_WM_QUERYDRAGICON()
 98     ON_BN_CLICKED(IDC_OPEN_FILE, &Cfft_dlgDlg::OnBnClickedOpenFile)
 99     
100 END_MESSAGE_MAP()
101 
102 
103 // Cfft_dlgDlg 消息处理程序
104 
105 BOOL Cfft_dlgDlg::OnInitDialog()
106 {
107     CDialogEx::OnInitDialog();
108 
109     // 将“关于...”菜单项添加到系统菜单中。
110 
111     // IDM_ABOUTBOX 必须在系统命令范围内。
112     ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX);
113     ASSERT(IDM_ABOUTBOX < 0xF000);
114 
115     CMenu* pSysMenu = GetSystemMenu(FALSE);
116     if (pSysMenu != NULL)
117     {
118         BOOL bNameValid;
119         CString strAboutMenu;
120         bNameValid = strAboutMenu.LoadString(IDS_ABOUTBOX);
121         ASSERT(bNameValid);
122         if (!strAboutMenu.IsEmpty())
123         {
124             pSysMenu->AppendMenu(MF_SEPARATOR);
125             pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);
126         }
127     }
128 
129     // 设置此对话框的图标。当应用程序主窗口不是对话框时，框架将自动
130     //  执行此操作
131     SetIcon(m_hIcon, TRUE);            // 设置大图标
132     SetIcon(m_hIcon, FALSE);        // 设置小图标
133 
134     // TODO: 在此添加额外的初始化代码
135 
136     return TRUE;  // 除非将焦点设置到控件，否则返回 TRUE
137 }
138 
139 void Cfft_dlgDlg::OnSysCommand(UINT nID, LPARAM lParam)
140 {
141     if ((nID & 0xFFF0) == IDM_ABOUTBOX)
142     {
143         CAboutDlg dlgAbout;
144         dlgAbout.DoModal();
145     }
146     else
147     {
148         CDialogEx::OnSysCommand(nID, lParam);
149     }
150 }
151 
152 // 如果向对话框添加最小化按钮，则需要下面的代码
153 //  来绘制该图标。对于使用文档/视图模型的 MFC 应用程序，
154 //  这将由框架自动完成。
155 
156 void Cfft_dlgDlg::OnPaint()
157 {
158     if (IsIconic())
159     {
160         CPaintDC dc(this); // 用于绘制的设备上下文
161 
162         SendMessage(WM_ICONERASEBKGND, reinterpret_cast<WPARAM>(dc.GetSafeHdc()), 0);
163 
164         // 使图标在工作区矩形中居中
165         int cxIcon = GetSystemMetrics(SM_CXICON);
166         int cyIcon = GetSystemMetrics(SM_CYICON);
167         CRect rect;
168         GetClientRect(&rect);
169         int x = (rect.Width() - cxIcon + 1) / 2;
170         int y = (rect.Height() - cyIcon + 1) / 2;
171 
172         // 绘制图标
173         dc.DrawIcon(x, y, m_hIcon);
174     }
175     else
176     {
177         CDialogEx::OnPaint();
178     }
179 }
180 
181 //当用户拖动最小化窗口时系统调用此函数取得光标
182 //显示。
183 HCURSOR Cfft_dlgDlg::OnQueryDragIcon()
184 {
185     return static_cast<HCURSOR>(m_hIcon);
186 }
187 
188 
189 
190 void Cfft_dlgDlg::OnBnClickedOpenFile()
191 {
192     // TODO: 在此添加控件通知处理程序代码
193     
194     /*-----------------------------------------------------------------------------
195      *  选择文件名
196      *-----------------------------------------------------------------------------*/
197     
198     CString strFileName,strszFilter,strtitle,strext;
199     strszFilter="位图文件(*.bmp)|*.bmp|位图文件(*.jpg)|*.jpg|全部文件(*.*)|*.*||";
200     CFileDialog bmpdlg(TRUE,NULL,NULL,OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT,strszFilter,NULL);
201     if(IDOK == bmpdlg.DoModal())
202     {
203         strFileName = bmpdlg.GetPathName();
204         strtitle=bmpdlg.GetFileTitle();
205         strext=bmpdlg.GetFileExt();
206 
207     }
208     if (strFileName.IsEmpty())
209     {
210         //MessageBox((LPCTSTR)strFileName);
211         MessageBox("请选择一副图像");
212 
213         return ;
214     }
215 
216     /*-----------------------------------------------------------------------------
217      *  opencv读入图像
218      *
219      *
220      *-----------------------------------------------------------------------------*/
221     Mat image=imread(strFileName.GetBuffer(0),0);
222 
223     //unsigned char*    lpSrc;
224     // 中间变量
225     double    dTemp;
226     // 循环变量
227     LONG    i;
228     LONG    j;
229     // 进行付立叶变换的宽度和高度（2的整数次方）
230     LONG    w;
231     LONG    h;
232     int        wp;
233     int        hp;
234     // 赋初值
235     w = 1;
236     h = 1;
237     wp = 0;
238     hp = 0;
239 
240     /*-----------------------------------------------------------------------------
241      *  填充到2的幂次方
242      *
243      *
244      *-----------------------------------------------------------------------------*/
245     // 计算进行付立叶变换的宽度和高度（2的整数次方）
246     while(w * 2 <= image.cols)
247     {
248         w *= 2;
249         wp++;
250     }
251 
252     while(h * 2 <= image.rows)
253     {
254         h *= 2;
255         hp++;
256     }
257     // 分配内存
258     complex<double> *TD = new complex<double>[w * h];
259     complex<double> *FD = new complex<double>[w * h];
260     // 行
261     for(i = 0; i < h; i++)
262     {
263         // 列
264         for(j = 0; j < w; j++)
265         {
266             // 给时域赋值
267             TD[j + w * i] = complex<double>(image.at<uchar>(i,j), 0); /* opencv函数读取图像像素到复数矩阵里 */
268         }
269     }
270     /*-----------------------------------------------------------------------------
271      *  把二维的FFT换成分别对行方向和列方向进行一维的FFT
272      *
273      *
274      *-----------------------------------------------------------------------------*/
275     for(i = 0; i < h; i++)
276     {
277         // 对y方向进行快速付立叶变换
278         FFT(&TD[w * i], &FD[w * i], wp);
279     }
280     // 保存变换结果
281     for(i = 0; i < h; i++)
282     {
283         for(j = 0; j < w; j++)
284         {
285             TD[i + h * j] = FD[j + w * i];
286         }
287     }
288 
289     for(i = 0; i < w; i++)
290     {
291         // 对x方向进行快速付立叶变换
292         FFT(&TD[i * h], &FD[i * h], hp);
293     }
294 
295     /*-----------------------------------------------------------------------------
296      *  计算功率谱
297      *    和归一化
298      *
299      *-----------------------------------------------------------------------------*/
300     // 行
301     for(i = 0; i < h; i++)
302     {
303         // 列
304         for(j = 0; j < w; j++)
305         {
306             // 计算频谱
307             dTemp = sqrt(FD[j * h + i].real() * FD[j * h + i].real() + 
308                 FD[j * h + i].imag() * FD[j * h + i].imag()) / 100;
309 
310             // 判断是否超过255
311             if (dTemp > 255)
312             {
313                 // 对于超过的，直接设置为255
314                 dTemp = 255;
315             }
316 
317             image.at<uchar>(i,j)=(uchar)dTemp;
318             // 更新源图像
319 
320         }
321     }
322 
323     /*-----------------------------------------------------------------------------
324      *  释放内存
325      *
326      *
327      *-----------------------------------------------------------------------------*/
328     // 删除临时变量
329     delete TD;
330     delete FD;
331     
332     /*-----------------------------------------------------------------------------
333      *  进行图像的中心化
334      *
335      *
336      *-----------------------------------------------------------------------------*/
337     int cy = image.rows/2;                      /* 中心点位置 cx ,cy */
338     int cx = image.cols/2;
339     uchar tmp13,tmp24;
340 
341     //imshow("未中心化的功率谱",image);
342 
343     IplImage  center_image=IplImage(image);
344     for( j = 0; j < cy; j++ ){  
345         for( i = 0; i < cx; i++ ){  
346             //中心化，将整体份成四块进行对角交换  
347             tmp13 = CV_IMAGE_ELEM( &center_image, uchar, j, i);  
348             CV_IMAGE_ELEM( &center_image, uchar, j, i) = CV_IMAGE_ELEM(  
349                 &center_image, uchar, j+cy, i+cx);  
350             CV_IMAGE_ELEM( &center_image, uchar, j+cy, i+cx) = tmp13;  
351 
352             tmp24 = CV_IMAGE_ELEM( &center_image, uchar, j, i+cx);  
353             CV_IMAGE_ELEM( &center_image, uchar, j, i+cx) =  
354                 CV_IMAGE_ELEM( &center_image, uchar, j+cy, i);  
355             CV_IMAGE_ELEM( &center_image, uchar, j+cy, i) = tmp24;  
356         }  
357     }  
358 
359 
360     /*-----------------------------------------------------------------------------
361      *  用于保存FFT图像
362      *
363      *
364      *-----------------------------------------------------------------------------*/
365     Mat img(&center_image,0);
366     if (BST_CHECKED == ::IsDlgButtonChecked(m_hWnd,IDC_CHECK_SAVE))
367     {
368 
369         strtitle="FFT_"+strtitle+"."+strext;
370         
371         //MessageBox(strtitle.GetBuffer(0));
372         imwrite(strtitle.GetBuffer(0),img);
373 
374     }
375     imshow("中心化后的功率谱",img);
376 
377     waitKey();
378     // 返回
379     //return TRUE;
380     
381 
382 
383 }
384 
385 /* 
386  * ===  FUNCTION  ======================================================================
387  *         Name:  FFT
388  *  Description:  计算一维FFT
389  * =====================================================================================
390  */
391 void Cfft_dlgDlg::FFT(complex<double> * TD, complex<double> * FD, int r)
392 {
393 
394     LONG    count;
395     // 循环变量
396     int        i,j,k;
397     // 中间变量
398     int        bfsize,p;
399     // 角度
400     double    angle;
401     complex<double> *W,*X1,*X2,*X;
402     // 计算付立叶变换点数
403     count = 1 << r;
404     // 分配运算所需存储器
405     W  = new complex<double>[count / 2];
406     X1 = new complex<double>[count];
407     X2 = new complex<double>[count];
408     // 计算加权系数
409     for(i = 0; i < count / 2; i++)
410     {
411         angle =-2*i*PI;
412         angle =angle/(double)count;
413         W[i] = complex<double> (cos(angle), sin(angle));
414     }
415     // 将时域点写入X1
416     memcpy(X1, TD, sizeof(complex<double>) * count);
417 
418     /*-----------------------------------------------------------------------------
419      *  
420      *                   采用蝶形算法进行快速付立叶变换
421      *
422      *-----------------------------------------------------------------------------*/
423      
424     for(k = 0; k < r; k++)
425     {
426         for(j = 0; j < 1 << k; j++)
427         {
428             bfsize = 1 << (r-k);
429             for(i = 0; i < bfsize / 2; i++)
430             {
431                 p = j * bfsize;
432                 X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
433                 X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];
434             }
435         }
436         X  = X1;
437         X1 = X2;
438         X2 = X;
439     }
440     // 重新排序
441     for(j = 0; j < count; j++)
442     {
443         p = 0;
444         for(i = 0; i < r; i++)
445         {
446             if (j&(1<<i))
447             {
448                 p+=1<<(r-i-1);
449             }
450         }
451         FD[j]=X1[p];
452     }
453     // 释放内存
454     delete W;
455     delete X1;
456     delete X2;
457 }

运行效果如下：

 

 

提供程序一份：

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

 

posted @ 2013-10-19 19:52  小菜鸟_yang  阅读(2531)  评论(1)  编辑  收藏  举报