1
2
3 """Computes a measure of the peakiness of the waveform.
4 This is large for pulse trains, and small for hiss-like
5 signals.
6 """
7
8
9 import sys
10 import math
11 from gmisclib import Num
12 from gmisclib import die
13
14 from gpk_voicing import RMSlib
15 from gpk_voicing import power
16
17
18
19 DT = 0.01
20 HPF = 40.0
21 Final_Smooth_Width = 0.030
22 SilenceTune = 0.85
23 DTSMOOTH1 = math.sqrt(Final_Smooth_Width**2 - DT**2/12.0)
24 DTSMOOTH2 = math.sqrt((SilenceTune*Final_Smooth_Width)**2 - DT**2/12.0)
25
26
27
29 EPS = 1e-6
30 avd = RMSlib.compute_intensity(data, dt, Dt, 'avd', Hpf,
31 DTSmooth=DTSMOOTH1)
32 aavd = Num.sum(avd)/avd.shape[0]
33 rms = RMSlib.compute_intensity(data, dt, Dt, 'rms', Hpf,
34 DTSmooth=DTSMOOTH2)
35 return Num.minimum(1.0, Num.log(Num.maximum(1.0, rms/(avd + EPS*aavd))))
36
37
38 if __name__ == '__main__':
39 import gpkimgclass
40 arglist = sys.argv[1:]
41 columns = None
42 signalfile = None
43 outfile = None
44 extrahdr = {}
45 while len(arglist)>0 and arglist[0][0]=='-':
46 arg = arglist.pop(0)
47 if arg == '-c':
48 if columns == None:
49 columns = []
50 tmp = arglist.pop(0)
51 try:
52 col = int(tmp)
53 except ValueError:
54 col = tmp
55 columns.append(col)
56 elif arg == '-f':
57 signalfile = arglist.pop(0)
58 die.note("signalfile", signalfile)
59 elif arg == '-o':
60 outfile = arglist.pop(0)
61 elif arg == '-dt':
62 DT = float(arglist.pop(0))
63 assert DT>0.0, 'Silly time step.'
64 elif arg == '-hpf':
65 HPF = float(arglist.pop(0))
66 assert HPF > 0.0, "Silly highpass frequency"
67 elif arg == '-write':
68 extrahdr = dict( [q.strip().split('=', 1)
69 for q in arglist.pop(0).split(';') ]
70 )
71 else:
72 die.die("Unrecognized flag: %s" % arg)
73 if len(arglist) > 0:
74 rootname = arglist.pop(0)
75 if signalfile == None:
76 signalfile = rootname + '.d'
77 if outfile == None:
78 outfile = '%s.pkns.dat' % rootname
79 if outfile == None:
80 outfile = 'pkns.dat'
81 signal = gpkimgclass.read(signalfile)
82 nch = signal.d.shape[1]
83 if columns == None:
84 columns = range(nch)
85 die.note("columns", str(columns))
86 for c in columns:
87 try:
88 signal.column(c)
89 except KeyError:
90 die.info("File has %d columns" % signal.n[1])
91 die.die("Bad column: %s" % str(c))
92
93 tmp = [
94 peakiness(signal.column(i), signal.dt(), DT, HPF)
95 for i in columns
96 ]
97 n = len(tmp[0])
98 o = Num.zeros((n,), Num.Float)
99 for t in tmp:
100 Num.add(o, t, o)
101 o = o / nch
102
103 hdr = signal.hdr.copy()
104 hdr['CDELT2'] = DT
105 hdr['CRPIX2'] = 0
106 hdr['CRVAL2'] = signal.start()
107 hdr['CDELT1'] = 1
108 hdr['TTYPE1'] = 'peakiness'
109 hdr['BITPIX'] = -32
110 hdr.update( extrahdr )
111 oo = gpkimgclass.gpk_img(hdr, Num.transpose( [o]))
112 oo.write(outfile)
113