diff --git a/msnoise/api.py b/msnoise/api.py index a45091b7..bb86bc34 100644 --- a/msnoise/api.py +++ b/msnoise/api.py @@ -226,7 +226,7 @@ def get_filters(session, all=False): def update_filter(session, ref, low, mwcs_low, high, mwcs_high, - rms_threshold, mwcs_wlen, mwcs_step, used): + rms_threshold, mwcs_wlen, mwcs_step, used, whitening_method): """Updates or Insert a new Filter in the database. .. seealso:: :class:`msnoise.msnoise_table_def.Filter` @@ -255,6 +255,8 @@ def update_filter(session, ref, low, mwcs_low, high, mwcs_high, :param mwcs_step: Step (in seconds) of the windowing procedure in MWCS :type used: bool :param used: Is the filter activated for the processing + :type whitening_method: str + :param whitening_method: Whitening method to use """ filter = session.query(Filter).filter(Filter.ref == ref).first() if filter is None: @@ -267,6 +269,7 @@ def update_filter(session, ref, low, mwcs_low, high, mwcs_high, filter.mwcs_wlen = mwcs_wlen filter.mwcs_step = mwcs_step filter.used = used + filter.whitening_method = whitening_method session.add(filter) else: filter.low = low @@ -277,6 +280,7 @@ def update_filter(session, ref, low, mwcs_low, high, mwcs_high, filter.mwcs_wlen = mwcs_wlen filter.mwcs_step = mwcs_step filter.used = used + filter.whitening_method = whitening_method session.commit() return diff --git a/msnoise/move2obspy.py b/msnoise/move2obspy.py index 08d7db4c..7ff43969 100644 --- a/msnoise/move2obspy.py +++ b/msnoise/move2obspy.py @@ -6,6 +6,8 @@ import scipy.optimize import scipy.signal from obspy.signal.invsim import cosine_taper +from obspy.signal.konnoohmachismoothing \ + import konno_ohmachi_smoothing_window as ko_window from scipy.fftpack.helper import next_fast_len try: @@ -71,7 +73,7 @@ def myCorr(data, maxlag, plot=False, nfft=None): return corr -def whiten(data, Nfft, delta, freqmin, freqmax, plot=False): +def whiten(data, Nfft, delta, freqmin, freqmax, plot=False, method='brutal'): """This function takes 1-dimensional *data* timeseries array, goes to frequency domain using fft, whitens the amplitude of the spectrum in frequency domain between *freqmin* and *freqmax* @@ -89,10 +91,33 @@ def whiten(data, Nfft, delta, freqmin, freqmax, plot=False): :param freqmax: The upper frequency bound :type plot: bool :param plot: Whether to show a raw plot of the action (default: False) + :type method: str + :param method: Whitening method to use, currently supported "brutal", + "smoothing", "smooth_whitening", "none". :rtype: :class:`numpy.ndarray` :returns: The FFT of the input trace, whitened between the frequency bounds -""" + """ + + def smoothen(spectrum, bandwidth, low, high): + """Smoothens the given spectrum between low and high using a + Konno Ohmachi window. Returns an array with the smoothened part only""" + + # half-width of the sliding window + wl = 90 + # create the filter + ko_filter = ko_window(freqVec[:2*wl+1], freqVec[wl], + bandwidth, normalize=False) + + # add zeros to the low end of the spectrum if filter doesn't fit + if low < wl: + offset = wl - low + spectrum = np.append(np.zeros(offset), spectrum) + low += offset + high += offset + + return np.array([np.dot(spectrum[x - wl:x + wl + 1], ko_filter) \ + for x in range(low, high)]) if plot: plt.subplot(411) @@ -109,8 +134,8 @@ def whiten(data, Nfft, delta, freqmin, freqmax, plot=False): if low <= 0: low = 1 - porte1 = J[0] - porte2 = J[-1] + p1 = J[0] + p2 = J[-1] high = J[-1] + Napod if high > Nfft / 2: high = int(Nfft // 2) @@ -124,33 +149,74 @@ def whiten(data, Nfft, delta, freqmin, freqmax, plot=False): plt.xlim(0, max(axis)) plt.title('FFTRawSign') - # Left tapering: + if method == 'brutal': + # Left tapering: amplitude spectrum is the cosine taper itself + FFTRawSign[low:p1] = np.cos( + np.linspace(np.pi / 2., np.pi, p1 - low)) ** 2 * np.exp( + 1j * np.angle(FFTRawSign[low:p1])) + + # Pass band: amplitude spectrum is 1 + FFTRawSign[p1:p2] = np.exp(1j * np.angle(FFTRawSign[p1:p2])) + + # Right tapering: amplitude spectrum is the cosine taper itself + FFTRawSign[p2:high] = np.cos( + np.linspace(0., np.pi / 2., high - p2)) ** 2 * np.exp( + 1j * np.angle(FFTRawSign[p2:high])) + + elif method in ('smoothing', 'smooth_whitening', 'none'): + if method == 'smoothing': + # Pass band: a smoothened spectrum for the pass band + FFTRawSign[low:high] = (smoothen(np.abs(FFTRawSign), 10, low, high) * + np.exp(1j * np.angle(FFTRawSign[low:high]))) + + elif method == 'smooth_whitening': + # compute a smooth amplitude spectrum + smooth = smoothen(np.abs(FFTRawSign), 10, low, high) + waterlevel = 0.001 * np.max(smooth) + # add waterlevel + smooth[smooth < waterlevel] = waterlevel + + # Pass band: the amplitudes are whitened using a smooth spectrum + FFTRawSign[low:high] = ((np.abs(FFTRawSign[low:high]) / smooth) * + np.exp(1j * np.angle(FFTRawSign[low:high]))) + + # Left tapering: amplitude spectrum is multiplied by cosine taper + FFTRawSign[low:p1] = (np.abs(FFTRawSign[low:p1]) * + np.cos(np.linspace(np.pi / 2., np.pi, p1 - low)) ** 2 * + np.exp(1j * np.angle(FFTRawSign[low:p1]))) + + # Right tapering: amplitude spectrum is multiplied by cosine taper + FFTRawSign[p2:high] = (np.abs(FFTRawSign[p2:high]) * + np.cos(np.linspace(0., np.pi / 2., high - p2)) ** 2 * + np.exp(1j * np.angle(FFTRawSign[p2:high]))) + + else: + raise ValueError('Unknown whitening type %s' % method) + + # set lower part of spectrum to zero FFTRawSign[0:low] *= 0 - FFTRawSign[low:porte1] = np.cos( - np.linspace(np.pi / 2., np.pi, porte1 - low)) ** 2 * np.exp( - 1j * np.angle(FFTRawSign[low:porte1])) - # Pass band: - FFTRawSign[porte1:porte2] = np.exp(1j * np.angle(FFTRawSign[porte1:porte2])) - # Right tapering: - FFTRawSign[porte2:high] = np.cos( - np.linspace(0., np.pi / 2., high - porte2)) ** 2 * np.exp( - 1j * np.angle(FFTRawSign[porte2:high])) + + # set upper part of spectrum to zero FFTRawSign[high:Nfft + 1] *= 0 # Hermitian symmetry (because the input is real) FFTRawSign[-(Nfft // 2) + 1:] = FFTRawSign[1:(Nfft // 2)].conjugate()[::-1] + # Normalize the amplitude spectrum + FFTRawSign = (np.abs(FFTRawSign) / np.amax(FFTRawSign) * + np.exp(1j * np.angle(FFTRawSign))) + if plot: plt.subplot(413) axis = np.arange(len(FFTRawSign)) plt.axvline(low, c='g') - plt.axvline(porte1, c='g') - plt.axvline(porte2, c='r') + plt.axvline(p1, c='g') + plt.axvline(p2, c='r') plt.axvline(high, c='r') plt.axvline(Nfft - high, c='r') - plt.axvline(Nfft - porte2, c='r') - plt.axvline(Nfft - porte1, c='g') + plt.axvline(Nfft - p2, c='r') + plt.axvline(Nfft - p1, c='g') plt.axvline(Nfft - low, c='g') plt.plot(axis, np.abs(FFTRawSign)) diff --git a/msnoise/msnoise_admin.py b/msnoise/msnoise_admin.py index 04c904e9..27bb12d0 100644 --- a/msnoise/msnoise_admin.py +++ b/msnoise/msnoise_admin.py @@ -182,18 +182,27 @@ def mwcs_step(form, field): if field.data > form.data['mwcs_wlen']: raise ValidationError("'mwcs_step' should be smaller or equal to" " 'mwcs_wlen'") + + def whitening_method(form, field): + whitening_methods = ('brutal', 'smoothing', 'smooth_whitening', 'none') + if field.data not in whitening_methods: + raise ValidationError( + "'whitening_method' should be one of %s" % whitening_methods) form_args = dict( mwcs_low=dict(validators=[mwcs_low]), mwcs_high=dict(validators=[mwcs_high]), high=dict(validators=[high]), mwcs_step=dict(validators=[mwcs_step]), + whitening_method=dict(validators=[whitening_method]), ) column_list = ('ref', 'low', 'mwcs_low', 'mwcs_high', 'high', - 'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used') + 'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used', + 'whitening_method') form_columns = ('low', 'mwcs_low', 'mwcs_high', 'high', - 'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used') + 'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used', + 'whitening_method') def __init__(self, session, **kwargs): # You can pass name and other parameters if you want to diff --git a/msnoise/msnoise_table_def.py b/msnoise/msnoise_table_def.py index 4f4845d2..6ab0d926 100644 --- a/msnoise/msnoise_table_def.py +++ b/msnoise/msnoise_table_def.py @@ -31,6 +31,8 @@ class Filter(Base): :param mwcs_step: Step (in seconds) of the windowing procedure in MWCS :type used: bool :param used: Is the filter activated for the processing + :type whitening_method: str + :param whitening_method: Whitening method to use """ __tablename__ = "filters" @@ -44,6 +46,7 @@ class Filter(Base): mwcs_wlen = Column(Float()) mwcs_step = Column(Float()) used = Column(Boolean(True)) + whitening_method = Column(String(50)) def __init__(self, **kwargs): """""" @@ -55,6 +58,7 @@ def __init__(self, **kwargs): # self.mwcs_wlen = mwcs_wlen # self.mwcs_step = mwcs_step # self.used = used + # self.whitening_method = whitening_method class Job(Base): diff --git a/msnoise/s001configurator.py b/msnoise/s001configurator.py index 6a8618a8..139378d5 100644 --- a/msnoise/s001configurator.py +++ b/msnoise/s001configurator.py @@ -56,6 +56,8 @@ class TFilter(HasTraits): mwcs_wlen = Float mwcs_step = Float used = CBool + whitening_method = Str + # XXX do we have to add whitening_method to View?? traits_view = View('low', 'high', 'rms_threshold', 'used', buttons=['OK', 'Cancel']) @@ -105,6 +107,7 @@ class TStation(HasTraits): StationColumn(name='rms_threshold'), StationColumn(name='mwcs_wlen'), StationColumn(name='mwcs_step'), + StationColumn(name='whitening_method'), ]) config_editor = TableEditor( @@ -225,7 +228,8 @@ def main(): TFilter( ref=f.ref, low=f.low, high=f.high, mwcs_low=f.mwcs_low, mwcs_high=f.mwcs_high, rms_threshold=f.rms_threshold, - mwcs_wlen=f.mwcs_wlen, mwcs_step=f.mwcs_step, used=f.used)) + mwcs_wlen=f.mwcs_wlen, mwcs_step=f.mwcs_step, used=f.used, + whitening_method=f.whitening_method)) configs = [] for name in default.keys(): @@ -263,7 +267,7 @@ def main(): for f in demo.company.filters: update_filter(db, f.ref, f.low, f.mwcs_low, f.high, f.mwcs_high, f.rms_threshold, f.mwcs_wlen, - f.mwcs_step, f.used) + f.mwcs_step, f.used, f.whitening_method) db.close() print("Done !") diff --git a/msnoise/s03compute_cc.py b/msnoise/s03compute_cc.py index 00fe9704..984aacbd 100644 --- a/msnoise/s03compute_cc.py +++ b/msnoise/s03compute_cc.py @@ -84,6 +84,10 @@ example, computing ZZ components for very close by stations (much closer than the wavelength sampled), leading to spatial autocorrelation issues. +The whitening method to use can be set in each filter's configuration as key +``whitening_method``. Currently supported options are "brutal" (MSNoise +default), "smoothing", "smooth_whitening" and "none". + When both traces are ready, the cross-correlation function is computed (see :ref:`mycorr`). The function returned contains data for time lags corresponding to ``maxlag`` in the acausal (negative lags) and causal @@ -390,6 +394,7 @@ def main(): low = float(filterdb.low) high = float(filterdb.high) rms_threshold = filterdb.rms_threshold + whitening_method = filterdb.whitening_method trames2hWb = np.zeros((2, int(nfft)), dtype=np.complex) skip = False @@ -399,7 +404,8 @@ def main(): #logging.debug("Whitening %s" % components) trames2hWb[i] = whiten(tmp[i].data, nfft, dt, low, high, - plot=False) + plot=False, + method=whitening_method) else: #logging.debug("Autocorr %s"%components) tmp[i].filter("bandpass", freqmin=low, diff --git a/msnoise/scripts/msnoise.py b/msnoise/scripts/msnoise.py index 2d7da7ec..5d1bee0f 100644 --- a/msnoise/scripts/msnoise.py +++ b/msnoise/scripts/msnoise.py @@ -173,7 +173,7 @@ def d(path): click.echo('') click.echo('Filters:') print('ID: [low:high] [mwcs_low:mwcs_high] mwcs_wlen mwcs_step' - ' used') + ' used whitening_method') for f in get_filters(db, all=True): data = (f.ref, f.low, @@ -182,8 +182,9 @@ def d(path): f.mwcs_high, f.mwcs_wlen, f.mwcs_step, - ['N', 'Y'][f.used]) - print('%02i: [%.03f:%.03f] [%.03f:%.03f] %.03i %.03i %s' % data) + ['N', 'Y'][f.used], + f.whitening_method) + print('%02i: [%.03f:%.03f] [%.03f:%.03f] %.03i %.03i %s %s' % data) click.echo('') click.echo('Stations:') diff --git a/msnoise/test/tests.py b/msnoise/test/tests.py index 85d668b1..4c287b27 100644 --- a/msnoise/test/tests.py +++ b/msnoise/test/tests.py @@ -63,6 +63,7 @@ def test_004_set_and_get_filters(self): f.mwcs_wlen = 10 f.mwcs_step = 5 f.used = True + f.whitening_method = 'brutal' filters.append(f) f = Filter() f.low = 0.1 @@ -73,16 +74,19 @@ def test_004_set_and_get_filters(self): f.mwcs_wlen = 10 f.mwcs_step = 5 f.used = True + f.whitening_method = 'smooth_whitening' filters.append(f) for f in filters: update_filter(db, f.ref, f.low, f.mwcs_low, f.high, f.mwcs_high, - f.rms_threshold, f.mwcs_wlen, f.mwcs_step, f.used) + f.rms_threshold, f.mwcs_wlen, f.mwcs_step, f.used, + f.whitening_method) dbfilters = get_filters(db) for i, filter in enumerate(dbfilters): for param in ['low', 'mwcs_low', 'high', 'mwcs_high', - 'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used']: + 'rms_threshold', 'mwcs_wlen', 'mwcs_step', 'used', + 'whitening_method']: self.failUnlessEqual(eval("filter.%s" % param), eval("filters[i].%s" % param))