有关傅立叶变化资料：
ok了。
http://www.academictutorials.com/graphics/graphics-2d-fourier-transform.asp
http://blog.csdn.net/lnstree/article/details/7412837
了解完傅立叶变化，如果想了解具体的算法实现的话，自行查找相关资料。
傅立叶变化对于图像处理的意义可以简单的描述成，一道阳光通过三棱镜分解成不同波段的光波。而傅立叶变换的功能就是三棱镜。
许多在空间域无法较好解决的问题，转换成频率域就可以很简单的解决。
以下是网上FFT（快速傅立叶变换的VC代码）（基于2的幂数）
/*************************************************************************
*
* 函数名称：
* 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;
}
以下是相应的二维傅立叶变换源码
/*************************************************************************
*
* 函数名称：
* Fourier()
*
* 参数:
* LPSTR lpDIBBits - 指向源DIB图像指针
* LONG lWidth - 源图像宽度（象素数）
* LONG lHeight - 源图像高度（象素数）
*
* 返回值:
* BOOL - 成功返回TRUE，否则返回FALSE。
*
* 说明:
* 该函数用来对图像进行付立叶变换。
*
************************************************************************/
BOOL WINAPI Fourier(LPSTR lpDIBBits, LONG lWidth, LONG lHeight)
{
// 指向源图像的指针
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>(*(lpSrc), 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;
// 判断是否超过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;
}
转成相应C#源码：
// 快速傅立叶变换
public static void FFT(COMPLEX[] t, COMPLEX[] f, int r) // t为时域，f为频域 r为2的幂数
{
long count;
int i, j, k, p, bsize;
COMPLEX[] W;
COMPLEX[] X1;
COMPLEX[] X2;
COMPLEX[] X;
COMPLEX comp;
double angle; // 计算加权时所需角度
count = 1 << r;
W = new COMPLEX[count / 2];
X1 = new COMPLEX[count];
X2 = new COMPLEX[count];
X = new COMPLEX[count];
for (i = 0; i < count / 2; i++)
{
angle = i * Math.PI * 2 / count;
W[i].re = (double)Math.Cos(angle);
W[i].im = -(double)Math.Sin(angle);
}
t.CopyTo(X1, 0);
for (k = 0; k < r; k++)
{
for (j = 0; j < 1 << k; j++)
{
bsize = 1 << (r - k);
for (i = 0; i < bsize/2; i++)
{
p = j * bsize;
X2[i + p].im = X1[i + p].im + X1[i + p + bsize / 2].im;
X2[i + p].re = X1[i + p].re + X1[i + p + bsize / 2].re;
comp.im = X1[i + p].im - X1[i + p + bsize / 2].im;
comp.re = X1[i + p].re - X1[i + p + bsize / 2].re;
X2[i + p + bsize / 2].re = comp.re * W[i * (1 << k)].re -
comp.im * W[i * (1 << k)].im;
X2[i + p + bsize / 2].im = comp.re * W[i * (1 << k)].im +
comp.im * W[i * (1 << k)].re; ;
}
}
X = X1;
X1 = X2;
X2 = X;
}
for (j = 0; j < count; j++)
{
p = 0;
for (i = 0; i < r; i++)
{
if ((j & (1 << i)) != 0)
{
p += 1 << (r - i - 1);
}
}
f[j].re = X1[p].re;
f[j].im = X1[p].im;
}
}// 二维傅立叶变换 public static Bitmap Fourier(Bitmap tp) { // 原图像的宽与高 int w = tp.Width; int h = tp.Height; // 傅立叶变换的实际宽高 long lw = 1; long lh = 1; // 迭代次数 int wp = 0; int hp = 0; long i, j; long n,m; double temp; byte [] ky = new byte[w*h]; ky = ChangeByte(tp); COMPLEX[] t; COMPLEX[] f; while (lw * 2 <= w) { lw *= 2; wp++; } while (lh * 2 <= h) { lh *= 2; hp++; } t = new COMPLEX[lw * lh]; f = new COMPLEX[lw * lh]; COMPLEX[] tw = new COMPLEX[lw]; COMPLEX[] th = new COMPLEX[lw]; for (i = 0; i < lh; i++) { for (j = 0; j < lw; j++) { t[i*lw + j].re = ky[ i*w+ j]; t[i *lw + j].im = 0; } } for (i = 0; i < lh; i++) // 垂直方向傅立叶变换 { Array.Copy(t, i * lw, tw, 0, lw); Array.Copy(f, i * lw, th, 0, lw); FFT(tw, th, wp); // Array.Copy(tw, 0, t, i * lw, lw); Array.Copy(th, 0, f, i * lw, lw); } for (i = 0; i < lh; i++) { for (j = 0; j < lw; j++) { t[j * lh + i].re = f[i * lw + j].re; t[j * lh + i].im = f[i * lw + j].im; } } COMPLEX[] ow = new COMPLEX[lh]; COMPLEX[] oh = new COMPLEX[lh]; for (i = 0; i < lw; i++) { Array.Copy(t, i * lh, ow, 0, lh); Array.Copy(f, i * lh, oh, 0, lh); FFT(ow, oh, hp); //Array.Copy(ow, 0, t, i * lh, lh); oh.CopyTo(f, i * lh); } for (i = 0; i < lh; i++) { for (j = 0; j < lw; j++) { temp = Math.Sqrt(f[j * lh + i].re * f[j * lh + i].re + f[j * lh + i].im * f[j * lh + i].im) / 100; if (temp > 255) { temp = 255; } n = i < lh / 2 ? i + lh / 2 : i - lh / 2; m = j < lw / 2 ? j + lw / 2 : j - lw / 2; ky[n* w + m] = (byte)(temp); } } tp = ChangeBitmap(ky, tp); return tp; }
测试图片：
| 原图 | 傅立叶变换图 |
 |  |
 |  |
2182
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
