㈠ 频率域快速数字滤波方法
1.频率域滤波的步骤
(1)对已知地震记录道进行频谱分析
设已知地震记录x(t)，如图9-2-1，包含了有效波s(t)和干扰波n(t)。对此地震记录道进行频谱分析，有效波频率成分在ω1～ω2范围，干扰波在ω3～ω4范围，两者基本上是分开的。见图9-2-2。
(2)设计合适的滤波器
为了滤去干扰波的频谱成分，需要设计一个带通滤波器(图9-2-3)，即在频率ω1～ω2范围|H(ω)|=1，在其他频率范围|H(ω)|=0，这个滤波器可表示如下:
物探数字信号分析与处理技术
(3)进行滤波运算
根据滤波方程，对地震记录道x(t)进行滤波，相当于令x(t)的谱X(ω)与滤波器的频率特性H(ω)相乘，得到 ,相乘后的谱 中消除了干扰波成分，见图9-2-4。
图9-2-1 滤波前地震记录道
图9-2-2 地震记录的频谱
(4)对输出信号谱 进行傅立叶反变换，得到滤波后的输出 ，见图9-2-5。频率滤波的过程可以归纳为以下的数学运算:
物探数字信号分析与处理技术
图9-2-3 带通滤波器
图9-2-4 滤波后地震记录道频谱
图9-2-5 滤波后地震记录道
可见，要进行频率滤波，必须进行两次傅立叶变换，即正、反傅立叶变换。由于采用了快速算法，运算时间大大减少，频率域滤波得到广泛应用。
2.用快速傅立叶变换进行滤波的几个问题
(1)周期性
已知正、反离散傅立叶变换(DFT)公式如式(9-2-2)和(9-2-3)
物探数字信号分析与处理技术
式中N是时间域抽样点个数，也是计算出的频率抽样个数，由连续傅立叶变换过渡到离散傅立叶变换时使用了
物探数字信号分析与处理技术
令
则(9-2-2)和(9-2-3)可以写成一种形式，即
物探数字信号分析与处理技术
(9-2-4)是完成一对DFT的条件，否则就不能进行正、反傅立叶变换的对应计算。可以看出，N就是傅立叶变换的频率抽样点周期，由(9-2-2)式可写出
物探数字信号分析与处理技术
由于
所以
物探数字信号分析与处理技术
(9-2-6)表示X(m)确是以N为频率抽样点数的周期，它表示应用(9-2-2)式计算X(m)时，如果给定的x(n)是N个值，那么只要计算N个X(m)就行了，再多计算就重复
了。见图9-2-6。例如N=50时，
X(0)=X(50)
X(1)=X(51)
……………………
X(49)=X(99)
图9-2-6 频谱图形
在m=0～49一段是计算出的X(m)值，由于以N=50为周期，m=50～99一段与m=0～49是重复的，这就出现了因离散而出现的伪门现象。因此公式(9-2-4)中的参数N，在编制程序时要选择好，应既是x(n)的抽样个数，也是计算X(m)的个数，又是频率抽样个数的周期。它必须满足条件 ，即在编制程序计算X(m)或x(n)时，选择参数Δt，Δf和N必须满足式(9-2-4)。同时周期性告诉我们，在进行快速傅立叶变换时，只要计算N个值就行了，再多计算就重复了。
(2)对称性
对称性是指当x(n)是实数序列时，计算出的频谱满足
物探数字信号分析与处理技术
证明:由式(9-2-2)可知
物探数字信号分析与处理技术
由于
所以得到
物探数字信号分析与处理技术
此式表明，N-m点处的频率对应的频谱值X(N-m)和m点处频率对应的频谱值是共轭关系，X(m)与X(N-m)共轭，其模是相等的
物探数字信号分析与处理技术
例如当N=50时，m=26～50一段的|X(m)|值与m=0～24一段的|X(m)|形状对称。这说明当x(n)取实数序列时，复变谱共轭，振幅谱对称于N/2点处，见图9-12的频谱图形。
3.用FFT算法实现频率域数字滤波的具体方法
1)首先确定理想滤波器的频率特性，起始频率ω1和终止频率ω2，对ω1和ω2要求是在频率间隔的整数倍处;
2)对给定的记录x(n)，(n=0，1，…，N-1)，取N=2m的离散点数做FFT，计算复变谱X(m)(m=1，1，2，…，N-1)，在内存中开辟两个区，一个区存入复变谱的实部，一个区存入复变谱的虚部;
3)按照滤波器的起始频率和频带宽度，对给定的复变谱实部和虚部将要滤去的频率成分充零，得到新的复变谱 的实部和虚部;
物探数字信号分析与处理技术
4)再对 做反傅立叶变换，得到滤波后的地震记录
物探数字信号分析与处理技术
下面举例说明以上步骤。例如，有一时间序列x(n)(n=0，0，1，1，1，1，0，0)，抽样间隔为Δt=10ms，N=8，要求用频率滤波滤去0，12.5Hz分量，求x^(n)。
①对x(n)做正傅立叶变换FFT，见表9-2-1。
表9-2-1 对x(n)做正变换数据
由于时间抽样间隔Δt的倒数和频率抽样间隔Δf相差N倍，所以此处重排时要被N除。由此得到xn的复变谱X(m)，见图9-2-7，由于N=8，Δt=0.01ms，所以Δf=12.5Hz。
②对复变谱X(m)进行频率滤波，为了滤去0、12.5Hz的频率分量，将0、12.5Hz及87.5Hz对应的X(m)值充零，得到滤波后的频谱 ，见表9-2-2。
表9-2-2 滤波后的频谱
根据表9-2-2计算出的振幅谱见图9-2-8。
图9-2-7 x(n)的离散复变谱
图9-2-8 滤波后的振幅谱
③对滤波后的频谱 做反傅立叶变换得到所要求的输出 ，见表9-2-3。
根据表9-2-3作出的振动图形见图9-2-9。
表9-2-3 输出 的数据
④为了验证 与 的对应关系，再对 做一次正傅立叶变换FFT，根据表9-3计算振幅谱作图9-2-9与图9-2-8相同。
由以上例子可以得到以下几点:
a.FFT全部是复数运算;
b.计算出的复变谱以N/2为中心，有共轭关系;
c.频率滤波时，对滤去的频率分量fk充0，同时对fN-k的频率分量也要充0，否则不能进行反变换。
图9-2-9 滤波后的振动图形
㈡ fft后怎么去掉不想要的频率成分
FFT是一种DFT的高效算法，称为快速傅立叶变换（fast Fourier transform）。FFT算法可分为按时间抽取算法和按频率抽取算法，先简要介绍FFT的基本原理。从DFT运算开始，说明FFT的基本原理。DFT的运算为：式中由这种方法计算DFT对于X（K）的每个K值