add new different and less brutal whitening options

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
 

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))

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

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):

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 !")

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,

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:')

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))