二维傅里叶变换滤波降噪

这一篇文章中说明了用"二维卷积"的方法进行滤波/降噪( 二维卷积滤波 )。本文主要介绍另一种滤波的方法：在二维傅里叶变换后的" 频振谱 "中，用" 滤波器 "进行滤波，并对比这两种滤波方法的优劣。

滤波器没那么复杂，就是一个函数式而已，只不过这个函数式有一些特别的功能。本文选用的是" 巴特沃斯滤波器 "；图像的噪声还是" 高斯噪声 "和" 椒盐噪声 "。巴特沃斯滤波器的函数式为：

其中 是截止频率(高于这个频率值，就被滤掉了)， 是阶次， 是" 中心化频振图( 中心化参考这里 ) "中各点" 距中心点的距离 "。非常简单的一个函数。而且注意到： 说明这就是在频域内的一个函数，所以它的用法就是直接和" F(u,v) F频域参看这里 "做" 矩阵点乘 "即可达到滤波~

下面我们就实 *** 一下，对一个原始图像做一下低通滤波看看(把图像" 变模糊 "，因为一些信号被滤掉了)，对应的Matlab程序如下：

效果如下，图1是原始图像，图2是低通滤波后图像(记得ifft2回到原始xy空间)：

利用这篇文章中 同样的噪声 (高斯随机噪声、椒盐噪声)，看看用频域滤波效果如何。这里我就设定截止频率 ，阶次 进行" 巴特沃斯低通滤波 "，加噪声后图像如下：

二维傅里叶变换后频域做低通滤波，效果说明：

最后，再补充一个" 加噪声-fft2-滤波-ifft2 "的完整流程的Matlab程序：

本文用到的" zxc.jpg "原图像， 在这里 。

在地球物理数据处理中，经常遇到处理二维实数据的情况。例如在地震勘探中，对面波勘探数据作频散分析解释时，要将时间－空间域的信息转换为频率－波数域频谱；在重磁异常的滤波或转换中，要将空间域的异常f（x，y）转换为波数域F（ω，υ）等。这些分析都需要进行二维的傅里叶变换（FFT）。

根据傅里叶变换的定义，对于连续二维函数f（x，y），其傅里叶变换对为

地球物理数据处理基础

对于离散的二维序列fjk（j＝0，1，…，M－1；k＝0，1，…，N－1），其傅里叶变换为

地球物理数据处理基础

1.二维复序列的FFT算法

对于M条测线，每条测线N个测点，构成复序列yjk（j＝0，1，…，M－1；k＝0，1，…，N－1），根据离散傅里叶公式（8－41），其傅里叶变换为

地球物理数据处理基础

于是，可以分两步套用一维复FFT完成二维复FFT的计算。

（1）沿测线方向计算

对于j＝0，1，…，M－1逐测线套用一维复FFT，执行式（8－43）。定义复数组 则算法为

1）对于j＝0，1，…，M－1，作2）～7）；

2）将yjk输入A1（k），即A1（k）＝yjk（k＝0，1，…，N－1）；

3）计算Wr，存入W（r），即

4）q＝1，2，…，p（p＝log2N），若q为偶数执行6），否则执行5）；

5）k＝0，1，2，…，（2p－q－1）和n＝0，1，2，…，（2q－1－1）循环，作

A2（k2q＋n）＝A1（k2q－1＋n）＋A1（k2q－1＋n＋2p－1）

A2（k2q＋n＋2q－1）＝［A1（k2q－1＋n）－A1（k2q－1＋n＋2p－1）］·W（k2q－1）

至k，n循环结束；

6）k＝0，1，2，…，（2p－q－1）和n＝0，1，2，…，（2q－1－1）循环，作

A1（k2q＋n）＝A2（k2q－1＋n）＋A2（k2q－1＋n＋2p－1）

A1（k2q＋n＋2q－1）＝［A2（k2q－1＋n）－A2（k2q－1＋n＋2p－1）］·W（k2q－1）

至k，n循环结束；

7）q循环结束，若p为偶数，将A1（n）输入到Yjn，否则将A2（n）输入到Yjn（n＝0，1，…，N－1）；

8）j循环结束，得到Yjn（j＝0，1，…，M－1；n＝0，1，…，N－1）。

（2）垂直测线方向计算

对于n＝0，1，…，N－1逐一套用一维复FFT，执行式（8－44）。即

1）对于n＝0，1，…，N－1，作2）～7）；

2）将Yjn输入A1（j），即A1（j）＝Yjn（j＝0，1，…，M－1）；

3）计算Wr存入W（r），即

4）q＝1，2，…，p（p＝log2M），若q为偶数执行6），否则执行5）；

5）j＝0，1，2，…，（2p－q－1）和m＝0，1，2，…，（2q－1－1）循环，作

A2（j2q＋m）＝A1（j2q－1＋m）＋A1（j2q－1＋m＋2p－1）

