㈠ 頻率域快速數字濾波方法
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值