在网上折腾了一阵子，终于把这个程序写好了，程序是基于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 // 判断是否超过25517 if (dTemp > 255)18 {19 // 对于超过的，直接设置为25520 dTemp = 255;21 }22 23 image.at (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 #include31 #include 32 #include 33 #include 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; // 除非将焦点设置到控件，否则返回 TRUE137 }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 else147 {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 (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 else176 {177 CDialogEx::OnPaint();178 }179 }180 181 //当用户拖动最小化窗口时系统调用此函数取得光标182 //显示。183 HCURSOR Cfft_dlgDlg::OnQueryDragIcon()184 {185 return static_cast (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 *TD = new complex [w * h];259 complex *FD = new complex [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 (image.at (i,j), 0); /* opencv函数读取图像像素到复数矩阵里 */268 }269 }270 /*-----------------------------------------------------------------------------271 * 把二维的FFT换成分别对行方向和列方向进行一维的FFT272 *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 // 判断是否超过255311 if (dTemp > 255)312 {313 // 对于超过的，直接设置为255314 dTemp = 255;315 }316 317 image.at (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( ¢er_image, uchar, j, i); 348 CV_IMAGE_ELEM( ¢er_image, uchar, j, i) = CV_IMAGE_ELEM( 349 ¢er_image, uchar, j+cy, i+cx); 350 CV_IMAGE_ELEM( ¢er_image, uchar, j+cy, i+cx) = tmp13; 351 352 tmp24 = CV_IMAGE_ELEM( ¢er_image, uchar, j, i+cx); 353 CV_IMAGE_ELEM( ¢er_image, uchar, j, i+cx) = 354 CV_IMAGE_ELEM( ¢er_image, uchar, j+cy, i); 355 CV_IMAGE_ELEM( ¢er_image, uchar, j+cy, i) = tmp24; 356 } 357 } 358 359 360 /*-----------------------------------------------------------------------------361 * 用于保存FFT图像362 *363 *364 *-----------------------------------------------------------------------------*/365 Mat img(¢er_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: FFT388 * Description: 计算一维FFT389 * =====================================================================================390 */391 void Cfft_dlgDlg::FFT(complex * TD, complex * 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 *W,*X1,*X2,*X;402 // 计算付立叶变换点数403 count = 1 << r;404 // 分配运算所需存储器405 W = new complex [count / 2];406 X1 = new complex [count];407 X2 = new complex [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 (cos(angle), sin(angle));414 }415 // 将时域点写入X1416 memcpy(X1, TD, sizeof(complex ) * 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<
运行效果如下：
提供程序一份：
http://pan.baidu.com/s/1ehvwy