关于FFT，书上已经给出了实现方法；曾在研2时也使用迭代法实现了自己的FFT，速度上要慢一些，但是理解起来要容易一些；
最近看书，发现了一些以前没有注意到的问题；比如，FFT产生是到底是什么呢？是频率的信息吗？完整吗？程序表现出来的结果到底正确吗？等等一些问题；以前没有考虑过。
今天来给出答案，当然是本人的一些个人理解，不一定正确!
一，FFT产生的到底是什么？
书上曾经把FFT后的信息用一幅图像来表示，其实这很勉强；目的只是让大家了解频域到底反应了一些什么东西，比如，边缘上在频域上就可以很直观的反应出来；别说一幅图像了，就是多幅图像也不一定能表示得了，因为频域变换出来的数据不是整型，如果强制为整型，则会造成数据的丢失；除此之外，更重要的一点，一幅图像的频域只反应了一个频谱，也就是相位并没有考虑进去，这样，通过FFT产生的频域来还原时域图像是不可能的。
二，这些频率信息完整吗？
如上所述，这些信息是不完整的，我曾试图使用一幅24位的bmp来表示频率的信息，但是最终因为图像的整数表示而放弃，当然也是可以做的，只是效果要差很多；以后有机会我会试一试的。频域应该是一个复数；在我看来，频域的的完整信息应该有两部分表示，一部分是实部，一部分是虚部；
三，以前图像处理程序中给出的实现正确吗？
首先，我应该说是正确的，因为频率上确实是这样的；但是，以前程序在为了更加细节描述频域的信息，并没有除M*N，取代除了个100；这样更能突出细节；
关于反向FFT的算法其实很简单，使用前向变换来计算就很方便!我试着用频域信息还原时域图像，有很少的地方失真，失真的原因是由于计算机在计算时浮点数的省略；
我还原时域图像的方法是这样的：使用FFT对时域图像进行变换，得到频域的信息，保存在两个txt文件当中，这两个文件分别是频域的实部和虚部!然后使用前向变换，用txt文件来构建时域图像。
代码如下：
.h文件
// FFT.h: interface for the CFFT class.
//
//
#if !defined(AFX_FFT_H__C3BE732F_8D82_4870_B422_25580A5BABC9__INCLUDED_)
#define AFX_FFT_H__C3BE732F_8D82_4870_B422_25580A5BABC9__INCLUDED_
#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
/**********************************************
*类名:FFT及反向FFT
*目的:时域到频域，频域到时域的变换
*作者:无忧狂澜
*日期:2010-02-19
**********************************************/
class CFFT
{
public:
BOOL ComputeConverseFFT(LONG lWidth,
LONG lHeight,
double *fFReal,
double *fFImag,
double *fTReal,
double *fTImag,
bool bCenter);
BOOL ConverseFourier(LPSTR lpDIBBits,
LONG lWidth,
LONG lHeight,
double *fFReal,
double *fFImag);
BOOL Fourier(LPSTR lpDIBBits, LONG lWidth, LONG lHeight,double *fReal,double *fImag,bool bCenter);
CFFT();
virtual ~CFFT();
};
#endif // !defined(AFX_FFT_H__C3BE732F_8D82_4870_B422_25580A5BABC9__INCLUDED_)
.cpp文件
[cpp] view plaincopyprint?
// FFT.cpp: implementation of the CFFT class.
//
//
#include "stdafx.h"
#include "mydip.h"
#include "FFT.h"
#include "math.h"
#include <complex>
using namespace std;
#define PI 3.1415926
#define TWOPI 6.2831852
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
//
// Construction/Destruction
//
/*************************************************************************
*
* 函数名称：
* FFT()
*
* 参数:
* complex<double> * TD - 指向时域数组的指针
* complex<double> * FD - 指向频域数组的指针
* r - 2的幂数，即迭代次数
*
* 返回值:
* 无。
*
* 说明:
* 用迭代法实现的快速付立叶变换。
* 速度比蝶形算法要慢一些，但理解起来容易
************************************************************************/
/*
void init(int r)
{
f=new complex<double>[1<<r];
int tempIndex,i,j;
double angle;
W=new complex<double>*[1<<r];
for(tempIndex=1;tempIndex<=r;tempIndex++)
{
W[1<<tempIndex]=new complex<double>[(1<<(r-1))-1];
}
for(i=2;i<=(1<<r);i=i<<1)
for(j=0;j<=(i>>1)-1;j++)
{
angle = -j*TWOPI/i;
W[i][j]=complex<double>(cos(angle),sin(angle));
}
}
complex<double>* WINAPI FFTMySelf(int n,int xr,int xi)
{
int u; //临时变量，用于循环计数
complex<double> *FMeven,*FModd; //用于指向上次FFT的运算结果
complex<double> *F=new complex<double>[n]; //用于保存此次FFT运算结果
if(n/2>=1)
{
FMeven=FFTMySelf(n/2,2*xr,xi); //迭代FFT
FModd=FFTMySelf(n/2,2*xr,xr+xi);
for(u=0;u<n/2;u++)
{
F[u]=(1/(double)2)*(FMeven[u]+FModd[u]*W[n][u]); //根据公式计算FFT
F[u+n/2]=(1/(double)2)*(FMeven[u]-FModd[u]*W[n][u]);
}
free(FMeven); //释放所占用的堆内存
free(FModd);
}
else if(n==1)
F[0]=f[xi]; //如果是一点变换，FFT就是本身
return F;
}
VOID WINAPI FFT(complex<double> * TD, complex<double> * FD, int r)
{
int tempCount=1<<r;
init(r);
// 将时域点写入f
complex<double> *tempTD=new complex<double>[tempCount];
memcpy(f, TD, sizeof(complex<double>)*(tempCount));
tempTD=FFTMySelf(tempCount,1,0);
int i;
for(i=0;i<tempCount;i++)
{
FD[i]=tempTD[i]*complex<double>(tempCount,0);
}
}
*/
/*************************************************************************
*
* 函数名称：
* FFT()
*
* 参数:
* complex<double> * TD - 指向时域数组的指针
* complex<double> * FD - 指向频域数组的指针
* r - 2的幂数，即迭代次数
*
* 返回值:
* 无。
*
* 说明:
* 该函数用来实现快速付立叶变换。
*
************************************************************************/
VOID WINAPI 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 = -i * PI * 2 / 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;
}
CFFT::CFFT()
{
}
CFFT::~CFFT()
{
}
/*************************************************************************
*功能:正向傅里叶变换
*最后两个参数表示变换后频域的实部和虚部
*注意:频域的实部和虚部我没有进行中心变换
*因此:在逆向傅里叶变换时不需要进行中心逆变换
*bCenter表示是否中心化
*************************************************************************/
BOOL CFFT::Fourier(LPSTR lpDIBBits, LONG lWidth, LONG lHeight, double *fReal, double *fImag,bool bCenter)
{
// 指向源图像的指针
unsigned char* lpSrc;
// 中间变量
double dTemp;
// 循环变量
LONG i;
LONG j;
// 进行付立叶变换的宽度和高度（2的整数次方）
LONG w;
LONG h;
int wp;
int hp;
// 图像每行的字节数
LONG lLineBytes;
// 计算图像每行的字节数
lLineBytes = WIDTHBYTES(lWidth * 8);
// 赋初值
w = 1;
h = 1;
wp = 0;
hp = 0;
// 计算进行付立叶变换的宽度和高度（2的整数次方）
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}
while(h * 2 <= lHeight)
{
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++)
{
// 指向DIB第i行，第j个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
// 给时域赋值
if(bCenter)
{//如果中心化
if((i+j)%2==0)
dTemp=*(lpSrc);
else
dTemp=*(lpSrc)*-1;
}
else
{
dTemp=*(lpSrc);
}
TD[j + w * i] = complex<double>(dTemp, 0);
}
}
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;
//注意:此处没有进行中心变换
fReal[j * h + i]=FD[j * h + i].real()/(lWidth*lHeight);
fImag[j * h + i]=FD[j * h + i].imag()/(lWidth*lHeight);
// 判断是否超过255
if (dTemp > 255)
{
// 对于超过的，直接设置为255
dTemp = 255;
}
// 指向DIB第(i<h/2 ? i+h/2 : i-h/2)行，第(j<w/2 ? j+w/2 : j-w/2)个象素的指针
// 此处不直接取i和j，是为了将变换后的原点移到中心
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
//lpSrc = (unsigned char*)lpDIBBits + lLineBytes *
//(lHeight - 1 - (i<h/2 ? i+h/2 : i-h/2)) + (j<w/2 ? j+w/2 : j-w/2);
// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
}
// 删除临时变量
delete TD;
delete FD;
// 返回
return TRUE;
}
/*************************************************************************
*功能:逆向向傅里叶变换(使用前向变换算法)
*最后两个参数表示变换后频域的实部和虚部
*注意:频域的实部和虚部我没有进行中心变换
*因此:在逆向傅里叶变换时不需要进行中心逆变换
*fFReal,fFImag，表示频域
*fTReal,fTImag，表示时域
*************************************************************************/
BOOL CFFT::ConverseFourier(LPSTR lpDIBBits,
LONG lWidth,
LONG lHeight,
double *fFReal,
double *fFImag)
{
// 指向源图像的指针
unsigned char* lpSrc;
// 中间变量
double dTemp;
// 循环变量
LONG i;
LONG j;
// 进行付立叶变换的宽度和高度（2的整数次方）
LONG w;
LONG h;
int wp;
int hp;
// 图像每行的字节数
LONG lLineBytes;
// 计算图像每行的字节数
lLineBytes = WIDTHBYTES(lWidth * 8);
// 赋初值
w = 1;
h = 1;
wp = 0;
hp = 0;
// 计算进行付立叶变换的宽度和高度（2的整数次方）
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}
while(h * 2 <= lHeight)
{
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++)
{
// 指向DIB第i行，第j个象素的指针
//lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
// 给频域赋值
TD[j + w * i] = complex<double>(fFReal[j*h + i], fFImag[j*h + i]*-1);
}
}
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());
// 判断是否超过255
if (dTemp > 255)
{
// 对于超过的，直接设置为255
dTemp = 255;
}
// 指向DIB第(i<h/2 ? i+h/2 : i-h/2)行，第(j<w/2 ? j+w/2 : j-w/2)个象素的指针
// 此处不直接取i和j，是为了将变换后的原点移到中心
//lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
lpSrc = (unsigned char*)lpDIBBits + lLineBytes *
(lHeight - 1 - i) + j;
// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
}
// 删除临时变量
delete TD;
delete FD;
// 返回
return TRUE;
}
//根据频域计算时域
//bCenter表示是否中心化
BOOL CFFT::ComputeConverseFFT(LONG lWidth,
LONG lHeight,
double *fFReal,
double *fFImag,
double *fTReal,
double *fTImag,
bool bCenter)
{
// 循环变量
LONG i;
LONG j;
// 进行付立叶变换的宽度和高度（2的整数次方）
LONG w;
LONG h;
int wp;
int hp;
double fTemp,fTemp1,fTemp2;
// 图像每行的字节数
LONG lLineBytes;
// 计算图像每行的字节数
lLineBytes = WIDTHBYTES(lWidth * 8);
// 赋初值
w = 1;
h = 1;
wp = 0;
hp = 0;
// 计算进行付立叶变换的宽度和高度（2的整数次方）
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}
while(h * 2 <= lHeight)
{
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++)
{
// 指向DIB第i行，第j个象素的指针
//lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
// 给频域赋值
if(fFImag==NULL)
fTemp=0;
else
fTemp=fFImag[j*h + i]*-1;
if(bCenter)
{//如果中心化
if((i+j)%2==0)
{
fTemp1=fFReal[j*h + i];
fTemp2=fTemp;
}
else
{
fTemp1=fFReal[j*h + i]*-1;
fTemp2=fTemp*-1;
}
}
else
{
fTemp1=fFReal[j*h + i];
fTemp2=fTemp;
}
TD[j + w * i] = complex<double>(fTemp1, fTemp2);
}
}
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++)
{
// 计算频谱
//FFT-1是不用除采样周期的哦
fTReal[j * h + i]=FD[j * h + i].real();
fTImag[j * h + i]=FD[j * h + i].imag()*-1;
}
}
// 删除临时变量
delete TD;
delete FD;
// 返回
return TRUE;
}
为了方便大家了解，现在把附加代码也粘上来
[cpp] view plaincopyprint?
//正向FFT
void CMyDIPView::OnFft()
{
// TODO: Add your command handler code here
// 获取文档
CMyDIPDoc* pDoc = GetDocument();
// 指向DIB的指针
LPSTR lpDIB;
// 指向DIB象素指针
LPSTR lpDIBBits;
double *fReal,*fImag;
int nWidth,nHeight,i,j;
CFFT fft;
// 锁定DIB
lpDIB = (LPSTR) ::GlobalLock((HGLOBAL) pDoc->GetHDIB());
// 找到DIB图像象素起始位置
lpDIBBits = ::FindDIBBits(lpDIB);
// 判断是否是8-bpp位图（这里为了方便，只处理8-bpp位图的付立叶变换，其它的可以类推）
if (::DIBNumColors(lpDIB) != 256)
{
// 提示用户
MessageBox("目前只支持256色位图的付立叶变换！", "系统提示" ,
MB_ICONINFORMATION | MB_OK);
// 解除锁定
::GlobalUnlock((HGLOBAL) pDoc->GetHDIB());
// 返回
return;
}
// 更改光标形状
BeginWaitCursor();
nWidth=::DIBWidth(lpDIB);
nHeight=::DIBHeight(lpDIB);
fReal=new double[nWidth*nHeight];
fImag=new double[nWidth*nHeight];
// 调用Fourier()函数进行付立叶变换
if (fft.Fourier(lpDIBBits, nWidth, nHeight,fReal,fImag))
{
// 设置脏标记
pDoc->SetModifiedFlag(TRUE);
// 更新视图
pDoc->UpdateAllViews(NULL);
}
else
{
// 提示用户
MessageBox("分配内存失败！", "系统提示" , MB_ICONINFORMATION | MB_OK);
}
// 解除锁定
::GlobalUnlock((HGLOBAL) pDoc->GetHDIB());
//把数据写入文件中
CString strTemp;
//写频域实部数据
m_file=fopen("频域实部.txt","w+");
//行
for(i = 0; i < nHeight; i++)
{
// 列
for(j = 0; j < nWidth; j++)
{
strTemp.Format("%f",fReal[j * nHeight + i]);
fwrite(strTemp,strTemp.GetLength(),1,m_file);
strTemp=" ";
fwrite(strTemp,strTemp.GetLength(),1,m_file);
}
}
fclose(m_file);
//写频域虚部数据
m_file=fopen("频域虚部.txt","w+");
//行
for(i = 0; i < nHeight; i++)
{
// 列
for(j = 0; j < nWidth; j++)
{
strTemp.Format("%f",fImag[j * nHeight + i]);
fwrite(strTemp,strTemp.GetLength(),1,m_file);
strTemp=" ";
fwrite(strTemp,strTemp.GetLength(),1,m_file);
}
}
fclose(m_file);
if(fReal)
{
delete []fReal;
fReal=NULL;
}
if(fImag)
{
delete []fImag;
fImag=NULL;
}
// 恢复光标
EndWaitCursor();
}
//逆向傅里叶变换
//通过频域信息还原时域图像
void CMyDIPView::OnConversefft()
{
// TODO: Add your command handler code here
// TODO: Add your command handler code here
// 获取文档
CMyDIPDoc* pDoc = GetDocument();
// 指向DIB的指针
LPSTR lpDIB;
// 指向DIB象素指针
LPSTR lpDIBBits;
double *fReal,*fImag;
int nWidth,nHeight,i,j;
CFFT fft;
// 锁定DIB
lpDIB = (LPSTR) ::GlobalLock((HGLOBAL) pDoc->GetHDIB());
// 找到DIB图像象素起始位置
lpDIBBits = ::FindDIBBits(lpDIB);
// 判断是否是8-bpp位图（这里为了方便，只处理8-bpp位图的付立叶变换，其它的可以类推）
if (::DIBNumColors(lpDIB) != 256)
{
// 提示用户
MessageBox("目前只支持256色位图的付立叶变换！", "系统提示" ,
MB_ICONINFORMATION | MB_OK);
// 解除锁定
::GlobalUnlock((HGLOBAL) pDoc->GetHDIB());
// 返回
return;
}
// 更改光标形状
BeginWaitCursor();
nWidth=::DIBWidth(lpDIB);
nHeight=::DIBHeight(lpDIB);
fReal=new double[nWidth*nHeight];
fImag=new double[nWidth*nHeight];
//从文件读取频域数据
CString strTemp;
//写频域实部数据
m_file=fopen("频域实部.txt","r");
//行
for(i = 0; i < nHeight; i++)
{
// 列
for(j = 0; j < nWidth; j++)
{
fscanf(m_file,"%lf ",&fReal[j * nHeight + i]);
}
}
fclose(m_file);
//写频域虚部数据
m_file=fopen("频域虚部.txt","r");
//行
for(i = 0; i < nHeight; i++)
{
// 列
for(j = 0; j < nWidth; j++)
{
fscanf(m_file,"%lf ",&fImag[j * nHeight + i]);
}
}
// 调用Fourier()函数进行付立叶变换
if (fft.ConverseFourier(lpDIBBits, nWidth, nHeight,fReal,fImag))
{
// 设置脏标记
pDoc->SetModifiedFlag(TRUE);
// 更新视图
pDoc->UpdateAllViews(NULL);
}
else
{
// 提示用户
MessageBox("分配内存失败！", "系统提示" , MB_ICONINFORMATION | MB_OK);
}
// 解除锁定
::GlobalUnlock((HGLOBAL) pDoc->GetHDIB());
if(fReal)
{
delete []fReal;
fReal=NULL;
}
if(fImag)
{
delete []fImag;
fImag=NULL;
}
// 恢复光标
EndWaitCursor();
}
349
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
