在地震波形分析时,时频图可以辅助识别信号,频谱图可以观察地震信号的频率分布范围,发现波形的优势频率,是地震类型判断和科学研究的重要工具。
下面这段基于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")