「Python与地震工程」滤波

Collapsed building (Generated via OpenAI)

「Python与地震工程」滤波

原理

利用 scipy.signal.filtfilt() 函数实现滤波。
滤波前需要构建滤波器。目前使用 Butterworth 滤波器，用 scipy.signal.butter() 构建。
滤波器截至频率需要输入相对频率（相对于奈奎斯特频率）。

程序代码

import numpy as np
from numpy.fft import fft, fftfreq
import matplotlib.pyplot as plt
from scipy import signal

plt.style.use("ggplot")
plt.rcParams['font.sans-serif'] = ['Microsoft Yahei']
plt.rcParams['axes.unicode_minus'] = False

next2pow = lambda x: \
    int(round(2**np.ceil(np.log(x)/np.log(2.))))

def filter_wave(x, dt, fl=0.1, fh=1, btype="lowpass", order=4):
    '''
    滤波\\
    参数含义：
    x: 信号序列
    dt: 信号的采样时间间隔
    fl: 滤波截止频率（低频）
    fh: 滤波截止频率（高频）
    btype: 滤波器类型
        "lowpass": 低通（保留低于fh的频率成分）
        "highpass": 高通（保留高于fl的频率成分）
        "bandpass": 带通（保留fl~fh之间的频率成分）
        "bandstop": 带阻（过滤fl~fh之间的频率成分）
    '''
    fn = 1.0/dt/2.0  # 奈奎斯特频率
    
    if btype=="lowpass":
        Wn = fh/fn
    elif btype=="highpass":
        Wn = fl/fn
    elif btype=="bandpass" or btype=="bandstop":
        Wn = (fl/fn, fh/fn)
    
    b, a = signal.butter(order, Wn, btype)
    y = signal.filtfilt(b, a, x)
    
    return y

if __name__ == '__main__':
    
    acc0 = np.loadtxt("EQ-S-3.txt") # 读取地震波
    dt = 0.005 # 时间间隔
    n = len(acc0)
    t = np.linspace(0.0,dt*(n-1),n)
    
    # 滤波
    btype = "bandstop" # 滤波类型："lowpass", "highpass", "bandpass", "bandstop"
    fl = 2
    fh = 10
    acc = filter_wave(acc0, dt, fl=fl, fh=fh, btype=btype)
    
    # 绘制时程曲线
    plt.figure(btype+r" 滤波",(10,7.5))
    plt.subplot(3,3,(1,3))
    plt.title("时程曲线")
    plt.plot(t,acc0,label=r"原始地震波")
    plt.plot(t,acc,label=btype+r" 滤波后 "+"（低频 %.1f Hz，高频 %.1f Hz）"%(fl, fh))
    plt.xlabel("时间")
    plt.ylabel("加速度")
    plt.ylim((-1,1))
    plt.grid(True)
    plt.legend()
    
    # 傅里叶变换
    n = len(acc0)
    Nfft = next2pow(n)*2
    af0 = abs(fft(acc0, Nfft))
    af = abs(fft(acc, Nfft))
    f = fftfreq(Nfft, d=dt)
    
    # 绘制傅里叶谱
    plt.subplot(3,3,(4,9))
    plt.title("傅里叶谱")
    plt.plot(f[1:Nfft//2], af0[1:Nfft//2], label="原始地震波")
    plt.plot(f[1:Nfft//2], af[1:Nfft//2], label=btype+r" 滤波后 "+"（低频 %.1f Hz，高频 %.1f Hz）"%(fl, fh))
    plt.xlim((1e-2,1.0/dt/2.0))
    plt.ylim(bottom=0)
    plt.xscale('log')
    plt.grid(True)
    plt.legend()
    plt.xlabel("频率")
    plt.ylabel("傅里叶谱")
    
    plt.tight_layout()
    plt.show()