Hello,
I'm trying to calculate ppsd values on an infrasound channel. The
same exact code works flawlessly on typical 100 Hz seismic data and
has a ppsd similar to what would be expected. Strangely, the ppsd
comes out to be exactly 50 dB everywhere below the Nyquist. If I
filter the data using a bandpass filter, then the ppsd has values that
are not 50 dB, but only outside the passband of the filter. The data
looks normal when plotting, and when removing the instrument response
with simulate. I've tried the ppsd with is_rotational_data=True
(which keeps the code from using the poles and zeroes) and the same
results occur. Here is some output:
print tr.stats
network: HV
station: AHUD
location: 03
channel: BDF
starttime: 2012-10-02T00:00:00.010000Z
endtime: 2012-10-02T01:00:01.010000Z
sampling_rate: 40.0
delta: 0.025
npts: 144041
calib: 1.0
paz
Out[117]:
{'digitizer_gain': 1.0,
'gain': 0.291944,
'poles': [(-1190+0j), (-0.157+3e-06j), (-0.157-3e-06j), (-0.157+0j)],
'seismometer_gain': 0.0097,
'sensitivity': 6104.47,
'zeros': [0j, 0j, 0j, (-4080+0j)]}
len(ppsd.times)
Out[122]: 25
print ppsd.get_percentile()
(array([ 3.69553048e-02, 4.03000457e-02, 4.39475115e-02,
4.79251011e-02, 5.22626933e-02, 5.69928712e-02,
6.21511667e-02, 6.77763279e-02, 7.39106097e-02,
8.06000914e-02, 8.78950229e-02, 9.58502022e-02,
1.04525387e-01, 1.13985742e-01, 1.24302333e-01,
1.35552656e-01, 1.47821219e-01, 1.61200183e-01,
1.75790046e-01, 1.91700404e-01, 2.09050773e-01,
2.27971485e-01, 2.48604667e-01, 2.71105312e-01,
2.95642439e-01, 3.22400366e-01, 3.51580092e-01,
3.83400809e-01, 4.18101547e-01, 4.55942970e-01,
4.97209334e-01, 5.42210623e-01, 5.91284878e-01,
6.44800731e-01, 7.03160183e-01, 7.66801617e-01,
8.36203093e-01, 9.11885939e-01, 9.94418668e-01,
1.08442125e+00, 1.18256976e+00, 1.28960146e+00,
1.40632037e+00, 1.53360323e+00, 1.67240619e+00,
1.82377188e+00, 1.98883734e+00, 2.16884249e+00,
2.36513951e+00, 2.57920292e+00, 2.81264073e+00,
3.06720647e+00, 3.34481237e+00, 3.64754376e+00,
3.97767467e+00, 4.33768499e+00, 4.73027902e+00,
5.15840585e+00, 5.62528147e+00, 6.13441294e+00,
6.68962474e+00, 7.29508751e+00, 7.95534934e+00,
8.67536997e+00, 9.46055804e+00, 1.03168117e+01,
1.12505629e+01, 1.22688259e+01, 1.33792495e+01,
1.45901750e+01, 1.59106987e+01, 1.73507399e+01,
1.89211161e+01, 2.06336234e+01, 2.25011259e+01,
2.45376518e+01, 2.67584990e+01, 2.91803501e+01,
3.18213974e+01, 3.47014799e+01, 3.78422322e+01,
4.12672468e+01, 4.50022517e+01, 4.90753035e+01,
5.35169980e+01, 5.83607001e+01, 6.36427947e+01,
6.94029598e+01, 7.56844643e+01, 8.25344936e+01,
9.00045035e+01, 9.81506070e+01, 1.07033996e+02,
1.16721400e+02, 1.27285589e+02, 1.38805920e+02,
1.51368929e+02, 1.65068987e+02, 1.80009007e+02,
1.96301214e+02, 2.14067992e+02, 2.33442800e+02,
2.54571179e+02, 2.77611839e+02, 3.02737857e+02,
3.30137974e+02, 3.60018014e+02, 3.92602428e+02,
4.28135984e+02, 4.66885601e+02, 5.09142358e+02,
5.55223678e+02]), array([-147. , -75.5, -55.5, -50. ,
-50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ,
-50. , -50. , -50. , -50. , -50. , -50. , -50. , -50. ]))
Here is the code:
from obspy.earthworm import Client
from obspy.core import UTCDateTime
from obspy.core import Stream
from obspy.xseed import Parser
from obspy.signal import PPSD, cornFreq2Paz, pazToFreqResp, cosTaper
import matplotlib.pyplot as plt
import numpy as np
import os
client = Client("hvo-wws0.wr.usgs.gov", 16022)
response = client.availability(network='HV', station='AHUD',
channel='BDF', location='03')
seedvol = 'ahud.seed'
parser = Parser(seedvol)
t1 = UTCDateTime(2012,10,1,0,0,0,0)
t2 = UTCDateTime(2012,10,2,0,0,0,0)
tNow = t1
while (t2-tNow >= 0):
# Get data
st = client.getWaveform(response[0][0], response[0][1],
response[0][2], response[0][3], tNow, tNow+3601)
if len(st) != 0:
tr = st[0]
paz = parser.getPAZ(tr.id)
if t1 == tNow:
# Initialize ppsd
ppsd = PPSD(tr.stats, paz, is_rotational_data=False)
ppsd.add(st)
tNow = tNow + 3600
# Plot ppsd
ppsd.plot(show_noise_models=False, show_percentiles=False )
Any insight into what I'm doing wrong would be much appreciated,
wes
Weston Thelen
Seismic Network Manager
Hawaiian Volcano Observatory
USGS HVO Bldg. 336
Hawaii National Park, HI 96718-0051
wthelen@usgs.gov
(808) 967-8803
