图像的傅立叶变换

实验三、图像的傅立叶变换一、实验目的1.了解图像变换的意义和手段；2.熟悉傅里叶变换的基本性质；3.熟练掌握FFT的方法及应用；4.通过实验了解二维频谱的分布特点；5.通过本实验掌握编程实现数字图像的傅立叶变换。

二、实验原理1.应用傅立叶变换进行图像处理傅里叶变换是线性系统分析的一个有力工具，它能够定量地分析诸如数字化系统、采样点、电子放大器、卷积滤波器、噪音和显示点等的作用。

通过实验培养这项技能，将有助于解决大多数图像处理问题。

对任何想在工作中有效应用数字图像处理技术的人来说，把时间用在学习和掌握博里叶变换上是很有必要的。

2.傅立叶（Fourier）变换的定义对于二维信号，二维Fourier变换定义为：二维离散傅立叶变换为：图像的傅立叶变换与一维信号的傅立叶变换变换一样，有快速算法，具体参见参考书目，有关傅立叶变换的快速算法的程序不难找到。

实际上，现在有实现傅立叶变换的芯片，可以实时实现傅立叶变换。

三、实验要求：1.实验之前要预习2.独立完成程序的编写3.写出实验报告4.实验每组1人5.实验设备：每组计算机一台四、实验内容1.实验程序和流程图1）程序整体流程图2）傅里叶算法分析通常的FFT算法：直接计算离散傅立叶变换具有n^2的复杂度，而cooley 和tukey在1965年发现了一种计算离散傅立叶变换的快速算法(即通常所说的FFT算法)，这个算法在计算变换长度n=2^k的离散傅立叶变换时，具有n*k的复杂度，即O(n)=n*log2(n),下面以此为例，讲讲快FFT的特点。

（1）复数运算：傅立叶变换是基于复数的，因此首先知道复数的运算规则，在FFT算法中，只涉及复数的加、减和乘法三种运算。

一个复数可表示为c=(x,yi),x 为复数的实部，y为复数的虚部，i为虚数单位，等于-1的平方根。

复数的运算规则是:若c1表示为(x1,y1),c2表示为(x2,y2),则(x1+x2,y1+y2)和(x1-x2,y1-y2)分别等于c1+c2的和，c1-c2的差，复数的乘法相对复杂一些，c1*c2的积为(x1*x2-y1*y2,x1*y2+x2*y1).（2）蝶形变换：普通的FFT算法称为基2的FFT算法，这种算法的核心是蝶形变换长度为n=2^k1的变换共需要做k1*n/2次蝶形变换，若需变换数据表示为一个复数数组c[],则每次蝶形变换有2个输入c[i],c[i+s],两个输出:c[i],c[i+s]，s成为翅间距。

每个变换的基本算法是：t=wr*c[i+s];c[i+s]=c[i]-t;c[i]=c[i]+t;前面说过，长度为n=2^k1的变换共需要做k1*n/2次变换，实际的程序是一个3层循环,共需要k1*k2*(k3/2)次变换(k2*k3/2=n/2)。

前面的wr是w的整数次方，w=e^(-2*PI/k3)(k3=2,4,8,16...n,PI是圆周率)，也成为旋转因子，例如n=32的变换需要log2(32)=5趟变换:第1趟变换需要16*1次变换，翅间距是1,若w=e^(-2*PI/2)，则wr=w^1 第2趟变换需要8*2次变换，翅间距是2,若w=e^(-2*PI/4)，则wr=w^1,w^2 第3趟变换需要4*2次变换，翅间距是4,若w=e^(-2*PI/8)，则wr=w^1,w^2,w^3,w^4第4趟变换需要2*8次变换，翅间距是8,若w=e^(-2*PI/16)，则wr=w^1,w^2,w^3,w^4,w^5,w^6,w^7,w^8第5趟变换需要16*1次变换，翅间距是16,若w=e^(-2*PI/32)，则wr=w^1,w^2,w^3,w^4,w^5...w^15,w^16（3）w数组，w的实部=cos(2*PI/k3),w的虚部=-sin(2*PI/k3),计算出w,则wr 数组就好求了，不断即相乘即可，当然也可以通过三角函数直接求。

