这一篇文章中说明了用"二维卷积"的方法进行滤波/降噪( 二维卷积滤波 )。本文主要介绍另一种滤波的方法:在二维傅里叶变换后的" 频振谱 "中,用" 滤波器 "进行滤波,并对比这两种滤波方法的优劣。
滤波器没那么复杂,就是一个函数式而已,只不过这个函数式有一些特别的功能。本文选用的是" 巴特沃斯滤波器 ";图像的噪声还是" 高斯噪声 "和" 椒盐噪声 "。巴特沃斯滤波器的函数式为:
其中 是截止频率(高于这个频率值,就被滤掉了), 是阶次, 是" 中心化频振图( 中心化参考这里 ) "中各点" 距中心点的距离 "。非常简单的一个函数。而且注意到: 说明这就是在频域内的一个函数,所以它的用法就是直接和" 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的值:
地球物理数据处理基础
地球物理数据处理基础
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)