def getPSD(self, ens):
"""
Calculate and return average power spectral density for ensemble
`ens` using Welch's method
"""
print('in getpsd, ens=', ens, self)
if not self.segyfile:
raise Exception("File not opened")
enstrcs = [t for t in self.segyfile.traces if t.header.original_field_record_number == ens]
psds = [sig.welch(t.data, fs=1/self.dt, scaling='spectrum') for t in enstrcs]
# psds is [(f,psd), (f,psd)...]
freqs = psds[0][0]
psdonly = [p for (f, p) in psds]
psdonly = np.array(psdonly).transpose()
psdavg = np.mean(psdonly, 1)
print(freqs.shape, psdavg.shape)
return(freqs, psdavg)
评论列表
文章目录