diff --git a/tests/test_bb_filters.py b/tests/test_bb_filters.py new file mode 100644 index 00000000..b89d27fa --- /dev/null +++ b/tests/test_bb_filters.py @@ -0,0 +1,233 @@ +"""Tests for the broadband source-filter correction in `bb_sim`. + +The correction exists because both solvers low-pass their *source time +functions* and `bb_sim` then low-passes the LF leg again, which leaves a hole in +the transition band. These tests pin the three things that can go wrong +silently: the analytic model of `qcore`'s filter drifting from the real one, the +correction not actually flattening the recombination, and the solver filters +being read wrongly out of a realisation. +""" + +import json +from pathlib import Path + +import numpy as np +import pytest + +from qcore import timeseries +from workflow.realisations import RealisationMetadata +from workflow.scripts.bb_sim import ( + _HIGHPASS_SHIFT, + _LOWPASS_SHIFT, + BB_FILTER_ORDER, + BB_FILTER_PASSES, + CorrectionLeg, + Solver, + SourceLowpass, + butterworth_gain, + recombination_gains, + solver_source_lowpass, +) + +DT = 0.005 +FLO = 1.0 +# The two configurations actually in use: SW4's v26.7.1Hz `prefilter` command, +# and EMOD3D's `bfilt=4` with `flo = min_vs / (5 * resolution) = 1.0 Hz`. +SOURCES = [ + SourceLowpass(order=2, passes=2, corner=1.0), + SourceLowpass(order=4, passes=2, corner=1.0), +] +# The band the correction is aimed at. Outside it the LF leg carries nothing +# (above) or the filters are flat (below), so neither says anything useful. +BAND = (0.2, 3.0) + + +def measured_response(band: timeseries.Band, frequencies: np.ndarray) -> np.ndarray: + """The magnitude response of `bwfilter` itself, from an impulse. + + The impulse sits in the middle of the record, not at the start: `bwfilter` + passes ``padtype=None`` to `sosfiltfilt`, so an impulse at index 0 has its + acausal half truncated and the measured response is meaningless. + """ + nt = 8192 + impulse = np.zeros(nt) + impulse[nt // 2] = 1.0 + response = timeseries.bwfilter(impulse, DT, FLO, band) + spectrum = np.abs(np.fft.rfft(response)) + return np.interp(frequencies, np.fft.rfftfreq(nt, DT), spectrum) + + +@pytest.mark.parametrize( + ("band", "shift", "btype"), + [ + (timeseries.Band.LOWPASS, _LOWPASS_SHIFT, "lowpass"), + (timeseries.Band.HIGHPASS, _HIGHPASS_SHIFT, "highpass"), + ], +) +def test_analytic_gain_matches_qcore( + band: timeseries.Band, shift: float, btype: str +) -> None: + """`butterworth_gain` must reproduce what `bwfilter` actually does. + + The shifts are mirrored from qcore rather than imported, so this is the + guard against the two definitions drifting apart. + """ + frequencies = np.linspace(*BAND, 200) + analytic = butterworth_gain( + frequencies, BB_FILTER_ORDER, FLO * shift, BB_FILTER_PASSES, DT, band=btype + ) + np.testing.assert_allclose( + analytic, measured_response(band, frequencies), atol=2e-3 + ) + + +def test_matched_pair_is_power_complementary() -> None: + """Without a source filter the legs already sum to one in power.""" + frequencies = np.linspace(*BAND, 400) + low = butterworth_gain( + frequencies, BB_FILTER_ORDER, FLO * _LOWPASS_SHIFT, BB_FILTER_PASSES, DT + ) + high = butterworth_gain( + frequencies, + BB_FILTER_ORDER, + FLO * _HIGHPASS_SHIFT, + BB_FILTER_PASSES, + DT, + band="highpass", + ) + # The pair is only approximately complementary by design; 0.05 in ln is the + # deviation the shifts leave behind, and it is what the correction restores. + assert np.abs(np.log(np.sqrt(low**2 + high**2))).max() < 0.05 + + +def test_no_source_filter_is_a_no_op() -> None: + frequencies = np.linspace(*BAND, 50) + lf_gain, hf_gain = recombination_gains(frequencies, DT, FLO, None, CorrectionLeg.LF) + np.testing.assert_array_equal(lf_gain, 1.0) + np.testing.assert_array_equal(hf_gain, 1.0) + + +@pytest.mark.parametrize("source", SOURCES, ids=["sw4", "emod3d"]) +def test_uncorrected_recombination_has_a_hole(source: SourceLowpass) -> None: + """The defect the correction exists for: without it there is a real deficit.""" + frequencies = np.linspace(*BAND, 400) + low = butterworth_gain( + frequencies, BB_FILTER_ORDER, FLO * _LOWPASS_SHIFT, BB_FILTER_PASSES, DT + ) + high = butterworth_gain( + frequencies, + BB_FILTER_ORDER, + FLO * _HIGHPASS_SHIFT, + BB_FILTER_PASSES, + DT, + band="highpass", + ) + source_gain = butterworth_gain( + frequencies, source.order, source.corner, source.passes, DT + ) + deficit = np.log(np.sqrt((source_gain * low) ** 2 + high**2)) + assert deficit.min() < -0.1, "expected a transition-band hole to correct" + + +def power_sum_error( + source: SourceLowpass, leg: CorrectionLeg, frequencies: np.ndarray +) -> np.ndarray: + """|ln| of the corrected power sum against the matched pair's own.""" + low = butterworth_gain( + frequencies, BB_FILTER_ORDER, FLO * _LOWPASS_SHIFT, BB_FILTER_PASSES, DT + ) + high = butterworth_gain( + frequencies, + BB_FILTER_ORDER, + FLO * _HIGHPASS_SHIFT, + BB_FILTER_PASSES, + DT, + band="highpass", + ) + source_gain = butterworth_gain( + frequencies, source.order, source.corner, source.passes, DT + ) + lf_gain, hf_gain = recombination_gains(frequencies, DT, FLO, source, leg) + corrected = (lf_gain * source_gain * low) ** 2 + (hf_gain * high) ** 2 + return np.abs(np.log(np.sqrt(corrected / (low**2 + high**2)))) + + +@pytest.mark.parametrize("source", SOURCES, ids=["sw4", "emod3d"]) +def test_both_restores_the_power_sum_exactly(source: SourceLowpass) -> None: + """Scaling the two legs together is exact everywhere, by construction.""" + frequencies = np.linspace(*BAND, 400) + assert power_sum_error(source, CorrectionLeg.BOTH, frequencies).max() < 1e-9 + + +@pytest.mark.parametrize("source", SOURCES, ids=["sw4", "emod3d"]) +def test_lf_restores_the_power_sum_where_the_lf_leg_matters( + source: SourceLowpass, +) -> None: + """Exact below the matching frequency; only clipped where the LF is spent.""" + below = np.linspace(BAND[0], FLO, 200) + assert power_sum_error(source, CorrectionLeg.LF, below).max() < 1e-9 + # Above it the boost is clipped, but the LF leg carries so little there that + # the total barely moves. This bound is the promise the docstring makes. + assert ( + power_sum_error(source, CorrectionLeg.LF, np.linspace(*BAND, 400)).max() < 0.01 + ) + + +@pytest.mark.parametrize("source", SOURCES, ids=["sw4", "emod3d"]) +def test_hf_cannot_restore_below_the_matching_frequency( + source: SourceLowpass, +) -> None: + """A known limitation, pinned so it cannot become a silent surprise. + + Filling an LF deficit from the HF leg below `flo` needs a boost of ten or + more on stochastic content with no valid long-period part, so the correction + clips and the power sum stays short. + """ + below = np.linspace(BAND[0], FLO, 200) + _, hf_gain = recombination_gains(below, DT, FLO, source, CorrectionLeg.HF) + assert hf_gain.max() > 5.0, "expected the HF correction to be ill-conditioned" + + +@pytest.mark.parametrize("source", SOURCES, ids=["sw4", "emod3d"]) +def test_both_is_the_best_conditioned_leg(source: SourceLowpass) -> None: + """`both` is the option that never needs a large boost.""" + frequencies = np.linspace(*BAND, 400) + lf_gain, hf_gain = recombination_gains( + frequencies, DT, FLO, source, CorrectionLeg.BOTH + ) + assert max(lf_gain.max(), hf_gain.max()) < 1.5 + + +def test_solver_filters_are_read_from_the_realisation(tmp_path: Path) -> None: + """The declared parameters must reach `SourceLowpass` unmangled. + + A realisation carrying only `metadata` and `resolution` falls back to the + defaults for everything else, which is what pins the shipped values: SW4's + `prefilter order=2 passes=2 fc2=1.0`, and EMOD3D's `bfilt=4` with a corner + derived the way `create_e3d_par` derives it. The EMOD3D corner coming out at + 1.0 Hz is the same number `e3d.par` carries as `flo`. + """ + realisation_ffp = tmp_path / "realisation.json" + realisation_ffp.write_text( + json.dumps( + { + "metadata": { + "name": "test", + "version": "1", + "defaults_version": "26.7.1Hz", + "tag": "gcmt", + }, + "resolution": {"resolution": 0.1}, + } + ) + ) + version = RealisationMetadata.read_from_realisation( + realisation_ffp + ).defaults_version + + assert solver_source_lowpass(Solver.SW4, realisation_ffp, version) == SourceLowpass( + order=2, passes=2, corner=1.0 + ) + assert solver_source_lowpass( + Solver.EMOD3D, realisation_ffp, version + ) == SourceLowpass(order=4, passes=2, corner=1.0) diff --git a/workflow/scripts/bb_sim.py b/workflow/scripts/bb_sim.py index 037c3955..b729d642 100644 --- a/workflow/scripts/bb_sim.py +++ b/workflow/scripts/bb_sim.py @@ -35,6 +35,8 @@ See the output of `bb-sim --help`. """ +import dataclasses +from enum import StrEnum from pathlib import Path from typing import Annotated @@ -48,7 +50,15 @@ from qcore import cli, timeseries from site_calculation import amplification from workflow import log_utils, realisations -from workflow.realisations import BroadbandParameters, RealisationMetadata +from workflow.realisations import ( + BroadbandParameters, + EMOD3DParameters, + RealisationMetadata, + Resolution, + SW4Parameters, + VelocityModelParameters, + find_command, +) from workflow.schemas import SiteAmpModel # Site amplification model -> (amplification function, model frequencies). @@ -69,6 +79,322 @@ G = 1 / 981.0 TARGET_CHUNK_BYTES = 256 * 2**20 +# `qcore.timeseries.bwfilter` is an order-4 Butterworth applied by `sosfiltfilt`, +# i.e. twice. These are the corner shifts it applies internally so that after +# both passes the response is exactly 1/sqrt(2) at the frequency it was handed. +# Mirrored here rather than imported because they are private to qcore; +# `tests/test_bb_filters.py` measures the real filter and fails if they drift. +BB_FILTER_ORDER = 4 +BB_FILTER_PASSES = 2 +_HIGHPASS_SHIFT = (np.sqrt(2) - 1) ** (1 / 8) +_LOWPASS_SHIFT = 1 / _HIGHPASS_SHIFT + +# EMOD3D's `tfilter` takes a `phase` argument and runs a reverse pass when it is +# zero (source.c:1989). Every call site hard-codes `int phase = 0`, so the +# EMOD3D source low-pass is zero-phase, not causal. +EMOD3D_SOURCE_PASSES = 2 + +# Ceiling on the boost any correction may apply. The two real configurations +# need at most ~2.4x; anything approaching this bound means the source filter +# rolls off faster than the target and the leg is being reconstructed from +# content that is not there. +MAX_BOOST = 10.0 +# Below this the high-frequency leg carries too little of the signal for the +# HF-side correction to be meaningful, and dividing by it amplifies nothing but +# numerical noise. +HF_GAIN_FLOOR = 1e-2 + + +class CorrectionLeg(StrEnum): + """Which leg of the recombination absorbs the source-filter correction.""" + + LF = "lf" + """Restore the low-frequency leg by dividing out its source filter. + + The physically direct choice: it puts the LF leg back to what the solver + would have produced without the extra filter, and touches nothing else. + Exact below the matching frequency. Above it the boost is clipped, which + costs under 0.01 in ln of the total because the LF leg carries almost + nothing there. + """ + HF = "hf" + """Leave the LF alone and fill the missing power from the HF leg. + + Only meaningful near and above the matching frequency. Below it, filling an + LF deficit from the HF leg needs a boost of ten or more applied to + stochastic content that has no valid long-period part, so the correction is + clipped and the power sum is *not* restored. `warn_if_ill_conditioned` says + so when it happens. + """ + BOTH = "both" + """Scale both legs together so the power sum is restored exactly. + + The best-conditioned option -- it needs a boost of under 1.4 anywhere -- but + it makes up part of the deficit from the high-frequency leg, which is a + different physical claim from restoring the low-frequency one. + """ + + +class Solver(StrEnum): + """The low-frequency solver that produced the LF waveforms.""" + + EMOD3D = "emod3d" + SW4 = "sw4" + + +@dataclasses.dataclass(frozen=True) +class SourceLowpass: + """A low-pass a solver has already applied to its source time functions. + + Both solvers filter the source rather than the output. The wave equation is + linear, so that is equivalent to filtering every trace, which is what lets + `bb_sim` divide the filter back out here instead of re-running the solver. + """ + + order: int + """Butterworth order.""" + passes: int + """1 for a causal filter, 2 for the zero-phase forward-and-back pair.""" + corner: float + """Corner frequency in Hz.""" + + +def solver_source_lowpass( + solver: Solver, realisation_ffp: Path, defaults_version: str +) -> SourceLowpass | None: + """The low-pass `solver` applied to its sources, read from the realisation. + + Parameters + ---------- + solver : Solver + The solver that produced the low-frequency waveforms. + realisation_ffp : Path + Path to the realisation file. + defaults_version : str + The realisation's defaults version. + + Returns + ------- + SourceLowpass | None + The filter, or None if the solver applied none. + + Raises + ------ + ValueError + If the solver applied a filter this function cannot describe, rather + than silently correcting for the wrong thing. + """ + if solver is Solver.SW4: + sw4_config = SW4Parameters.read_from_realisation_or_defaults( + realisation_ffp, defaults_version + ) + prefilter = find_command(sw4_config.commands, "prefilter") + if prefilter is None: + return None + parameters = prefilter.parameters + if parameters.get("type") != "lowpass": + raise ValueError( + f"SW4 prefilter is type={parameters.get('type')!r}; only 'lowpass' " + "can be corrected for here" + ) + return SourceLowpass( + order=int(parameters["order"]), + passes=int(parameters["passes"]), + corner=float(parameters["fc2"]), + ) + + emod3d_config = EMOD3DParameters.read_from_realisation_or_defaults( + realisation_ffp, defaults_version + ) + if not emod3d_config.bfilt: + # `tfilter` is guarded by `if(bfilt)`, so zero means no filter at all. + return None + if emod3d_config.fhi: + raise ValueError( + f"EMOD3D fhi={emod3d_config.fhi} applies a source high-pass as well as " + "the low-pass; correcting for the low-pass alone would be wrong" + ) + # EMOD3D's flo is not stored in the realisation or the LF file. It is + # derived in `create_e3d_par.create_duration_parameters`, and is repeated + # here from the same two inputs. + velocity_model = VelocityModelParameters.read_from_realisation_or_defaults( + realisation_ffp, defaults_version + ) + resolution = Resolution.read_from_realisation_or_defaults( + realisation_ffp, defaults_version + ) + return SourceLowpass( + order=emod3d_config.bfilt, + passes=EMOD3D_SOURCE_PASSES, + corner=velocity_model.min_vs / (5 * resolution.resolution), + ) + + +def butterworth_gain( + frequencies: np.ndarray, + order: int, + corner: float, + passes: int, + dt: float, + band: str = "lowpass", +) -> np.ndarray: + """Magnitude response of a Butterworth filter applied `passes` times. + + Parameters + ---------- + frequencies : np.ndarray + Frequencies (Hz) to evaluate the response at. + order : int + Butterworth order. + corner : float + Corner frequency in Hz. + passes : int + Number of times the filter is applied. Two passes square the magnitude. + dt : float + Sample interval in seconds. + band : str + Either 'lowpass' or 'highpass'. + + Returns + ------- + np.ndarray + The magnitude response at `frequencies`. + """ + sos = sp.signal.butter(order, corner, btype=band, output="sos", fs=1.0 / dt) + _, response = sp.signal.sosfreqz(sos, worN=2 * np.pi * frequencies * dt) + return np.abs(response) ** passes + + +def warn_if_ill_conditioned( + dt: float, flo: float, source: SourceLowpass | None, leg: CorrectionLeg +) -> None: + """Report where the correction is clipped and cannot restore the power sum. + + Called once, rather than per chunk, so the message is not repeated for every + block dask schedules. + + Parameters + ---------- + dt : float + Broadband sample interval. + flo : float + The LF/HF matching frequency. + source : SourceLowpass | None + The solver's source filter. + leg : CorrectionLeg + Which leg absorbs the correction. + """ + if source is None: + return + frequencies = np.logspace(np.log10(flo / 10), np.log10(flo * 5), 400) + lf_gain, hf_gain = recombination_gains(frequencies, dt, flo, source, leg) + clipped = np.isclose(lf_gain, MAX_BOOST) | np.isclose(hf_gain, MAX_BOOST) + if not clipped.any(): + return + logger = log_utils.get_logger(__name__) + logger.warning( + "The source-filter correction is clipped, so the power sum is not fully " + "restored across the whole band." + + ( + " Below the matching frequency the high-frequency leg has no valid " + "content to supply; consider --filter lf or --filter both." + if leg is CorrectionLeg.HF + else " This is above the matching frequency, where the " + "low-frequency leg carries almost nothing, so the effect on the " + "total is small." + ), + correction_leg=str(leg), + max_boost=MAX_BOOST, + clipped_from_hz=float(frequencies[clipped].min()), + clipped_to_hz=float(frequencies[clipped].max()), + ) + + +def recombination_gains( + frequencies: np.ndarray, + dt: float, + flo: float, + source: SourceLowpass | None, + leg: CorrectionLeg, +) -> tuple[np.ndarray, np.ndarray]: + """Per-leg correction for the source filter already baked into the LF. + + The LF and HF legs are independent realisations, so they add in *power*, and + `bwfilter`'s shifted pair is designed to be power-complementary: + ``|H_lp|^2 + |H_hp|^2 = 1``. A source low-pass multiplies the LF leg a second + time, which breaks that identity and leaves a hole in the transition band. + These gains put the power sum back where the pair intended. + + Parameters + ---------- + frequencies : np.ndarray + Frequencies (Hz) the gains are sampled at. + dt : float + Broadband sample interval. + flo : float + The LF/HF matching frequency. + source : SourceLowpass | None + The solver's source filter. None applies no correction. + leg : CorrectionLeg + Which leg absorbs the correction. + + Returns + ------- + tuple[np.ndarray, np.ndarray] + Multiplicative gains for the LF and HF legs, both ones if `source` is None. + """ + ones = np.ones_like(frequencies) + if source is None: + return ones, ones + + lowpass = butterworth_gain( + frequencies, BB_FILTER_ORDER, flo * _LOWPASS_SHIFT, BB_FILTER_PASSES, dt + ) + highpass = butterworth_gain( + frequencies, + BB_FILTER_ORDER, + flo * _HIGHPASS_SHIFT, + BB_FILTER_PASSES, + dt, + band="highpass", + ) + source_gain = butterworth_gain( + frequencies, source.order, source.corner, source.passes, dt + ) + + # What the matched pair is supposed to deliver, and what it delivers once + # the source filter has been applied to the LF leg a second time. + target_power = lowpass**2 + highpass**2 + current_power = (source_gain * lowpass) ** 2 + highpass**2 + + match leg: + case CorrectionLeg.LF: + # Dividing the source filter back out restores the LF leg exactly. + # The denominator is floored rather than the quotient clipped: the + # two are identical for a non-negative gain, but this does not go + # via an intermediate inf where the source filter has rolled off + # to zero. + return 1.0 / np.maximum(source_gain, 1.0 / MAX_BOOST), ones + case CorrectionLeg.HF: + # Make the shortfall up from the HF side, leaving the LF as it is. + shortfall = np.sqrt( + np.clip(target_power - (source_gain * lowpass) ** 2, 0.0, None) + ) + hf_gain = np.where( + highpass > HF_GAIN_FLOOR, + np.clip( + shortfall / np.maximum(highpass, HF_GAIN_FLOOR), None, MAX_BOOST + ), + 1.0, + ) + return ones, hf_gain + case CorrectionLeg.BOTH: + # One factor on both legs restores the total without either leg + # carrying the whole boost. Well conditioned everywhere, because + # `current_power` is bounded below by the high-pass leg. + shared = np.clip(np.sqrt(target_power / current_power), None, MAX_BOOST) + return shared, shared + def align_datasets( lf: xr.Dataset, hf: xr.Dataset, dt: float @@ -235,6 +561,8 @@ def _process_bb_chunk( fhightop: float, fmax: float, site_amp_model: SiteAmpModel, + source_lowpass: SourceLowpass | None = None, + correction_leg: CorrectionLeg = CorrectionLeg.LF, ) -> xr.Dataset: """Compute broadband waveforms for a chunk of stations. @@ -242,6 +570,10 @@ def _process_bb_chunk( waveforms, then merges them with the low-frequency waveforms using a matched pair of high-pass and low-pass Butterworth filters. + Where `source_lowpass` is given, the solver's own source filter is divided + back out of the recombination first, so the LF leg is filtered once rather + than twice and the matched pair's power sum is restored. + Parameters ---------- dset : xr.Dataset @@ -267,6 +599,11 @@ def _process_bb_chunk( (highpass end of the amplification band). site_amp_model : SiteAmpModel The site amplification model to apply. + source_lowpass : SourceLowpass | None + The low-pass the solver already applied to its source time functions. + None leaves the recombination exactly as it was. + correction_leg : CorrectionLeg + Which leg absorbs the correction for `source_lowpass`. Returns ------- @@ -285,6 +622,16 @@ def _process_bb_chunk( # sampled at. n_fft = pyfftw.next_fast_len(nt) fft_freqs = np.fft.rfftfreq(n_fft, dt) + lf_gain, hf_gain = recombination_gains( + fft_freqs, dt, flo, source_lowpass, correction_leg + ) + correcting = source_lowpass is not None + if correcting: + # `amplify_waveform` takes one gain curve per station, and these are the + # same curve for every station in the chunk. + stations = lf_waveform.shape[1] + lf_gain = np.tile(lf_gain, (stations, 1)) + hf_gain = np.tile(hf_gain, (stations, 1)) # The amplification models require float64 inputs. vs30 = dset["vs30"].values.astype(np.float64) @@ -313,6 +660,11 @@ def _process_bb_chunk( lf_filtered = timeseries.bwfilter( lf_waveform[i], dt, flo, timeseries.Band.LOWPASS ) + # Applied after the matched pair rather than in place of it, so that + # with no correction the arithmetic is bit-for-bit what it always was. + if correcting: + lf_filtered = amplification.amplify_waveform(lf_filtered, lf_gain, n_fft) + hf_filtered = amplification.amplify_waveform(hf_filtered, hf_gain, n_fft) bb_waveform[i] = (hf_filtered + lf_filtered) * G return dset.drop_vars(["lf_waveform", "hf_waveform", "vs30"]).assign( @@ -332,6 +684,12 @@ def combine_hf_and_lf( Path, typer.Argument(exists=True, dir_okay=False) ], output_ffp: Annotated[Path, typer.Argument(dir_okay=False, writable=True)], + solver: Annotated[ + Solver | None, typer.Option("--solver", case_sensitive=False) + ] = None, + filter_leg: Annotated[ + CorrectionLeg, typer.Option("--filter", case_sensitive=False) + ] = CorrectionLeg.LF, ) -> None: """Combine low-frequency and high-frequency seismic waveforms. @@ -347,11 +705,34 @@ def combine_hf_and_lf( File containing high-frequency waveform data. output_ffp : Path Path to the output file where the combined broadband waveforms will be saved. + solver : Solver | None + The solver that produced the low-frequency waveforms. Given, the source + low-pass it already applied (SW4's `prefilter` command, or EMOD3D's + `bfilt` and `flo`) is divided back out of the recombination, so the LF + leg is filtered once rather than twice. Omitted, the recombination is + unchanged. + filter_leg : CorrectionLeg + Which leg absorbs that correction: `lf` restores the low-frequency leg, + `hf` fills the missing power from the high-frequency leg instead, and + `both` scales the two together so neither is boosted hard. Ignored + without `--solver`. """ metadata = RealisationMetadata.read_from_realisation(realisation_ffp) broadband_config = BroadbandParameters.read_from_realisation_or_defaults( realisation_ffp, metadata.defaults_version ) + source_lowpass = ( + None + if solver is None + else solver_source_lowpass(solver, realisation_ffp, metadata.defaults_version) + ) + if solver is not None and source_lowpass is None: + log_utils.get_logger(__name__).warning( + "The solver applied no source low-pass in this realisation, so the " + "recombination is unchanged.", + solver=str(solver), + ) + # Open lazily (no dask) and select the common stations *before* chunking. # The LF and HF files store stations in different orders, so selecting after # chunking is an all-to-all dask shuffle in which every output chunk depends @@ -375,6 +756,7 @@ def combine_hf_and_lf( hf = hf.sel(station=common_stations).chunk(chunking) bb_dt = min(lf.attrs["dt"], hf.attrs["dt"]) + warn_if_ill_conditioned(bb_dt, broadband_config.flo, source_lowpass, filter_leg) if not np.isclose(lf.attrs["dt"], bb_dt): lf = resample_signal(lf, bb_dt) @@ -437,6 +819,8 @@ def combine_hf_and_lf( fhightop=broadband_config.fhightop, fmax=broadband_config.fmax, site_amp_model=broadband_config.site_amp_version, + source_lowpass=source_lowpass, + correction_leg=filter_leg, ), template=template, ) @@ -450,6 +834,17 @@ def combine_hf_and_lf( fmax=broadband_config.fmax, site_amp_model=str(broadband_config.site_amp_version), ) + # Recorded so a broadband file says whether it was corrected and for what. + # Without this the correction is invisible downstream, and two files that + # differ by it look identical. + if source_lowpass is not None: + attributes |= { + "source_filter_solver": str(solver), + "source_filter_order": source_lowpass.order, + "source_filter_passes": source_lowpass.passes, + "source_filter_corner": source_lowpass.corner, + "source_filter_correction": str(filter_leg), + } # Attributes, unlike station coordinates, are *not* carried through # map_blocks: `template` above only has `combined`'s. The LF file's # supergrid width describes the run that produced the waveforms, and