New parameter: spectral_smooth_width_decades

Fixes #2

Author: claudiodsf

CHANGELOG.md

@@ -22,6 +22,8 @@ Earthquake source parameters from S-wave displacement spectra
     source radius, Brune stress drop, source radius, quality factor
   - Colored console output for log messages! (Not supported on Windows)
 - Processing:
+  - New parameter for setting the width of the spectral smoothing window
+    in terms of frequency decades: `spectral_smooth_width_decades`
   - Compute spectral weights after spectral correction (when a station
     residuals file is specified via `residuals_filepath`)
   - New approach for trace clipping detection (requires just one configuration

sourcespec/configspec.conf

@@ -119,6 +119,20 @@ taper_halfwidth = float(0, 0.5, default=0.05)
 # S-wave window as spectral window length.
 spectral_win_length = float(default=None)
 
+# Spectral smoothing window width in frequency decades
+# (i.e., log10 frequency scale).
+# Example:
+#  spectral_smooth_width_decades=1 means a width of 1 decade
+#  (generally, too large, producing a spectrum which is too smooth).
+#  spectrum(f0) is smoothed using values between f1 and f2, so that
+#  log10(f1)=log10(f0)-0.5 and log10(f2)=log10(f0)+0.5
+#    i.e.,
+#  f1=f0/(10^0.5) and f2=f0*(10^0.5)
+#    or,
+#  f2/f1=10 (1 decade width)
+# Default value of 0.2 is generally a good choice
+spectral_smooth_width_decades = float(default=0.2)
+
 # Residuals file path
 # (a pickle file with the mean residuals per station,
 # used for station correction):

sourcespec/ssp_build_spectra.py

@@ -209,16 +209,29 @@ def _displacement_to_moment(stats, config):
     return coeff
 
 
-def _smooth_spectrum(spec, npts=5):
-    """Smooth spectrum in a log_freq space."""
+def _smooth_spectrum(spec, smooth_width_decades=0.2):
+    """Smooth spectrum in a log10-freq space."""
+    # 1. Generate log10-spaced frequencies
     freq = spec.get_freq()
-    f = interp1d(freq, spec.data, fill_value='extrapolate')
-    freq_log = np.logspace(np.log10(freq[0]), np.log10(freq[-1]))
+    _log_freq = np.log10(freq)
+    # frequencies in logarithmic spacing
+    freq_log = np.logspace(_log_freq[0], _log_freq[-1], len(_log_freq))
     spec.freq_log = freq_log
-    spec.data_log = f(freq_log)
-    spec.data_log = smooth(spec.data_log, window_len=npts)
+    # 2. Reinterpolate data using log10 frequencies
+    # make sure that extrapolation does not create negative values
+    f = interp1d(freq, spec.data, fill_value='extrapolate')
+    _data_log = f(freq_log)
+    _data_log[_data_log <= 0] = np.min(spec.data)
+    # 3. Smooth log10-spaced data points
+    log_df = np.log10(freq_log[1]) - np.log10(freq_log[0])
+    npts = max(1, int(round(smooth_width_decades/log_df)))
+    spec.data_log = smooth(_data_log, window_len=npts)
+    # 4. Reinterpolate to linear frequencies
+    # make sure that extrapolation does not create negative values
     f = interp1d(freq_log, spec.data_log, fill_value='extrapolate')
-    spec.data = f(freq)
+    _data = f(freq)
+    _data[_data <= 0] = np.min(spec.data)
+    spec.data = _data
 
 
 def _build_spectrum(config, trace):
@@ -246,11 +259,11 @@ def _build_spectrum(config, trace):
     # for radiated_energy()
     spec.coeff = coeff
     # smooth
-    _smooth_spectrum(spec, npts=5)
+    _smooth_spectrum(spec, config.spectral_smooth_width_decades)
     return spec
 
 
-def _build_weight(spec, specnoise):
+def _build_weight(config, spec, specnoise):
     weight = spec.copy()
     if specnoise is not None:
         weight.data /= specnoise.data
@@ -259,7 +272,8 @@ def _build_weight(spec, specnoise):
         # The inversion is done in magnitude units,
         # so let's take log10 of weight
         weight.data = np.log10(weight.data)
-        _smooth_spectrum(weight, npts=11)
+        # Weight spectrum is smoothed once more
+        _smooth_spectrum(weight, config.spectral_smooth_width_decades)
         weight.data /= np.max(weight.data)
         # slightly taper weight at low frequencies, to avoid overestimating
         # weight at low frequencies, in cases where noise is underestimated
@@ -277,7 +291,7 @@ def _build_weight(spec, specnoise):
     return weight
 
 
-def _build_weight_st(spec_st, specnoise_st):
+def _build_weight_st(config, spec_st, specnoise_st):
     """Build the weight spectrum."""
     weight_st = Stream()
     spec_ids = set(sp.id[:-1] for sp in spec_st if not sp.stats.ignore)
@@ -287,7 +301,7 @@ def _build_weight_st(spec_st, specnoise_st):
             specnoise_h = _select_spectra(specnoise_st, specid + 'H')[0]
         except Exception:
             continue
-        weight = _build_weight(spec_h, specnoise_h)
+        weight = _build_weight(config, spec_h, specnoise_h)
         weight_st.append(weight)
     return weight_st
 
@@ -323,7 +337,7 @@ def _build_H(spec_st, specnoise_st=None, wave_type='S'):
 
 
 def _check_spectral_sn_ratio(config, spec, specnoise):
-    weight = _build_weight(spec, specnoise)
+    weight = _build_weight(config, spec, specnoise)
     if config.spectral_sn_freq_range is not None:
         sn_fmin, sn_fmax = config.spectral_sn_freq_range
         freqs = weight.get_freq()
@@ -396,7 +410,7 @@ def build_spectra(config, st):
     spec_st = station_correction(spec_st, config)
 
     # build the weight spectrum
-    weight_st = _build_weight_st(spec_st, specnoise_st)
+    weight_st = _build_weight_st(config, spec_st, specnoise_st)
 
     logger.info('Building spectra: done')
     if config.weighting == 'noise':