Merge branch 'signal-to-noise' into 'main'

FIX: Normalize spectra by duration

See merge request ghsc/esi/groundmotion-processing!1110

CHANGELOG.md

@@ -1,5 +1,6 @@
 ## main
 
+- Modified the signal-to-noise-ratio calculation to normalize the spectra by duration.
 - Bugfix for how period arrays are defined when using start/stop values in the config file along with the logspace option.
 
 ## 1.2.2 / 2022-11-08

src/gmprocess/metrics/reduction/duration.py

@@ -37,7 +37,7 @@ def __init__(
             smoothing (string):
                 Smoothing type. Default is None.
             interval (list):
-                List of length 2 with the quantiles (0-1) for duration interval
+                List of length 2 with the percentiles (0-100) for duration interval
                 calculation.
             config (dict):
                 Config dictionary.

src/gmprocess/waveform_processing/fft.py

@@ -39,9 +39,11 @@ def compute_and_smooth_spectrum(tr, bandwidth, section, window=None, nfft=None):
         nfft = next_pow_2(tr.stats.npts)
     if window is None:
         window = tr
-
+    lowest_usable_freq = 1 / tr.stats.delta / tr.stats.npts
     spec_raw, freqs_raw = compute_fft(window, nfft)
+    spec_raw[freqs_raw < lowest_usable_freq] = np.nan
     spec_smooth, freqs_smooth = smooth_spectrum(spec_raw, freqs_raw, nfft, bandwidth)
+    spec_smooth[freqs_smooth < lowest_usable_freq] = np.nan
 
     raw_dict = {"spec": spec_raw, "freq": freqs_raw}
     smooth_dict = {"spec": spec_smooth, "freq": freqs_smooth}

src/gmprocess/waveform_processing/snr.py

@@ -4,6 +4,7 @@
 import numpy as np
 from obspy.signal.util import next_pow_2
 
+from gmprocess.metrics.reduction.duration import Duration
 from gmprocess.waveform_processing.fft import compute_and_smooth_spectrum
 from gmprocess.waveform_processing.spectrum import brune_f0, moment_from_magnitude
 from gmprocess.waveform_processing.processing_step import ProcessingStep
@@ -170,17 +171,28 @@ def compute_snr_trace(tr, bandwidth=20.0):
         tr.setCached(
             "smooth_signal_spectrum",
             {
-                "spec": (
-                    tr.getCached("smooth_signal_spectrum")["spec"]
-                    - tr.getCached("smooth_noise_spectrum")["spec"]
-                ),
+                "spec": (tr.getCached("smooth_signal_spectrum")["spec"]),
                 "freq": tr.getCached("smooth_signal_spectrum")["freq"],
             },
         )
 
         smooth_signal_spectrum = tr.getCached("smooth_signal_spectrum")["spec"]
         smooth_noise_spectrum = tr.getCached("smooth_noise_spectrum")["spec"]
-        snr = smooth_signal_spectrum / smooth_noise_spectrum
+
+        # Need to normalize for duration (D). Assumptions:
+        #  - Noise amplitudes scale as sqrt(D)
+        #  - Ground motions at these frequencies are approximately white noise, and thus
+        #    also scale as sqrt(D)
+        #  - Noise is stationary and so D is the full length of the window
+        #  - Ground motion is not stationary, so we'll use the 5-95% significant
+        #    duration.
+        dur_noise = noise.stats.endtime - noise.stats.starttime
+        d595 = Duration([signal], interval=[5.0, 95.0]).result
+        dur_signal = d595[signal.stats.channel]
+        smooth_signal_spectrum = smooth_signal_spectrum / np.sqrt(dur_noise)
+        smooth_signal_spectrum = smooth_signal_spectrum / np.sqrt(dur_signal)
+
+        snr = (smooth_signal_spectrum - smooth_noise_spectrum) / smooth_noise_spectrum
 
         snr_dict = {"snr": snr, "freq": tr.getCached("smooth_signal_spectrum")["freq"]}
         tr.setCached("snr", snr_dict)

tests/data/asdf/asdf_layout.txt

@@ -177,12 +177,8 @@ AuxiliaryData/TraceProcessingParameters/NZ.WTMC/NZ.WTMC.--.HN1_us1000778i_ptest
 AuxiliaryData/TraceProcessingParameters/NZ.WTMC/NZ.WTMC.--.HN2_us1000778i_ptest
 AuxiliaryData/TraceProcessingParameters/NZ.WTMC/NZ.WTMC.--.HNZ_us1000778i_ptest
 AuxiliaryData/WaveFormMetrics
-AuxiliaryData/WaveFormMetrics/NZ.HSES
-AuxiliaryData/WaveFormMetrics/NZ.HSES/NZ.HSES.--.HN_us1000778i_ptest
 AuxiliaryData/WaveFormMetrics/NZ.THZ
 AuxiliaryData/WaveFormMetrics/NZ.THZ/NZ.THZ.--.HN_us1000778i_ptest
-AuxiliaryData/WaveFormMetrics/NZ.WTMC
-AuxiliaryData/WaveFormMetrics/NZ.WTMC/NZ.WTMC.--.HN_us1000778i_ptest
 Provenance
 Provenance/NZ.HSES.--.HN1_us1000778i_ptest
 Provenance/NZ.HSES.--.HN1_us1000778i_unprocessed
@@ -206,11 +202,8 @@ QuakeML
 Waveforms
 Waveforms/NZ.HSES
 Waveforms/NZ.HSES/NZ.HSES.--.HN1__2016-11-13T11:02:20__2016-11-13T11:07:47__us1000778i_unprocessed
-Waveforms/NZ.HSES/NZ.HSES.--.HN1__2016-11-13T11:03:00__2016-11-13T11:05:24__us1000778i_ptest
 Waveforms/NZ.HSES/NZ.HSES.--.HN2__2016-11-13T11:02:20__2016-11-13T11:07:47__us1000778i_unprocessed
-Waveforms/NZ.HSES/NZ.HSES.--.HN2__2016-11-13T11:03:00__2016-11-13T11:05:24__us1000778i_ptest
 Waveforms/NZ.HSES/NZ.HSES.--.HNZ__2016-11-13T11:02:20__2016-11-13T11:07:47__us1000778i_unprocessed
-Waveforms/NZ.HSES/NZ.HSES.--.HNZ__2016-11-13T11:03:00__2016-11-13T11:05:24__us1000778i_ptest
 Waveforms/NZ.HSES/StationXML
 Waveforms/NZ.THZ
 Waveforms/NZ.THZ/NZ.THZ.--.HN1__2016-11-13T11:02:33__2016-11-13T11:07:28__us1000778i_unprocessed
@@ -222,9 +215,6 @@ Waveforms/NZ.THZ/NZ.THZ.--.HNZ__2016-11-13T11:03:13__2016-11-13T11:06:06__us1000
 Waveforms/NZ.THZ/StationXML
 Waveforms/NZ.WTMC
 Waveforms/NZ.WTMC/NZ.WTMC.--.HN1__2016-11-13T11:02:19__2016-11-13T11:07:46__us1000778i_unprocessed
-Waveforms/NZ.WTMC/NZ.WTMC.--.HN1__2016-11-13T11:02:57__2016-11-13T11:05:10__us1000778i_ptest
 Waveforms/NZ.WTMC/NZ.WTMC.--.HN2__2016-11-13T11:02:19__2016-11-13T11:07:46__us1000778i_unprocessed
-Waveforms/NZ.WTMC/NZ.WTMC.--.HN2__2016-11-13T11:02:57__2016-11-13T11:05:10__us1000778i_ptest
 Waveforms/NZ.WTMC/NZ.WTMC.--.HNZ__2016-11-13T11:02:19__2016-11-13T11:07:46__us1000778i_unprocessed
-Waveforms/NZ.WTMC/NZ.WTMC.--.HNZ__2016-11-13T11:02:57__2016-11-13T11:05:10__us1000778i_ptest
 Waveforms/NZ.WTMC/StationXML

tests/gmprocess/waveform_processing/corner_frequencies_test.py

@@ -71,14 +71,16 @@ def test_corner_frequencies():
             cfdict = stream[0].getParameter("corner_frequencies")
             lp.append(cfdict["lowpass"])
             hp.append(cfdict["highpass"])
-    np.testing.assert_allclose(np.sort(hp), [0.024313, 0.024378, 0.025132], atol=1e-5)
+    np.testing.assert_allclose(np.sort(hp), [0.02437835, 0.06223441], atol=1e-5)
 
-    st = processed_streams.select(station="HSES")[0]
+    st = processed_streams.select(station="THZ")[0]
     lps = [tr.getParameter("corner_frequencies")["lowpass"] for tr in st]
     hps = [tr.getParameter("corner_frequencies")["highpass"] for tr in st]
-    np.testing.assert_allclose(np.sort(lps), [100.0, 100.0, 100.0], atol=1e-6)
     np.testing.assert_allclose(
-        np.sort(hps), [0.02431315, 0.02431315, 0.02431315], atol=1e-5
+        np.sort(lps), [21.76376408, 21.76376408, 42.04482076], atol=1e-6
+    )
+    np.testing.assert_allclose(
+        np.sort(hps), [0.02437835, 0.02437835, 0.0345267], atol=1e-5
     )
 
     lp = []
@@ -92,16 +94,17 @@ def test_corner_frequencies():
             lp.append(cfdict["lowpass"])
             hp.append(cfdict["highpass"])
 
-    np.testing.assert_allclose(
-        np.sort(hp), [0.02431315, 0.02437835, 0.02513194], atol=1e-5
-    )
+    np.testing.assert_allclose(np.sort(hp), [0.02437835, 0.06223441], atol=1e-5)
 
     st = processed_streams.select(station="HSES")[0]
     lps = [tr.getParameter("corner_frequencies")["lowpass"] for tr in st]
     hps = [tr.getParameter("corner_frequencies")["highpass"] for tr in st]
-    np.testing.assert_allclose(np.sort(lps), [100.0, 100.0, 100.0], atol=1e-6)
+
+    np.testing.assert_allclose(
+        np.sort(lps), [35.35533906, 36.6021424, 36.6021424], atol=1e-6
+    )
     np.testing.assert_allclose(
-        np.sort(hps), [0.02431315, 0.02431315, 0.02431315], atol=1e-5
+        np.sort(hps), [0.04882812, 0.06223441, 0.06223441], atol=1e-5
     )
 
 

tests/gmprocess/waveform_processing/processing_test.py

@@ -51,7 +51,7 @@ def test_process_streams():
     )
 
     np.testing.assert_allclose(
-        trace_maxes, np.array([157.812449, 240.379521, 263.601519]), rtol=1e-5
+        trace_maxes, np.array([158.989615, 239.48121, 258.440584]), rtol=1e-5
     )