w^p的实部=cos(2*PI/K3*p),虚部=-sin(2*PI/k3*p)（4）复数数组排序，在基2的蝶形变换中，复数数组需要重新排序，c[i]要放置到数组c的第reverse(c[i])的位置，m=reverse(n)函数的算法是这样的，若n的k位2进制的为b[],b[k-1],B[k-2],...b[2],b[1],b[0],(b[i]等于1或者0,b[0]为最低bit).则m=reverse(n)的2进制的为b[0],b[1],b[2],b[3],...b[k-1](b[k-1]为最低bit).程序代码// MyDIPDoc.cpp : implementation of the CMyDIPDoc class//#include "stdafx.h"#include "MyDIP.h"#include "MyDIPDoc.h"#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif///////////////////////////////////////////////////////////////////////////// // CMyDIPDocIMPLEMENT_DYNCREATE(CMyDIPDoc, CDocument)BEGIN_MESSAGE_MAP(CMyDIPDoc, CDocument)//{{AFX_MSG_MAP(CMyDIPDoc)// NOTE - the ClassWizard will add and remove mapping macros here.// DO NOT EDIT what you see in these blocks of generated code!//}}AFX_MSG_MAPEND_MESSAGE_MAP()///////////////////////////////////////////////////////////////////////////// // CMyDIPDoc construction/destructionCMyDIPDoc::CMyDIPDoc(){// TODO: add one-time construction code herem_sizeDoc.cx = 800;m_sizeDoc.cy = 600;m_palDIB = NULL;m_hDIB = NULL;}CMyDIPDoc::~CMyDIPDoc(){if (m_hDIB != NULL){::GlobalFree((HGLOBAL) m_hDIB);m_hDIB=NULL;}if (m_palDIB != NULL){delete m_palDIB;m_palDIB=NULL;}}BOOL CMyDIPDoc::OnNewDocument(){if (!CDocument::OnNewDocument())return FALSE;// TODO: add reinitialization code here// (SDI documents will reuse this document)m_hDIB = NewDIB(m_sizeDoc.cx, m_sizeDoc.cy,8);InitDIBData();return TRUE;}///////////////////////////////////////////////////////////////////////////// // CMyDIPDoc serializationvoid CMyDIPDoc::Serialize(CArchive& ar){if (ar.IsStoring()){// TODO: add storing code here}else{// TODO: add loading code here}}///////////////////////////////////////////////////////////////////////////// // CMyDIPDoc diagnostics#ifdef _DEBUGvoid CMyDIPDoc::AssertValid() const{CDocument::AssertValid();}void CMyDIPDoc::Dump(CDumpContext& dc) const{CDocument::Dump(dc);}#endif //_DEBUG/////////////////////////////////////////////////////////////////////////////// CMyDIPDoc commandsBOOL CMyDIPDoc::ReadImgFile(CString sName){if (m_hDIB != NULL){::GlobalFree((HGLOBAL) m_hDIB);m_hDIB=NULL;}if (m_palDIB != NULL){delete m_palDIB;m_palDIB=NULL;}// replace calls to Serialize with ReadDIBFile function CFile nFile;if(!nFile.Open(sName,CFile::modeRead))return false;m_hDIB = ReadDIBFile(nFile);nFile.Close();InitDIBData();if (m_hDIB == NULL)return FALSE;return TRUE;}BOOL CMyDIPDoc::InitDIBData(){if (m_palDIB != NULL){delete m_palDIB;m_palDIB = NULL;}if (m_hDIB == NULL){return FALSE;}m_palDIB = new CPalette;if (m_palDIB == NULL){// we must be really low on memory::GlobalFree((HGLOBAL) m_hDIB);m_hDIB = NULL;return FALSE;}if (::CreateDIBPalette(m_hDIB, m_palDIB) == NULL){// DIB may not have a palettedelete m_palDIB;m_palDIB = NULL;}return TRUE;}BOOL CMyDIPDoc::OnOpenDocument(LPCTSTR lpszPathName){if (!CDocument::OnOpenDocument(lpszPathName))return FALSE;// TODO: Add your specialized creation code hereBOOL b = ReadImgFile(lpszPathName);if(b){LPSTR lpDIB = (LPSTR) ::GlobalLock((HGLOBAL) m_hDIB);m_sizeDoc.cx = (int) ::DIBWidth(lpDIB); // Size of DIB - xm_sizeDoc.cy = (int) ::DIBHeight(lpDIB); // Size of DIB - y::GlobalUnlock((HGLOBAL) m_hDIB);}POSITION pos;pos = GetFirstViewPosition();while(pos != NULL){CScrollView * pView = (CScrollView *)GetNextView(pos);if(pView != NULL){pView->SetScrollSizes(MM_TEXT, m_sizeDoc);}}return b;}BOOL CMyDIPDoc::OnSaveDocument(LPCTSTR lpszPathName){// TODO: Add your specialized code here and/or call the base classCFile nFile;if(m_hDIB == NULL)return false;nFile.Open(lpszPathName,CFile::modeWrite | CFile::modeCreate);SaveDIB(m_hDIB, nFile);nFile.Close();return true;//return CDocument::OnSaveDocument(lpszPathName);}3）furious算法程序实现傅里叶算法流程图furious算法流程图程序代码// MyDIPView.cpp : implementation of the CMyDIPView class//#include "stdafx.h"#include "MyDIP.h"#include "MyDIPDoc.h"#include "MyDIPView.h"#include "math.h"#include<complex>using namespace std;#define PI 3.14159#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif/************************************************************************** 函数名称：* FFT()* 参数:* complex<double> * TD - 指向时域数组的指针* complex<double> * FD - 指向频域数组的指针* r －2的幂数，即迭代次数* 返回值:* 无。