用AS3实现一维傅里叶变换，一维快速傅里叶变换，二维快速傅里叶变换。以及反变换。
参考：[精通Visual.Cpp数字图像处理典型算法及实现(第2版)].张宏林.扫描版
下载地址链接: http://pan.baidu.com/s/1bnCtVjD 密码: kbhu
这本书写的不错，讲解的详细，代码也看的懂。我只是把里边的c代码人工翻译成了AS3代码。里边还有好多算法，有时间都用AS3实现一遍。
Fourier.as代码如下：
package
{
/**
* @date 2015/7/27 15:50
* @author hanyeah
*/
public class Fourier
{
/**
*
*/
public function Fourier()
{
//构造函数
throw(new Error("静态类，直接通过类名来调用类中的方法，如：Fourier.FFT(fre,fim,r)"));
}
/**
* 一维离散傅里叶变换，直接由定义计算。
* @param fre 输入信号的实部。
* @param fim 输入信号的虚部。
* @param dir 1/-1:正变换/反变换，默认是1。
* @return Object { "re":gre, "im":gim }。
*/
public static function DFT(fre:Vector.<Number>, fim:Vector.<Number>, dir:int = 1):Object
{
var n = fre.length;
var gre:Vector.<Number> = new Vector.<Number>(n);
var gim:Vector.<Number> = new Vector.<Number>(n);
var i:int, j:int;
var sqrtN:Number = Math.sqrt(n);
for (i = 0; i < n; i++)
{
gre[i] = 0;
gim[i] = 0;
var arg = -dir * 2 * Math.PI * i / n;
for (j = 0; j < n; j = j++)
{
var cos:Number = Math.cos(j * arg);
var sin:Number = Math.sin(j * arg);
gre[i] += fre[j] * cos - fim[j] * sin;
gim[i] += fre[j] * sin + fim[j] * cos;
}
if (dir == 1)
{
//do nothing;
}
else
{
gre[i] = gre[i] / sqrtN;
gim[i] = gim[i] / sqrtN;
}
}
return {"re": gre, "im": gim};
}
/**
* 一维快速傅里叶变换。
* @param fre 输入信号的实部。
* @param fim 输入信号的虚部。
* @param r 信号长度应该是2^r
* @return Object { "re":gre, "im":gim }。
*/
public static function FFT(fre:Vector.<Number>, fim:Vector.<Number>, r:int):Object
{
var n:int = 1 << r;
var hn:int = n >> 1;
var wre:Vector.<Number> = new Vector.<Number>(hn);
var wim:Vector.<Number> = new Vector.<Number>(hn);
var gre:Vector.<Number> = new Vector.<Number>(n);
var gim:Vector.<Number> = new Vector.<Number>(n);
var lre:Vector.<Number>;
var lim:Vector.<Number>;
var i:int, j:int, k:int, p:int;
var sqrtN:Number = Math.sqrt(n);
//计算加权系数
for (i = 0; i < hn; i++)
{
var angle:Number = -2 * Math.PI * i / n;
wre[i] = Math.cos(angle);
wim[i] = Math.sin(angle);
}
//蝶形运算
for (k = 0; k < r; k++)
{
for (j = 0; j < (1 << k); j++)
{
var bfsize:int = 1 << (r - k);
var hbfsize:int = bfsize >> 1;
for (i = 0; i < hbfsize; i++)
{
p = j * bfsize;
var i1:int = i + p;
var i2:int = i + p + hbfsize;
gre[i1] = fre[i1] + fre[i2];
gim[i1] = fim[i1] + fim[i2];
var kk:int = i << k;
var xre:Number = fre[i1] - fre[i2];
var xim:Number = fim[i1] - fim[i2];
gre[i2] = xre * wre[kk] - xim * wim[kk];
gim[i2] = xre * wim[kk] + xim * wre[kk];
}
}
//交换
lre = fre;
lim = fim;
fre = gre;
fim = gim;
gre = lre;
gim = lim;
}
//重新排序
for (j = 0; j < n; j++)
{
p = 0;
for (i = 0; i < r; i++)
{
if (j & (1 << i))
{
p += 1 << (r - i - 1);
}
}
gre[j] = fre[p]/sqrtN;
gim[j] = fim[p]/sqrtN;
}
return {"re": gre, "im": gim};
}
/**
* 一维快速傅里叶变换,反变换。
* @param fre 输入信号的实部。
* @param fim 输入信号的虚部。
* @param r 信号长度应该是2^r
* @return Object { "re":gre, "im":gim }。
*/
public static function IFFT(fre:Vector.<Number>, fim:Vector.<Number>, r:int):Object
{
var n:int = 1 << r;
var i:int;
for (i = 0; i < n; i++)
{
fim[i] = -fim[i];
}
var g:Object = FFT(fre, fim, r);
for (i = 0; i < n; i++)
{
g.re[i] = g.re[i];
g.im[i] = -g.im[i];
}
return g;
}
/**
* 二维快速傅里叶变换。
* @param fre 输入信号的实部。
* @param fim 输入信号的虚部。
* @param wr 宽=2^wr。
* @param hr 高=2^hr。
* @return Object { "re":gre, "im":gim }。
*/
public static function FFT2(fre:Vector.<Number>, fim:Vector.<Number>, wr:int, hr:int):Object
{
var n:int = 1 << wr;
var m:int = 1 << hr;
var gre:Vector.<Number> = new Vector.<Number>(m * n);
var gim:Vector.<Number> = new Vector.<Number>(m * n);
var i:int, j:int, index:int;
var lre:Vector.<Number>;
var lim:Vector.<Number>;
var g:Object;
lre = new Vector.<Number>(n);
lim = new Vector.<Number>(n);
for (i = 0; i < m; i++)
{
for (j = 0; j < n; j++)
{
index = i * n + j;
lre[j] = fre[index];
lim[j] = fim[index];
}
g = FFT(lre, lim, wr);
for (j = 0; j < n; j++)
{
index = i * n + j;
gre[index] = g.re[j];
gim[index] = g.im[j];
}
}
lre = new Vector.<Number>(m);
lim = new Vector.<Number>(m);
for (j = 0; j < n; j++)
{
for (i = 0; i < m; i++)
{
index = i * n + j;
lre[i] = gre[index];
lim[i] = gim[index];
}
g = FFT(lre, lim, hr);
for (i = 0; i < m; i++)
{
index = i * n + j;
gre[index] = g.re[i];
gim[index] = g.im[i];
}
}
/*
var sqrtMul:Number = Math.sqrt( m * n);
for (i = 0; i < m; i++)
{
for (j = 0; j < n; j++)
{
index = i * n + j;
gre[index] = gre[index] / sqrtMul;
gre[index] = -gre[index] / sqrtMul;
}
}
*/
return {"re": gre, "im": gim};
}
/**
* 二维快速傅里叶变换,反变换。
* @param fre 输入信号的实部。
* @param fim 输入信号的虚部。
* @param wr 宽=2^wr。
* @param hr 高=2^hr。
* @return Object { "re":gre, "im":gim }。
*/
public static function IFFT2(fre:Vector.<Number>, fim:Vector.<Number>, wr:int, hr:int):Object
{
var n:int = 1 << wr;
var m:int = 1 << hr;
var i:int, j:int, index:int;
for (i = 0; i < m; i++)
{
for (j = 0; j < n; j++)
{
index = i * n + j;
fim[index] = -fim[index];
}
}
var g:Object = FFT2(fre, fim, wr, hr);
for (i = 0; i < m; i++)
{
for (j = 0; j < n; j++)
{
index = i * n + j;
g.re[index] = g.re[index];
g.im[index] = -g.im[index];
}
}
return g;
}
//end function
}
}我们知道，离散傅里叶变换，F(i)=(1/N)*求和项，反变换是f(x)=求和项；物理上，通常是F(i)=（根号(1/N)）*求和项，f(x)=（根号(1/N)）*求和项。
参考书中采用前者，上面的代码采用后者。
用上面的代码，两次傅里叶变换，会的到图像的等大的倒立的实像。听着很耳熟吧，没错，就是透镜成像。不难推出，4次傅里叶变化，图像就会还原。（这就是为什么正负变换都要乘以（根号(1/N)）的原因了，符合物理规律，功率保持不变）
同样，也用js实现了一下，js实现代码如下：
function DFT(fre,fim,dir){
var n=fre.length;
var gre=[];
var gim=[];
for(var i=0;i<n;i++){
gre[i]=0;
gim[i]=0;
var arg=-dir*2*Math.PI*i/n;
for(var j=0;j<n;j++){
var cos=Math.cos(j*arg);
var sin=Math.sin(j*arg);
gre[i]+=fre[j]*cos-fim[j]*sin;
gim[i]+=fre[j]*sin+fim[j]*cos;
}
if(dir==1){
//do nothing
}
else{
gre[i]=gre[i]/n;
gim[i]=gim[i]/n;
}
}
return {re:gre,im:gim};
}
function FFT(fre,fim,r){
var n=1<<r;
var wre=[];
var wim=[];
var hn=n>>1;
var gre=[];
var gim=[];
var lre=[];
var lim=[];
//计算加权系数
for(var i=0;i<hn;i++){
var angle=-2*Math.PI*i/n;
wre[i]=Math.cos(angle);
wim[i]=Math.sin(angle);
}
//蝶形运算
for(var k=0;k<r;k++){
for(var j=0;j<(1<<k);j++){
var bfsize=1<<(r-k);
var hbfsize=bfsize>>1;
for(var i=0;i<hbfsize;i++){
var p=j*bfsize;
var i1=i+p;
var i2=i+p+hbfsize;
gre[i1]=fre[i1]+fre[i2];
gim[i1]=fim[i1]+fim[i2];
var kk=i<<k;
var xre=fre[i1]-fre[i2];
var xim=fim[i1]-fim[i2];
gre[i2]=xre*wre[kk]-xim*wim[kk];
gim[i2]=xre*wim[kk]+xim*wre[kk];
}
}
//交换
lre=fre;
lim=fim;
fre=gre;
fim=gim;
gre=lre;
gim=lim;
}
//重新排序
for(var j=0;j<n;j++){
var p=0;
for(var i=0;i<r;i++){
if(j&(1<<i)){
p+=1<<(r-i-1);
}
}
gre[j]=fre[p];
gim[j]=fim[p];
}
return {re:gre,im:gim};
}
function IFFT(fre,fim,r){
var n=1<<r;
for(var i=0;i<n;i++){
fim[i]=-fim[i];
}
var g=FFT(fre,fim,r);
for(var i=0;i<n;i++){
g.re[i]=g.re[i]/n;
g.im[i]=-g.im[i]/n;
}
return g;
}
function FFT2(fre,fim,wr,hr){
var n=1<<wr;
var m=1<<hr;
var gre=[];
var gim=[];
for(var i=0;i<m;i++){
var lre=[];
var lim=[];
var index;
for(var j=0;j<n;j++){
index=i*n+j;
lre[j]=fre[index];
lim[j]=fim[index];
}
var g=FFT(lre,lim,wr);
for(j=0;j<n;j++){
index=i*n+j;
gre[index]=g.re[j];
gim[index]=g.im[j];
}
}
for(var j=0;j<n;j++){
var lre=[];
var lim=[];
var index;
for(var i=0;i<m;i++){
index=i*n+j;
lre[i]=gre[index];
lim[i]=gim[index];
}
var g=FFT(lre,lim,hr);
for(i=0;i<m;i++){
index=i*n+j;
gre[index]=g.re[i];
gim[index]=g.im[i];
}
}
return {re:gre,im:gim};
}
function IFFT2(fre,fim,wr,hr){
var n=1<<wr;
var m=1<<hr;
var g;
var mul=m*n;
for(var i=0;i<m;i++){
for(var j=0;j<n;j++){
index=i*n+j;
fim[index]=-fim[index];
}
}
var g=FFT2(fre,fim,wr,hr);
for(var i=0;i<m;i++){
for(var j=0;j<n;j++){
index=i*n+j;
g.re[index]=g.re[index]/mul;
g.im[index]=-g.im[index]/mul;
}
}
return g;
}看一个demo吧。（不要点DFT，会卡15秒钟，然后flashplayer可能会崩溃，不是BUG，是因为DFT太慢了）
发表评论:
◎欢迎参与讨论，请在这里发表您的看法、交流您的观点。