关于FFT，书上已经给出了实现方法；曾在研2时也使用迭代法实现了自己的FFT，速度上要慢一些，但是理解起来要容易一些；
最近看书，发现了一些以前没有注意到的问题；比如，FFT产生是到底是什么呢？是频率的信息吗？完整吗？程序表现出来的结果到底正确吗？等等一些问题；以前没有考虑过。
今天来给出答案，当然是本人的一些个人理解，不一定正确!
一，FFT产生的到底是什么？
书上曾经把FFT后的信息用一幅图像来表示，其实这很勉强；目的只是让大家了解频域到底反应了一些什么东西，比如，边缘上在频域上就可以很直观的反应出来；别说一幅图像了，就是多幅图像也不一定能表示得了，因为频域变换出来的数据不是整型，如果强制为整型，则会造成数据的丢失；除此之外，更重要的一点，一幅图像的频域只反应了一个频谱，也就是相位并没有考虑进去，这样，通过FFT产生的频域来还原时域图像是不可能的。
二，这些频率信息完整吗？
如上所述，这些信息是不完整的，我曾试图使用一幅24位的bmp来表示频率的信息，但是最终因为图像的整数表示而放弃，当然也是可以做的，只是效果要差很多；以后有机会我会试一试的。频域应该是一个复数；在我看来，频域的的完整信息应该有两部分表示，一部分是实部，一部分是虚部；
三，以前图像处理程序中给出的实现正确吗？
首先，我应该说是正确的，因为频率上确实是这样的；但是，以前程序在为了更加细节描述频域的信息，并没有除M*N，取代除了个100；这样更能突出细节；
关于反向FFT的算法其实很简单，使用前向变换来计算就很方便!我试着用频域信息还原时域图像，有很少的地方失真，失真的原因是由于计算机在计算时浮点数的省略；
我还原时域图像的方法是这样的：使用FFT对时域图像进行变换，得到频域的信息，保存在两个txt文件当中，这两个文件分别是频域的实部和虚部!然后使用前向变换，用txt文件来构建时域图像。
代码如下：
.h文件
//FFT.h:interfacefortheCFFTclass.
//
//
#if!defined(AFX_FFT_H__C3BE732F_8D82_4870_B422_25580A5BABC9__INCLUDED_)
#defineAFX_FFT_H__C3BE732F_8D82_4870_B422_25580A5BABC9__INCLUDED_
#if_MSC_VER>1000
#pragmaonce
#endif//_MSC_VER>1000
/**********************************************
*类名:FFT及反向FFT
*目的:时域到频域，频域到时域的变换
*作者:无忧狂澜
*日期:-02-19
**********************************************/
classCFFT
{
public:
BOOLComputeConverseFFT(LONGlWidth,
LONGlHeight,
double*fFReal,
double*fFImag,
double*fTReal,
double*fTImag,
boolbCenter);
BOOLConverseFourier(LPSTRlpDIBBits,
LONGlWidth,
LONGlHeight,
double*fFReal,
double*fFImag);
BOOLFourier(LPSTRlpDIBBits,LONGlWidth,LONGlHeight,double*fReal,double*fImag,boolbCenter);
CFFT();
virtual~CFFT();
};
#endif//!defined(AFX_FFT_H__C3BE732F_8D82_4870_B422_25580A5BABC9__INCLUDED_)
.cpp文件
[cpp] view plaincopyprint?
//FFT.cpp:implementationoftheCFFTclass.
//
//
#include"stdafx.h"
#include"mydip.h"
#include"FFT.h"
#include"math.h"
#include<complex>
usingnamespacestd;
#definePI3.1415926
#defineTWOPI6.2831852
#ifdef_DEBUG
#undefTHIS_FILE
staticcharTHIS_FILE[]=__FILE__;
#definenewDEBUG_NEW
#endif
//
//Construction/Destruction
//
/*************************************************************************
*
*函数名称：
*FFT()
*
*参数:
*complex<double>*TD-指向时域数组的指针
*complex<double>*FD-指向频域数组的指针
*r-2的幂数，即迭代次数
*
*返回值:
*无。
*
*说明:
*用迭代法实现的快速付立叶变换。
*速度比蝶形算法要慢一些，但理解起来容易
************************************************************************/
/*
voidinit(intr)
{
f=newcomplex<double>[1<<r];
inttempIndex,i,j;
doubleangle;
W=newcomplex<double>*[1<<r];
for(tempIndex=1;tempIndex<=r;tempIndex++)
{
W[1<<tempIndex]=newcomplex<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>*WINAPIFFTMySelf(intn,intxr,intxi)
{
intu;//临时变量，用于循环计数
complex<double>*FMeven,*FModd;//用于指向上次FFT的运算结果
complex<double>*F=newcomplex<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);
}
elseif(n==1)
F[0]=f[xi];//如果是一点变换，FFT就是本身
returnF;
}
VOIDWINAPIFFT(complex<double>*TD,complex<double>*FD,intr)
{
inttempCount=1<<r;
init(r);
//将时域点写入f
complex<double>*tempTD=newcomplex<double>[tempCount];
memcpy(f,TD,sizeof(complex<double>)*(tempCount));
tempTD=FFTMySelf(tempCount,1,0);
inti;
for(i=0;i<tempCount;i++)
{
FD[i]=tempTD[i]*complex<double>(tempCount,0);
}
}
*/
/*************************************************************************
*
*函数名称：
*FFT()
*
*参数:
*complex<double>*TD-指向时域数组的指针
*complex<double>*FD-指向频域数组的指针
*r-2的幂数，即迭代次数
*
*返回值:
*无。
*
*说明:
*该函数用来实现快速付立叶变换。
*
************************************************************************/
VOIDWINAPIFFT(complex<double>*TD,complex<double>*FD,intr)
{
//付立叶变换点数
LONGcount;
//循环变量
inti,j,k;
//中间变量
intbfsize,p;
//角度
doubleangle;
complex<double>*W,*X1,*X2,*X;
//计算付立叶变换点数
count=1<<r;
//分配运算所需存储器
W=newcomplex<double>[count/2];
X1=newcomplex<double>[count];
X2=newcomplex<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];
}
//释放内存
deleteW;
deleteX1;
deleteX2;
}
CFFT::CFFT()
{
}
CFFT::~CFFT()
{
}
/*************************************************************************
*功能:正向傅里叶变换
*最后两个参数表示变换后频域的实部和虚部
*注意:频域的实部和虚部我没有进行中心变换
*因此:在逆向傅里叶变换时不需要进行中心逆变换
*bCenter表示是否中心化
*************************************************************************/
BOOLCFFT::Fourier(LPSTRlpDIBBits,LONGlWidth,LONGlHeight,double*fReal,double*fImag,boolbCenter)
{
//指向源图像的指针
unsignedchar*lpSrc;
//中间变量
doubledTemp;
//循环变量
LONGi;
LONGj;
//进行付立叶变换的宽度和高度（2的整数次方）
LONGw;
LONGh;
intwp;
inthp;
//图像每行的字节数
LONGlLineBytes;
//计算图像每行的字节数
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=newcomplex<double>[w*h];
complex<double>*FD=newcomplex<double>[w*h];
//行
for(i=0;i<h;i++)
{
//列
for(j=0;j<w;j++)
{
//指向DIB第i行，第j个象素的指针
lpSrc=(unsignedchar*)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=(unsignedchar*)lpDIBBits+lLineBytes*(lHeight-1-i)+j;
//lpSrc=(unsignedchar*)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);
}
}
//删除临时变量
deleteTD;
deleteFD;
//返回
returnTRUE;
}
/*************************************************************************
*功能:逆向向傅里叶变换(使用前向变换算法)
*最后两个参数表示变换后频域的实部和虚部
*注意:频域的实部和虚部我没有进行中心变换
*因此:在逆向傅里叶变换时不需要进行中心逆变换
*fFReal,fFImag，表示频域
*fTReal,fTImag，表示时域
*************************************************************************/
BOOLCFFT::ConverseFourier(LPSTRlpDIBBits,
LONGlWidth,
LONGlHeight,
double*fFReal,
double*fFImag)
{
//指向源图像的指针
unsignedchar*lpSrc;
//中间变量
doubledTemp;
//循环变量
LONGi;
LONGj;
//进行付立叶变换的宽度和高度（2的整数次方）
LONGw;
LONGh;
intwp;
inthp;
//图像每行的字节数
LONGlLineBytes;
//计算图像每行的字节数
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=newcomplex<double>[w*h];
complex<double>*FD=newcomplex<double>[w*h];
//行
for(i=0;i<h;i++)
{
//列
for(j=0;j<w;j++)
{
//指向DIB第i行，第j个象素的指针
//lpSrc=(unsignedchar*)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=(unsignedchar*)lpDIBBits+lLineBytes*(lHeight-1-i)+j;
lpSrc=(unsignedchar*)lpDIBBits+lLineBytes*
(lHeight-1-i)+j;
//更新源图像
*(lpSrc)=(BYTE)(dTemp);
}
}
//删除临时变量
deleteTD;
deleteFD;
//返回
returnTRUE;
}
//根据频域计算时域
//bCenter表示是否中心化
BOOLCFFT::ComputeConverseFFT(LONGlWidth,
LONGlHeight,
double*fFReal,
double*fFImag,
double*fTReal,
double*fTImag,
boolbCenter)
{
//循环变量
LONGi;
LONGj;
//进行付立叶变换的宽度和高度（2的整数次方）
LONGw;
LONGh;
intwp;
inthp;
doublefTemp,fTemp1,fTemp2;
//图像每行的字节数
LONGlLineBytes;
//计算图像每行的字节数
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=newcomplex<double>[w*h];
complex<double>*FD=newcomplex<double>[w*h];
//行
for(i=0;i<h;i++)
{
//列
for(j=0;j<w;j++)
{
//指向DIB第i行，第j个象素的指针
//lpSrc=(unsignedchar*)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;
}
}
//删除临时变量
deleteTD;
deleteFD;
//返回
returnTRUE;
}
为了方便大家了解，现在把附加代码也粘上来
[cpp] view plaincopyprint?
//正向FFT
voidCMyDIPView::OnFft()
{
//TODO:Addyourcommandhandlercodehere
//获取文档
CMyDIPDoc*pDoc=GetDocument();
//指向DIB的指针
LPSTRlpDIB;
//指向DIB象素指针
LPSTRlpDIBBits;
double*fReal,*fImag;
intnWidth,nHeight,i,j;
CFFTfft;
//锁定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=newdouble[nWidth*nHeight];
fImag=newdouble[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());
//把数据写入文件中
CStringstrTemp;
//写频域实部数据
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();
}
//逆向傅里叶变换
//通过频域信息还原时域图像
voidCMyDIPView::OnConversefft()
{
//TODO:Addyourcommandhandlercodehere
//TODO:Addyourcommandhandlercodehere
//获取文档
CMyDIPDoc*pDoc=GetDocument();
//指向DIB的指针
LPSTRlpDIB;
//指向DIB象素指针
LPSTRlpDIBBits;
double*fReal,*fImag;
intnWidth,nHeight,i,j;
CFFTfft;
//锁定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=newdouble[nWidth*nHeight];
fImag=newdouble[nWidth*nHeight];
//从文件读取频域数据
CStringstrTemp;
//写频域实部数据
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();
}
![[推荐]中国达人秀：刘伟的1500字作文](https://1500zi.500zi.com/uploadfile/img/2024/08/28/d8eef0d8c8c197a59e7d5cc3ca1e9f0f.jpg)
