Add functions for broadband TS spectrum - #1664
LOCEANlloydizard wants to merge 32 commits into
Conversation
- Add public `compute_Sv_f()` and `compute_TS_f()` for the API - Add `_cal_complex_samples_f()` prototype - Expose `frequency_resolution` and `range_step` parameters for future spectral calibration - Return prototype `Sv_f` xarray output structure: `(channel, ping_time, svf_range, frequency)`
Use lazy method getattr in calibration dispatcher to avoid access to unsupported compute_Sv_f / compute_TS_f methods on non-FM calibrators
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1664 +/- ##
==========================================
- Coverage 85.58% 85.56% -0.02%
==========================================
Files 79 78 -1
Lines 6998 7074 +76
==========================================
+ Hits 5989 6053 +64
- Misses 1009 1021 +12
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
Hey @LOCEANlloydizard, I thought of a few other params after our discussion the other day so thought I'd lay them out and see what you think. Feel free to ping for more discussions! def compute_Sv_spectrum(
self,
NFFT: int,
range_step: float = 0.5,
window: str = "hann", # TODO: should this be default to None?
frequency_resolution: Optional[float] = None,
) -> xr.Dataset:
"""Compute broadband EK80 Sv spectrum, or Sv(f).
Parameters
----------
NFFT : int
Number of points in the FFT window for computing spectrum.
range_step : float, default 0.5
Range spacing in meters for the output Sv(f) window centres.
window : str, default "hann"
Type of window function to apply to each FFT segment.
Options include "hann", "hamming", etc.
frequency_resolution : optional[float], default None
Frequency spacing in Hz for the output spectral grid.
When this is provided, the output Sv(f) will be interpolated to this frequency grid.
If None, the output Sv(f) will be on the original FFT frequency grid.
Returns
-------
xr.Dataset
A Dataset containing frequency-dependent volume backscattering strength.
"""
def compute_TS_spectrum(
self,
point_locations: xr.DataSet,
NFFT: int,
split_front: float = 0.25,
window: str = "hann",
frequency_resolution: Optional[float] = None,
) -> xr.Dataset:
"""Compute broadband EK80 Sv spectrum, or Sv(f).
Parameters
----------
point_locations : xr.DataSet
Locations of the points for which to compute the spectrum.
This is required to be an xarray.DataSet with dimensions ("target_id") and
data variables "ping_time" and "depth" for each ``target_id``.
NFFT : int
Number of points in the FFT window for computing spectrum.
split_front : float, default 0.25
Each echo spectrum is computed from a segment of the complex echo signal.
This parameter specifies how to identify the segment around the specified point location.
For example, if split_front=0.25, then 25% the NFFT length will be before the point location,
and the remaining 75% will be after the point location.
window : str, default "hann"
Type of window function to apply to each FFT segment.
Options include "hann", "hamming", etc.
frequency_resolution : optional[float], default None
Frequency spacing in Hz for the output spectral grid.
When this is provided, the output Sv(f) will be interpolated to this frequency grid.
If None, the output Sv(f) will be on the original FFT frequency grid.
Returns
-------
xr.Dataset
A Dataset containing frequency-dependent point backscattering strength.
""" |
LOCEANlloydizard
changed the title
- Add compute_TS_spectrum() for EK80 FM complex data - Add broadband TS spectrum processing pipeline and calibration helpers - Add CRIMAC reference dataset and validation tests for TS(f) - Add frequency-dependent absorption validation against CRIMAC - Add broadband Sp and Sv(f) integration tests - Add pooch support for ts_spectrum_example_data bundle - Refactor TS/Sp calibration workflow and API exports - Add cached test-data handling to avoid unnecessary downloads
try to fix error occuring only in ubuntu job
LOCEANlloydizard
changed the title
change the test to a NotImplementedError test
- Remove tqdm import and use - Deprecate BB waveform_mode alias (keep backward compatibility) - Change test_fm_equals_bb test to mark "BB" it as deprecated (others would need it too)
| stacklevel=2, | ||
| ) | ||
|
|
||
| waveform_mode = "BB" if waveform_mode in ("FM", "BB") else waveform_mode |
There was a problem hiding this comment.
This basically makes BB the default instead of FM. I think it'll be better to just wire everything to be based on FM right now, so that when we think it's time to remove support for the BB syntax, we just need to remove this line and the warning.
|
|
||
| if compute_method is None: | ||
| raise ValueError( | ||
| f"{cal_type} calibration is not supported for " f"{echodata.sonar_model} data." |
There was a problem hiding this comment.
| f"{cal_type} calibration is not supported for " f"{echodata.sonar_model} data." | |
| f"{cal_type} computation is not supported for " f"{echodata.sonar_model} data." |
| "Sp": "compute_Sp", | ||
| "TS": "compute_TS", | ||
| "Sv": "compute_Sv", | ||
| # add Sp_spectrum?? |
There was a problem hiding this comment.
Yes I think there should be compute_Sp_spectrum and compute_Sv_spectrum.
updates: from our chat sounds like they will get added after the single target detection PR gets updated (since it's tied to this one in some specific aspects)?
| if "filter_time" in cal_ds_iteration: | ||
| cal_ds_iteration = cal_ds_iteration.drop_vars("filter_time") |
There was a problem hiding this comment.
What are the cases that cal_ds_iteration would not have filter_time? I have forgotten since I review this section last time, but is it from dataset where there is only 1 filter_time?
| if "range_sample" in ds: | ||
| ds["range_sample"].attrs = {"long_name": "Along-range sample number, base 0"} | ||
|
|
||
| if "echo_range" in ds: | ||
| ds["echo_range"].attrs = { | ||
| "long_name": "Range distance", | ||
| "units": "m", | ||
| } | ||
|
|
||
| if "frequency" in ds: | ||
| ds["frequency"].attrs = { | ||
| "long_name": "Frequency", | ||
| "units": "Hz", | ||
| } | ||
|
|
||
| ds[cal_type].attrs = { | ||
| "long_name": { | ||
| "Sv": "Volume backscattering strength (Sv re 1 m-1)", | ||
| "Sp": "Point scattering strength (Sp re 1 m^2)", | ||
| "TS": "Target strength (TS re 1 m^2)", | ||
| "Sv": "Volume backscattering strength (Sv re 1 m-1)", | ||
| "TS_spectrum": "Frequency-dependent target strength spectrum (TS(f) re 1 m^2)", | ||
| }[cal_type], | ||
| "units": "dB", | ||
| } |
There was a problem hiding this comment.
I think it'll be nice for us to have a consolidated place to specify all the attributes for variables. Right now we define these in an ad-hoc way and missed many of them. It seems would be a good idea to just have a YAML or similar to list the corresponding long_name, units, and other attributes for specific variables, and in functions that generates organized datasets, always pass the dataset through a function that calls the YAML to add attributes -- so like a decorator? I'll add an issue for this.
There was a problem hiding this comment.
one more take on this:
since _add_attrs is only used once and these vars are "under our control," wouldn't it be better to expect that ALL of these vars would exist so should error out when they don't exist (and therefore there shouldn't be if-else statement)? Or, if some of these only exist for some cal_type, I think it would be better to have these assignments to go as expected behaviors instead of if-else.
| For CW data, Sp is computed from received power samples on the range grid. | ||
| For EK80 broadband/FM complex data, Sp is computed after pulse compression | ||
| and represents a band-averaged point-scattering-strength echogram. |
There was a problem hiding this comment.
I am not exactly sure what you mean between "on the range grid" vs "point-scattering-strength echogram". Aren't both still on the range grid, either it's the original one or the one after pulse compression?
| try: | ||
| method_name = compute_methods[cal_type] | ||
| except KeyError: | ||
| raise ValueError(f"Unsupported calibration type: {cal_type}") from None | ||
|
|
||
| compute_method = getattr(cal_obj, method_name, None) | ||
|
|
||
| if compute_method is None: | ||
| raise ValueError( | ||
| f"{cal_type} calibration is not supported for " f"{echodata.sonar_model} data." |
There was a problem hiding this comment.
I am not following the logic here - wouldn't it be better to just do a check if cal_type exists in compute_methods instead of have it error out? And also from None suppresses the trace back (I just learned this!) - is this intentional?
And then there's a second getattr check that seems to have a very similar intention.
Maybe these can be consolidated to something simpler, that checks if cal_type is not None nor valid, and if otherwise spit out clear error messages?
| def _get_transducer_halves( | ||
| pc: xr.DataArray, | ||
| ) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray]: | ||
| """Calculate half-transducer pulse-compressed signals. | ||
|
|
||
| Equivalent to CRIMAC ``calcTransducerHalves`` for 4-sector transducers. | ||
| """ | ||
| if pc.sizes["beam"] != 4: | ||
| raise NotImplementedError( | ||
| "Transducer halves are only defined for 4-sector split-beam data." | ||
| ) | ||
|
|
||
| pc_fore = 0.5 * (pc.isel(beam=2) + pc.isel(beam=3)) | ||
| pc_aft = 0.5 * (pc.isel(beam=0) + pc.isel(beam=1)) | ||
| pc_star = 0.5 * (pc.isel(beam=0) + pc.isel(beam=3)) | ||
| pc_port = 0.5 * (pc.isel(beam=1) + pc.isel(beam=2)) | ||
|
|
||
| return pc_fore, pc_aft, pc_star, pc_port | ||
|
|
||
|
|
||
| def _get_splitbeam_angles( | ||
| pc: xr.DataArray, | ||
| gamma_alongship, | ||
| gamma_athwartship, | ||
| ) -> tuple[xr.DataArray, xr.DataArray]: | ||
| """Calculate raw split-beam physical angles before angle-offset correction. | ||
|
|
||
| For 4-sector data this follows CRIMAC ``calcAngles``. For 3-sector data, | ||
| the sector geometry follows the same convention used by ``add_splitbeam_angle``. | ||
| Angle offsets are not applied here because TS(f) beam compensation applies | ||
| frequency-dependent offsets later. | ||
| """ | ||
| if pc.sizes["beam"] == 4: | ||
| pc_fore, pc_aft, pc_star, pc_port = _get_transducer_halves(pc) | ||
|
|
||
| y_theta = pc_fore * np.conj(pc_aft) | ||
| y_phi = pc_star * np.conj(pc_port) | ||
|
|
||
| theta = np.rad2deg( | ||
| np.arcsin(np.arctan2(np.imag(y_theta), np.real(y_theta)) / gamma_alongship) | ||
| ) | ||
|
|
||
| phi = np.rad2deg(np.arcsin(np.arctan2(np.imag(y_phi), np.real(y_phi)) / gamma_athwartship)) | ||
| else: | ||
| raise NotImplementedError( | ||
| f"Split-beam angle calculation is not implemented for {pc.sizes['beam']} sectors." | ||
| ) | ||
|
|
||
| theta.name = "angle_alongship" | ||
| phi.name = "angle_athwartship" | ||
|
|
||
| return theta, phi |
There was a problem hiding this comment.
Maybe these could be combined and refactored with what was implemented in add_splitbeam_angle so that we have a single place to maintain for these calculations? Also because the code there handles 3-sector transducers as well from a recent PR.
There was a problem hiding this comment.
Just to clarify: I saw your notes under _get_beam_compensated_gain and the differences. What I meant is that the underlying functions can be merged, so that there is only 1 copy of the operation for split beam angle calculation from complex signals. The additional compensation related to angle offset can be a flag and hence optional for the add_splitbeam_angle function, for example.
| def _get_average_signal( | ||
| signal: xr.DataArray, | ||
| ) -> xr.DataArray: | ||
| """Average complex signal over transducer sectors. | ||
|
|
||
| Equivalent to CRIMAC ``calcAverageSignal``. | ||
| """ | ||
| out = signal.mean(dim="beam") | ||
| out.name = "average_signal" | ||
|
|
||
| return out |
There was a problem hiding this comment.
do we need to keep this as a function? Seems too small and easier to keep track of if the .mean operation is used directly instead.
| def _get_pulse_compressed_signal( | ||
| beam: xr.Dataset, | ||
| matched_filter: Dict, | ||
| ) -> xr.DataArray: | ||
| """Calculate pulse-compressed complex samples for each transducer sector. | ||
|
|
||
| Equivalent to CRIMAC ``calcPulseCompressedSignals``. | ||
| """ | ||
| pc = compress_pulse( | ||
| backscatter=beam["backscatter_r"] + 1j * beam["backscatter_i"], | ||
| chirp=matched_filter, | ||
| ) | ||
| pc = pc / get_norm_fac(chirp=matched_filter) | ||
| pc.name = "pulse_compressed_signal" | ||
|
|
||
| return pc |
There was a problem hiding this comment.
The same operation is used in split_beam_angle.get_angle_complex_samples() - maybe use this function there too. This would be part of the merging/refactoring to have a single place to host split-beam angle calculation code as commented here.
| gain_db = np.interp( | ||
| frequency, | ||
| self.vend["cal_frequency"].values, | ||
| self.vend["gain"].sel(cal_channel_id=channel).values, | ||
| ) | ||
|
|
||
| angle_offset_alongship = np.interp( | ||
| frequency, | ||
| self.vend["cal_frequency"].values, | ||
| self.vend["angle_offset_alongship"].sel(cal_channel_id=channel).values, | ||
| ) | ||
|
|
||
| angle_offset_athwartship = np.interp( | ||
| frequency, | ||
| self.vend["cal_frequency"].values, | ||
| self.vend["angle_offset_athwartship"].sel(cal_channel_id=channel).values, | ||
| ) | ||
|
|
||
| beamwidth_alongship = np.interp( | ||
| frequency, | ||
| self.vend["cal_frequency"].values, | ||
| self.vend["beamwidth_alongship"].sel(cal_channel_id=channel).values, | ||
| ) | ||
|
|
||
| beamwidth_athwartship = np.interp( | ||
| frequency, | ||
| self.vend["cal_frequency"].values, | ||
| self.vend["beamwidth_athwartship"].sel(cal_channel_id=channel).values, | ||
| ) | ||
|
|
||
| beam_correction_db = self._get_beam_correction( | ||
| theta=theta, | ||
| phi=phi, | ||
| angle_offset_alongship=angle_offset_alongship, | ||
| angle_offset_athwartship=angle_offset_athwartship, | ||
| beamwidth_alongship=beamwidth_alongship, | ||
| beamwidth_athwartship=beamwidth_athwartship, | ||
| ) |
There was a problem hiding this comment.
This part of the code that does the scaling of parameters based on frequency overlaps with what's in cal_params.get_cal_params_EK(). Similar to the split-beam angle calculation, maybe it makes sense for us to come up with a module that does the scaling, and just use it in both here and cal parameter preparation?
| return mf_auto[idx_start:idx_stop] | ||
|
|
||
|
|
||
| def _extract_target_from_range_gate( |
There was a problem hiding this comment.
I feel this function should probably go to a separate module and not in ek80_complex.py since it's not related to the complex data operations. Not exactly sure how to call it.
There was a problem hiding this comment.
Looking at the ops in this function, I am not sure I agree with this approach to change everything to pure numpy array and then use pure indexing to find the values corresponding with a specific range gate. I think the same operations can be completed in xarray and the intention would be clearer because the dimensions/coordinates are explicit.
| except ValueError: | ||
| continue |
There was a problem hiding this comment.
This is very weird - silent failure without any warning???
| if mf_auto_red.size < n_fft: | ||
| mf_pad = np.zeros(n_fft, dtype=complex) | ||
| mf_pad[: mf_auto_red.size] = mf_auto_red | ||
| mf_auto_red = mf_pad | ||
| elif mf_auto_red.size > n_fft: | ||
| mf_auto_red = mf_auto_red[:n_fft] |
There was a problem hiding this comment.
This padding or truncation operation is very important. I think we should have it in the docstring and also the example notebook, to demonstrate what would happen with the different length specified when computing the spectrum.
| def _compute_ts_spectrum_power( | ||
| normalized_spectrum: np.ndarray, | ||
| n_beams: int, | ||
| z_et: float, | ||
| z_er: float, | ||
| ): | ||
| """Convert normalised TS(f) spectrum to received power spectrum. | ||
|
|
||
| Equivalent to CRIMAC ``calcPowerFreqTS``. | ||
| """ | ||
| return _compute_complex_power( | ||
| normalized_spectrum=normalized_spectrum, | ||
| n_beams=n_beams, | ||
| z_et=z_et, | ||
| z_er=z_er, | ||
| ) |
There was a problem hiding this comment.
I think this is over factoring, since _compute_complex_power is only used once here. _compute_complex_power also overlaps with what's in _compute_power_from_complex_signal. I think these should all be combined and refactored to use the same function.
| ts = _compute_ts_spectrum_calibrated( | ||
| power_spectrum=power_spectrum, | ||
| target_range=target_range, | ||
| frequency=frequency, | ||
| sound_speed=sound_speed, | ||
| absorption_f=absorption_f, | ||
| transmit_power=transmit_power, | ||
| gain_f=gain_f, | ||
| ) |
There was a problem hiding this comment.
I'd suggest keeping the "gut" exposed here for readability. The docstring for the function can just become inline comments.
| def compute_TS_spectrum(echodata: EchoData, **kwargs) -> xr.Dataset: | ||
| """ | ||
| Compute broadband frequency-dependent target strength spectrum, TS(f), | ||
| from EK80 broadband/FM complex data. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| point_locations : xr.Dataset | ||
| Locations of targets for which TS(f) should be computed. | ||
| Must contain ``channel``, ``ping_time``, and ``target_range`` for each | ||
| ``target_id``. If ``target_range_min`` and ``target_range_max`` are | ||
| provided, they define the target echo segment. Otherwise, the segment | ||
| is built around ``target_range`` using ``NFFT`` and ``split_front``. | ||
|
|
||
| NFFT : int, optional | ||
| Number of FFT points used to compute the target spectrum. If not | ||
| provided, a value is inferred from the output frequency grid. | ||
|
|
||
| n_f_points : int, optional | ||
| Number of frequency points in the output TS(f) spectrum. Used when | ||
| ``frequency_resolution`` is not provided. | ||
|
|
||
| split_front : float, default 0.25 | ||
| Each echo spectrum is computed from a segment of the complex echo signal. | ||
| This parameter specifies how to position that segment around the target location | ||
| when only ``target_range`` is provided. For example, if ``split_front=0.25``, | ||
| then 25% of the NFFT window is placed before ``target_range`` and the remaining | ||
| 75% after it. | ||
|
|
||
| window : str, tuple, float or None, default None | ||
| Window passed directly to ``scipy.signal.get_window``. If ``None``, | ||
| a rectangular/boxcar window is used. | ||
|
|
||
| frequency_resolution : float, optional | ||
| Desired spacing of the output frequency grid in Hz. Used to define the | ||
| frequency grid on which TS(f) is evaluated. |
There was a problem hiding this comment.
I think allowing n_f_points and frequency_resolution is confusing when they are together with NFFT and split_front -- these two are messy enough. I would suggest that we show people how to do interpolation outside by operating on the output of this function, instead of doing it silently inside the function. Only the selected signal, NFFT, and the related padding really matters, the interpolation is cosmetic. I worry that by giving all these options people could get lost in what is really important for processing the data.
I also suggest split_front to be remained as fraction_before_point_location since "split" could be confused with split-beam.
| if not out_by_channel: | ||
| raise ValueError("No valid TS_spectrum targets produced.") | ||
|
|
||
| return xr.concat(out_by_channel, dim="channel") |
There was a problem hiding this comment.
What happens to the frequency coordinate (specific to the channel) when you do xr.concat? Does it concat the multiple frequency vectors sequentially?
| if not ts_list: | ||
| continue |
There was a problem hiding this comment.
Better to not have silent skip like this since the user may end up having an empty dataset and it can be hard to keep track of why.
|
|
||
| raise ValueError(f"Unsupported calibration type: {cal_type}") | ||
|
|
||
| def _compute_cal(self, cal_type, **kwargs) -> xr.Dataset: |
There was a problem hiding this comment.
I think this was in place before this PR already, but looking back at this I think we can actually do without the _compute_cal structure and just dispatch the following 3 functions directly in compute_Sv/Sp/TS/TS_spectrum:
_cal_power_samples_cal_complex_samples_cal_complex_samples_f
Is there an advantage in keeping this _compute_cal function?
|
|
||
| return xr.concat(out_by_channel, dim="channel") | ||
|
|
||
| def _cal_complex_samples_f( |
There was a problem hiding this comment.
Seems like you can just absorb what's in _cal_complex_samples_TS_spectrum here? I don't see much advantage of separating these 2. But I will make some comments about what I think could be factor out in _cal_complex_samples_TS_spectrum.
| return theta, phi | ||
|
|
||
|
|
||
| def _compute_power_from_complex_signal( |
There was a problem hiding this comment.
This overlaps with _compute_complex_power
There was a problem hiding this comment.
Hey @LOCEANlloydizard : Thanks for the PR! Sorry it took a very long time for me to review this!
Below are a few higher-level things that may make the inline comments easier to read:
- There are some redundancy within the new code and with the old code that I think would benefit from combining and refactoring, such as:
- frequency-based scaling of calibration parameters
- split-beam angle estimation between beampattern compensation and
add_splitbeam_angle - computing power from complex signals with the right scaling related to impedance
- I'd suggest removing
n_f_pointsandfrequency_resolutionas input arguments when computing spectrum, so that the function output would just be what people would get based onNFFTortarget_range_min/max; any downstream interpolation the users will do outside of the function, so that there's no potential hidden/unintended errors. - refactoring and commenting code in
_cal_complex_samples_TS_spectrum- it was hard to follow the logic because the code jumps between doing different things -- I think those can be better bunched up to be functions called by this method, to make things easier to follow
- there is very minimal inline comments in this function -- it would be better to add the intention for the operations and any important caveats as comments
- would it make sense to just always first compute the split-beam angles before estimating the TS? That is the logical step and we can enforce it in the code, so that the
_cal_complex_samples_TS_spectrumfunction for example does not need to jump between computing the spectrum and the angles. It'll then just be computing the spectrum and grab the corresponding angles to get the right beampattern compensation.
Maybe we can talk through the function when we meet to make a plan for the above changes? I think it would be great to make this PR a one-time thing to get clean code, since EK80 is stable now and the code likely won't change much later, so it'll be a worthwhile investment.
This branch has not been deployed
3 participants
Implement a CRIMAC-equivalent broadband TS(f) processing pipeline, including pulse compression, matched-filter autocorrelation, normalized spectra, frequency-dependent calibration, and beam compensation (Test broadband processing against CRIMAC outputs #1033)
Add public compute_TS_spectrum() and compute_Sp() APIs, and extend the calibration framework to support frequency-dependent TS(f) products (Add broadband Sv and TS code #1034)
Update compute_TS() for EK80 complex data to return a beam-compensated band-averaged TS product, while keeping compute_TS_spectrum() as the frequency-dependent TS(f) implementation (Add broadband Sv and TS code #1034, discussion from Add variable Sv attributes to output dataset from compute_Sv #582)
Return TS(f) as an xarray Dataset indexed by (channel, target_id, frequency) (Add broadband Sv and TS code #1034)
Deprecate the public "BB" waveform mode alias in favor of "FM" while preserving backward compatibility (Change
waveform_modedefault vocabulary to "FM" for broadband EK80 data #1651)Refs: #582, #1033, #1034, #1651