600字范文,内容丰富有趣,生活中的好帮手!
600字范文 > python 时频图_python 地震波形 时频图 频谱图计算和显示软件

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

时间:2023-11-05 18:01:22

相关推荐

python 时频图_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")

本内容不代表本网观点和政治立场,如有侵犯你的权益请联系我们处理。
网友评论
网友评论仅供其表达个人看法,并不表明网站立场。