diff --git a/CHANGELOG.md b/CHANGELOG.md index 32b935142..93cb6accb 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/src/gmprocess/metrics/reduction/duration.py b/src/gmprocess/metrics/reduction/duration.py index 7829670f9..7651b82fe 100644 --- a/src/gmprocess/metrics/reduction/duration.py +++ b/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. diff --git a/src/gmprocess/waveform_processing/fft.py b/src/gmprocess/waveform_processing/fft.py index 2a5fabe56..1374e8ba1 100644 --- a/src/gmprocess/waveform_processing/fft.py +++ b/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} diff --git a/src/gmprocess/waveform_processing/snr.py b/src/gmprocess/waveform_processing/snr.py index c5718982c..5126de37e 100644 --- a/src/gmprocess/waveform_processing/snr.py +++ b/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) diff --git a/tests/data/asdf/asdf_layout.txt b/tests/data/asdf/asdf_layout.txt index e9e43213f..20fe28c27 100644 --- a/tests/data/asdf/asdf_layout.txt +++ b/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 \ No newline at end of file diff --git a/tests/gmprocess/waveform_processing/corner_frequencies_test.py b/tests/gmprocess/waveform_processing/corner_frequencies_test.py index eb2942817..aab750f32 100755 --- a/tests/gmprocess/waveform_processing/corner_frequencies_test.py +++ b/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 ) diff --git a/tests/gmprocess/waveform_processing/processing_test.py b/tests/gmprocess/waveform_processing/processing_test.py index 35565b68e..050c8ad5e 100755 --- a/tests/gmprocess/waveform_processing/processing_test.py +++ b/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 )