A2（j2q＋m＋2q－1）＝［A1（j2q－1＋m）－A1（j2q－1＋m＋2p－1）］·W（j2q－1）

至j，m循环结束；

6）j＝0，1，2，…，（2p－q－1）和m＝0，1，2，…，（2q－1－1）循环，作

A1（j2q＋m）＝A2（j2q－1＋m）＋A2（j2q－1＋m＋2p－1）

A1（j2q＋m＋2q－1）＝［A2（j2q－1＋m）－A2（j2q－1＋m＋2p－1）］·W（j2q－1）

至j，m循环结束；

7）q循环结束，若p为偶数，将A1（m）输入到Ymn，否则将A2（m）输入到Ymn（m＝0，1，…，M－1）；

8）n循环结束，得到二维复序列的傅氏变换Ymn（m＝0，1，…，M－1；n＝0，1，…，N－1），

所求得的Ymn是复数值，可以写为

Ymn＝Rmn＋iImn （m＝0，1，…，M－1；n＝0，1，…，N－1）

其中，Rmn，Imn的值也是已知的。

2.二维实序列的FFT算法

对于二维的实序列，我们把其看作是虚部为零的复序列，套用上述的二维复序列FFT方法来求其频谱算法上也是可行的，但势必会增加大量的无功运算。因此，有必要研究二维实序列FFT的实用算法，同一维实序列FFT的实现思路一样，同样把二维实序列按一定的规律构造成二维复序列，调用二维复序列FFT，然后通过分离和加工得到原实序列的频谱。

例如采样区域有2 M条测线，每条测线有N个点，并且M，N都是2的整数幂，需要计算实样本序列xjk（j＝0，1，2，…，2 M－1；k＝0，1，2，…，N－1）的傅氏变换：

地球物理数据处理基础

类似于一维实序列FFT的思想，直接建立下面的二维实序列FFT算法：

（1）将一个二维实序列按偶、奇线号分为两个二维子实序列，分别作为实部和虚部组合为一个二维复序列。即令

地球物理数据处理基础

（2）调用二维复FFT过程，求出yjk的二维傅氏变换Ymn的复数值：

地球物理数据处理基础

式中：Rmn，Imn是Ymn的实部和虚部。

（3）利用Rmn，Imn换算Xmn的值。

前两步容易实现，下面分析第（3）步的实现。

记hjk，gjk的傅氏变换为Hmn，Gmn。根据傅里叶变换的定义，我们导出Xmn与Hmn，Gmn的关系式：

地球物理数据处理基础

式中，Hmn，Gmn为复数，我们用上标r和i表示其实部和虚部，将上式右端实部、虚部分离

地球物理数据处理基础

其中：

地球物理数据处理基础

下面的任务是将Hmn，Gmn各分量与通过二维复FFT求出的Rmn，Imn值联系起来。为此先给出奇、偶分解性质和类似于一维情况的三个二维傅氏变换性质：

（1）奇偶分解性

任何一个正负对称区间定义的函数，均可唯一地分解为如下偶（even）、奇（odd）函数之和：

地球物理数据处理基础

（2）周期性

地球物理数据处理基础

（3）复共轭性

地球物理数据处理基础

现在我们来建立Rmn，Imn与Hmn，Gmn的关系。对Ymn作奇偶分解：

地球物理数据处理基础

根据线性性质

地球物理数据处理基础

对照式（8－54）和式（8－55），得

地球物理数据处理基础

由于hjk，gjk是实函数，根据复共轭性质，上面两式对应的奇偶函数相等。即

地球物理数据处理基础

再由奇偶分解性和周期性，得

地球物理数据处理基础

将式（8－57）代入式（8－50），得

地球物理数据处理基础

再利用Hmn，Gmn周期性及复共轭性，可以得到m＝M/2＋1，…，M－1；n＝0，1，…，N－1的傅氏变换，即

地球物理数据处理基础

将式（8－50）中M，N改为M－m，N－n，并将上式代入，得

地球物理数据处理基础

由式（8－58）、式（8－59）和式（8－61）即可得到原始序列xjk（j＝0，1，…，2M－1；n＝0，1，…，N－1）在m＝0，1，…，M－1；n＝0，1，…，N－1区间的傅氏变换Xmn。

具体二维实序列的FFT算法如下：

（1）令hjk＝x2j，k，gjk＝x2j＋1，k，形成

yjk＝hjk＋igjk （j＝0，1，…，2 M－1；n＝0，1，…，N－1）

（2）调用二维复序列FFT过程，即从两个方向先后调用一维复FFT算法式（8－43）和式（8－44），求得yjk的二维傅氏变换Ymn的复数值：

Ymn＝Rmn＋iImn （m＝0，1，…，M－1；n＝0，1，…，N－1）

（3）用下列公式由Rmn，Imn的值换算Xmn的值：

地球物理数据处理基础

地球物理数据处理基础

---