1# http://pyrocko.org - GPLv3
2#
3# The Pyrocko Developers, 21st Century
4# ---|P------/S----------~Lg----------
6from ..snuffling import Snuffling, Param, Switch
7import numpy as num
8from matplotlib import cm
11def window(freqs, fc, b):
12 w = num.zeros(len(freqs))
13 if fc == 0.:
14 w[freqs == 0.0] = 1.
15 return w
16 mask = num.logical_and(freqs != 0.0, freqs != fc)
17 T = num.log10(freqs[mask]/fc)*b
18 w[mask] = (num.sin(T)/T)**4
19 w[freqs == fc] = 1.
20 w /= num.sum(w)
21 return w
24class AmpSpec(Snuffling):
25 '''
26 <html>
27 <head>
28 <style type="text/css">
29 body { margin-left:10px };
30 </style>
31 <body>
32 <h1 align='center'>Plot Amplitude Spectrum</h1>
33 <p>
34 When smoothing is activated, a smoothing algorithm is applied as proposed
35 by Konno and Ohmachi, (1998). </p>
36 <p style='font-family:courier'>
37 Konno, K. and Omachi, T. (1998). Ground-motion characteristics
38 estimated from spectral ratio between horizontal and vertical
39 components of microtremor, Bull. Seism. Soc. Am., 88, 228-241
40 </p>
41 </body>
42 </html>
43 '''
44 def setup(self):
45 '''Customization of the snuffling.'''
47 self.set_name('Plot Amplitude Spectrum')
48 self.add_parameter(Switch('Smoothing', 'want_smoothing', False))
49 self.add_parameter(Param('Smoothing band width', 'b', 40., 1., 100.))
50 self.set_live_update(False)
51 self._wins = {}
53 def call(self):
54 '''Main work routine of the snuffling.'''
56 all = []
57 for traces in self.chopper_selected_traces(fallback=True):
58 for trace in traces:
59 all.append(trace)
61 colors = iter(cm.Accent(num.linspace(0., 1., len(all))))
62 if self.want_smoothing:
63 alpha = 0.2
64 additional = 'and smoothing'
65 else:
66 alpha = 1.
67 additional = ''
69 flimits = []
70 alimits = []
71 pblabel = 'Calculating amplitude spectra %s' % additional
72 pb = self.get_viewer().parent().get_progressbars()
73 pb.set_status(pblabel, 0)
74 num_traces = len(all)
75 maxval = float(num_traces)
76 plot_data = []
77 plot_data_supplement = []
78 for i_tr, tr in enumerate(all):
79 val = i_tr/maxval*100.
80 pb.set_status(pblabel, val)
82 tr.ydata = tr.ydata.astype(float)
83 tr.ydata -= tr.ydata.mean()
84 f, a = tr.spectrum()
85 flimits.append(f[1])
86 flimits.append(f[-1])
87 absa = num.abs(a)
88 labsa = num.log(absa)
89 stdabsa = num.std(labsa)
90 meanabsa = num.mean(labsa)
91 lamin, lamax = meanabsa - 5*stdabsa, meanabsa + 5*stdabsa
92 alimits.append(num.exp(lamin))
93 alimits.append(num.exp(lamax))
94 c = next(colors)
95 plot_data.append((f, num.abs(a)))
96 plot_data_supplement.append((c, alpha, '.'.join(tr.nslc_id)))
97 if self.want_smoothing:
98 smoothed = self.konnoohmachi(num.abs(a), f, self.b)
99 plot_data.append((f, num.abs(smoothed)))
100 plot_data_supplement.append((c, 1., '.'.join(tr.nslc_id)))
101 self.get_viewer().update()
103 pb.set_status(pblabel, 100.)
105 fig = self.figure(name='Amplitude Spectra')
106 p = fig.add_subplot(111)
107 args = ('c', 'alpha', 'label')
108 for d, s in zip(plot_data, plot_data_supplement):
109 p.plot(*d, **dict(zip(args, s)))
111 p.set_xscale('log')
112 p.set_yscale('log')
114 amin, amax = num.min(alimits), num.max(alimits)
115 p.set_ylim(amin, amax)
117 fmin, fmax = num.min(flimits), num.max(flimits)
118 p.set_xlim(fmin, fmax)
120 p.set_xlabel('Frequency [Hz]')
121 p.set_ylabel('Counts')
123 handles, labels = p.get_legend_handles_labels()
124 leg_dict = dict(zip(labels, handles))
125 if num_traces > 1:
126 p.legend(list(leg_dict.values()), list(leg_dict.keys()),
127 loc=2,
128 borderaxespad=0.,
129 bbox_to_anchor=((1.05, 1.)))
130 fig.subplots_adjust(right=0.8,
131 left=0.1)
132 else:
133 p.set_title(list(leg_dict.keys())[0], fontsize=16)
134 fig.subplots_adjust(right=0.9,
135 left=0.1)
136 fig.canvas.draw()
138 def konnoohmachi(self, amps, freqs, b=20):
139 smooth = num.zeros(len(freqs), dtype=freqs.dtype)
140 amps = num.array(amps)
141 for i, fc in enumerate(freqs):
142 fkey = tuple((b, fc, freqs[0], freqs[1], freqs[-1]))
143 if fkey in self._wins.keys():
144 win = self._wins[fkey]
145 else:
146 win = window(freqs, fc, b)
147 self._wins[fkey] = win
148 smooth[i] = num.sum(win*amps)
150 return smooth
153def __snufflings__():
154 '''
155 Returns a list of snufflings to be exported by this module.
156 '''
158 return [AmpSpec()]