python,地震波形、时频图、频谱图计算和显示软件

在地震波形分析时，时频图可以辅助识别信号，频谱图可以观察地震信号的频率分布范围，发现波形的优势频率，是地震类型判断和科学研究的重要工具。

下面这段基于python开发的小程序只需要给定地震数据的SAC格式文件，自动进行傅丽叶变换，画出原始波形、时频图、频谱图，非常方便地实现了地震信号分析。

import struct
import pylab
import numpy as np
from scipy.fftpack import fft,ifft


class sacfile_wave:
    def read(self, sFile):
        # sFile='e:\\10.sac'

        f = open(sFile, 'rb')
        hdrBin = f.read(632)

        sfmt = 'f' * 70 + 'I ' * 40 + '8s ' * 22 + '16s';
        hdrFmt = struct.Struct(sfmt)
        self.m_header = hdrFmt.unpack(hdrBin)
        self.dt = self.m_header[0]

        npts = int(self.m_header[79])
        fmt_data = 'f' * npts
        dataFmt = struct.Struct(fmt_data)
        dataBin = f.read(4 * npts)
        f.close()
        self.m_data = dataFmt.unpack(dataBin)

    def draw(self, sImageFile):
        npts = len(self.m_data)
        xd = range(1, npts + 1)
        dt = self.dt
        NFFT = 250
        Fs = int(1.0 / dt)
        x0 = np.arange(0, npts/Fs, np.round(dt,2))

        pylab.subplot(311)
        pylab.plot(x0, self.m_data, linewidth=.6)
        # pylab.xlabel('Time / second')
        # pylab.xlabel('时间 / 秒', fontproperties = 'FangSong', fontsize = 20)
        pylab.subplot(312)
        Pxx, freqs, bins, imm = pylab.specgram(self.m_data, NFFT=NFFT, Fs=Fs,
                                                noverlap=NFFT * 19 / 20,cmap=pylab.get_cmap('jet'))


        pylab.subplot(313) # FFT快速傅丽叶变换，显示频谱曲线
        n = len(self.m_data)
        k = np.arange(n)
        T = n / Fs
        frq = k / T
        frq1 = frq[range(int(n / 2))]
        YY = np.fft.fft(self.m_data)
        Y = np.fft.fft(self.m_data) / n
        Y1 = Y[range(int(n / 2))]
        Y1[0] = Y1[1]
        pylab.plot(frq1, abs(Y1))

        pylab.savefig(sImageFile, dpi = 1200)
        pylab.show()


    def exportAsc(self, sAscFile):
        f2 = open(sAscFile, "wt")
        sdataAsc = [str(x) for x in self.m_data]
        sDataAsc = '\n'.join(sdataAsc)
        f2.writelines(sDataAsc)
        f2.close()



if __name__ == "__main__":
    print('输入文件名必须为1.sac')
    print('输出文件名为1.png')
    sacfile = '1.sac'
    sac = sacfile_wave()
    sac.read(sacfile)
    sac.draw("1.png")