1
2
3 """RMS or peak amplitude time series."""
4
5
6 import sys
7 import numpy
8 from gmisclib import die
9 from gpk_voicing import RMSlib
10
11
12
13
14
15 DR = 'rms'
16 DTSMOOTH = 0.015
17 DT = 0.01
18 HPF = 40.0
19 PKWDT = 0.020
20
21
22
23 if __name__ == '__main__':
24 import gpkimgclass
25 arglist = sys.argv[1:]
26 columns = None
27 signalfile = None
28 outfile = None
29 extrahdr = {}
30 while len(arglist)>0 and arglist[0][0]=='-':
31 arg = arglist.pop(0)
32 if arg == '-c':
33 if columns == None:
34 columns = []
35 tmp = arglist.pop(0)
36 try:
37 col = int(tmp)
38 except ValueError:
39 col = tmp
40 columns.append(col)
41 elif arg == '-pkpk':
42 DR = 'pkpk'
43 elif arg == '-f':
44 signalfile = arglist.pop(0)
45 die.note("signalfile", signalfile)
46 elif arg == '-o':
47 outfile = arglist.pop(0)
48 elif arg == '-dt':
49 DT = float(arglist.pop(0))
50 assert DT>0.0, 'Silly time step.'
51 elif arg == '-hpf':
52 HPF = float(arglist.pop(0))
53 assert HPF > 0.0, "Silly highpass frequency"
54 elif arg == '-write':
55 extrahdr = dict( [q.strip().split('=', 1)
56 for q in arglist.pop(0).split(';') ]
57 )
58 else:
59 die.die("Unrecognized flag: %s" % arg)
60 if len(arglist) > 0:
61 rootname = arglist.pop(0)
62 if signalfile == None:
63 signalfile = rootname + '.d'
64 if outfile == None:
65 outfile = '%s.%s.dat' % (rootname, DR)
66 if outfile == None:
67 outfile = '%s.dat' % DR
68 signal = gpkimgclass.read(signalfile)
69 nch = signal.d.shape[1]
70 if columns == None:
71 columns = range(nch)
72 die.note("columns", str(columns))
73 for c in columns:
74 try:
75 signal.column(c)
76 except KeyError:
77 die.info("File has %d columns" % signal.n[1])
78 die.die("Bad column: %s" % str(c))
79
80 tmp = []
81 t00 = None
82 for i in columns:
83 tmpR, t0 = RMSlib.compute_intensity(signal.column(i),
84 signal.dt(),
85 DT, DR, HPF,
86 DTSmooth=DTSMOOTH)
87 tmp.append( tmpR )
88 if t00 == None:
89 t00 = t0
90 else:
91 assert abs(t00-t0)/DT < 0.001
92 n = len(tmp[0])
93 o = numpy.zeros((n,))
94 for t in tmp:
95 numpy.add(o, t, o)
96 o = o / nch
97
98 hdr = signal.hdr.copy()
99 hdr['CDELT2'] = DT
100 hdr['CRPIX2'] = 0
101 hdr['CRVAL2'] = signal.start() + t0
102 hdr['CDELT1'] = 1
103 hdr['TTYPE1'] = DR
104 hdr['BITPIX'] = -32
105 hdr.update( extrahdr )
106 oo = gpkimgclass.gpk_img(hdr, numpy.transpose( [o]))
107 oo.write(outfile)
108