From 83c1310b1c9d5e4cdfd0e5cdcb8d4581f38ce0b0 Mon Sep 17 00:00:00 2001 From: Peter Hollender Date: Tue, 31 Mar 2026 10:46:45 -0400 Subject: [PATCH 1/9] add frequency-dependent sensitivity Addresses #439 --- src/openlifu/plan/solution.py | 9 +- src/openlifu/sim/kwave_if.py | 2 +- src/openlifu/xdc/element.py | 97 +++++++++--------- src/openlifu/xdc/transducer.py | 148 +++++++++++++++------------- src/openlifu/xdc/transducerarray.py | 4 +- tests/test_transducer.py | 91 ++++++++++++++++- 6 files changed, 235 insertions(+), 116 deletions(-) diff --git a/src/openlifu/plan/solution.py b/src/openlifu/plan/solution.py index ffe49b5c..a206378e 100644 --- a/src/openlifu/plan/solution.py +++ b/src/openlifu/plan/solution.py @@ -256,7 +256,14 @@ def analyze(self, apodization = self.apodizations[focus_index] origin = self.transducer.get_effective_origin(apodizations=apodization, units=options.distance_units) - output_signal_Pa = self.transducer.calc_output(input_signal_V, dt, delays=self.delays[focus_index, :], apod=self.apodizations[focus_index, :]) + output_signal_Pa = self.transducer.calc_output( + input_signal_V, + cycles=self.pulse.duration * self.pulse.frequency, + frequency=self.pulse.frequency, + dt=dt, + delays=self.delays[focus_index, :], + apod=self.apodizations[focus_index, :], + ) p0_Pa = np.max(output_signal_Pa, axis=1) mainlobe_mask = get_mask( diff --git a/src/openlifu/sim/kwave_if.py b/src/openlifu/sim/kwave_if.py index 3c54a6ec..6687ea03 100644 --- a/src/openlifu/sim/kwave_if.py +++ b/src/openlifu/sim/kwave_if.py @@ -156,7 +156,7 @@ def run_simulation(arr: xdc.Transducer, if _source is not None: source = _source else: - source_mat = arr.calc_output(input_signal, kgrid.dt, delays, apod) + source_mat = arr.calc_output(input_signal, cycles=cycles, frequency=freq, dt=kgrid.dt, delays=delays, apod=apod) if arr.crosstalk_frac != 0: # Simulate crosstalk by adding additional elements to the array for each element that # is within the crosstalk distance, with the signal scaled by the crosstalk fraction. diff --git a/src/openlifu/xdc/element.py b/src/openlifu/xdc/element.py index 1fea76d2..b961a3ae 100644 --- a/src/openlifu/xdc/element.py +++ b/src/openlifu/xdc/element.py @@ -10,6 +10,33 @@ from openlifu.util.units import getunitconversion +def sensitivity_at_frequency(sensitivity: float | dict[float, float], frequency: float) -> float: + if isinstance(sensitivity, dict): + freqs = np.array(list(sensitivity.keys()), dtype=np.float64) + values = np.array(list(sensitivity.values()), dtype=np.float64) + return float(np.interp(frequency, freqs, values, left=values[0], right=values[-1])) + return float(sensitivity) + + +def generate_drive_signal(input_signal, cycles: float, frequency: float, dt: float) -> np.ndarray: + """Generate a drive signal with duration constrained by cycles/frequency.""" + if dt <= 0: + raise ValueError("dt must be positive.") + if frequency <= 0: + raise ValueError("frequency must be positive.") + if cycles <= 0: + raise ValueError("cycles must be positive.") + n_samples = max(1, int(np.round(cycles / (frequency * dt)))) + if np.isscalar(input_signal): + t = np.arange(n_samples, dtype=np.float64) * dt + return float(input_signal) * np.sin(2 * np.pi * frequency * t) + base = np.asarray(input_signal, dtype=np.float64).reshape(-1) + drive_signal = np.zeros(n_samples, dtype=np.float64) + n_copy = min(n_samples, len(base)) + drive_signal[:n_copy] = base[:n_copy] + return drive_signal + + def matrix2xyz(matrix): x = matrix[0, 3] y = matrix[1, 3] @@ -43,15 +70,9 @@ class Element: size: Annotated[np.ndarray, OpenLIFUFieldData("Size", "Size of the element in 2D")] = field(default_factory=lambda: np.array([1., 1.])) """ Size of the element in 2D as a numpy array [width, length].""" - sensitivity: Annotated[float | None, OpenLIFUFieldData("Sensitivity", "Sensitivity of the element (Pa/V)")] = None + sensitivity: Annotated[float | dict[float, float], OpenLIFUFieldData("Sensitivity", "Sensitivity of the element (Pa/V), scalar or {frequency_hz: sensitivity}")] = 1.0 """Sensitivity of the element (Pa/V)""" - impulse_response: Annotated[np.ndarray | None, OpenLIFUFieldData("Impulse response", "Impulse response of the element")] = None - """Impulse response of the element, can be a single value or an array of values. If an array, `impulse_dt` must be set to the time step of the impulse response. Is convolved with the input signal.""" - - impulse_dt: Annotated[float | None, OpenLIFUFieldData("Impulse response timestep", """Impulse response timestep""")] = None - """Impulse response timestep. If `impulse_response` is an array, this is the time step of the impulse response.""" - pin: Annotated[int, OpenLIFUFieldData("Pin", "Channel pin to which the element is connected")] = -1 """Channel pin to which the element is connected. 1-(64*number of modules).""" @@ -68,14 +89,18 @@ def __post_init__(self): self.size = np.array(self.size, dtype=np.float64) if self.size.shape != (2,): raise ValueError("Size must be a 2-element array.") - if self.impulse_response is not None: - if isinstance(self.impulse_response, int | float): - self.impulse_response = np.array([self.impulse_response]) - self.impulse_response = np.array(self.impulse_response, dtype=np.float64) - if self.impulse_response.ndim != 1: - raise ValueError("Impulse response must be a 1-dimensional array.") - if len(self.impulse_response)>1 and self.impulse_dt is None: - raise ValueError("Impulse response timestep must be set if impulse response is an array.") + if self.sensitivity is None: + self.sensitivity = 1.0 + elif isinstance(self.sensitivity, dict): + if len(self.sensitivity) == 0: + raise ValueError("Sensitivity dictionary must not be empty.") + mapping = {float(k): float(v) for k, v in self.sensitivity.items()} + freqs = np.array(sorted(mapping.keys()), dtype=np.float64) + if np.any(np.diff(freqs) <= 0): + raise ValueError("Sensitivity dictionary frequencies must be strictly increasing.") + self.sensitivity = {float(f): mapping[float(f)] for f in freqs} + else: + self.sensitivity = float(self.sensitivity) @property def x(self): @@ -141,17 +166,12 @@ def length(self): def length(self, value): self.size[1] = value - def calc_output(self, input_signal, dt): - if self.impulse_response is None: - filtered_signal = input_signal * 1 - elif len(self.impulse_response) == 1: - filtered_signal = input_signal * self.impulse_response[0] - else: - impulse = self.interp_impulse_response(dt) - filtered_signal = np.convolve(input_signal, impulse, mode='full') - if self.sensitivity is not None: - filtered_signal = filtered_signal * self.sensitivity - return filtered_signal + def get_sensitivity(self, frequency: float) -> float: + return sensitivity_at_frequency(self.sensitivity, frequency) + + def calc_output(self, input_signal, cycles: float, frequency: float, dt: float): + drive_signal = generate_drive_signal(input_signal, cycles=cycles, frequency=frequency, dt=dt) + return drive_signal * self.get_sensitivity(frequency) def copy(self): return copy.deepcopy(self) @@ -225,17 +245,6 @@ def get_angle(self, units="rad"): roll = np.degrees(self.roll) return el, az, roll - def interp_impulse_response(self, dt=None): - if dt is None: - dt = self.impulse_dt - n0 = len(self.impulse_response) - t0 = self.impulse_dt * np.arange(n0) - t1 = np.arange(0, t0[-1] + dt, dt) - impulse_response = np.interp(t1, t0, self.impulse_response) - impulse_t = np.arange(len(impulse_response)) * dt - impulse_t = impulse_t - np.mean(impulse_t) - return impulse_response, impulse_t - def distance_to_point(self, point, units=None, matrix=np.eye(4)): units = self.units if units is None else units pos = np.concatenate([self.get_position(units=units), [1]]) @@ -273,10 +282,7 @@ def to_dict(self): "size": self.size.tolist(), "pin": self.pin, "units": self.units} - if self.impulse_response is not None: - d["impulse_response"] = self.impulse_response.tolist() - if self.impulse_dt is not None: - d["impulse_dt"] = self.impulse_dt + d["sensitivity"] = self.sensitivity return d @staticmethod @@ -286,8 +292,9 @@ def from_dict(d): d["position"] = np.array([d.pop('x'), d.pop('y'), d.pop('z')]) d["orientation"] = np.array([d.pop('az'), d.pop('el'), d.pop('roll')]) d["size"] = np.array([d.pop('w'), d.pop('l')]) - if "impulse_response" in d and d["impulse_response"] is not None: - d["impulse_response"] = np.array(d["impulse_response"]) - if "impulse_dt" in d and d["impulse_dt"] is not None: - d["impulse_dt"] = float(d["impulse_dt"]) + # Backward compatibility: legacy impulse keys are ignored. + d.pop("impulse_response", None) + d.pop("impulse_dt", None) + if "sensitivity" not in d or d["sensitivity"] is None: + d["sensitivity"] = 1.0 return Element(**d) diff --git a/src/openlifu/xdc/transducer.py b/src/openlifu/xdc/transducer.py index a07d50c3..ebeaaaca 100644 --- a/src/openlifu/xdc/transducer.py +++ b/src/openlifu/xdc/transducer.py @@ -11,11 +11,34 @@ from openlifu.util.annotations import OpenLIFUFieldData from openlifu.util.units import getunitconversion -from openlifu.xdc.element import Element +from openlifu.xdc.element import ( + Element, + generate_drive_signal, + sensitivity_at_frequency, +) DIMS = ['x', 'y', 'z'] LDIMS = Literal['x','y','z'] + +def _combine_sensitivities( + base_sensitivity: float | dict[float, float], + scale_sensitivity: float | dict[float, float], +) -> float | dict[float, float]: + if isinstance(base_sensitivity, dict) and isinstance(scale_sensitivity, dict): + base_keys = tuple(base_sensitivity.keys()) + scale_keys = tuple(scale_sensitivity.keys()) + if base_keys != scale_keys: + raise ValueError("Cannot combine sensitivity dictionaries with different frequency keys.") + return {frequency: base_sensitivity[frequency] * scale_sensitivity[frequency] for frequency in base_keys} + if isinstance(base_sensitivity, dict): + factor = float(scale_sensitivity) + return {frequency: value * factor for frequency, value in base_sensitivity.items()} + if isinstance(scale_sensitivity, dict): + factor = float(base_sensitivity) + return {frequency: factor * value for frequency, value in scale_sensitivity.items()} + return float(base_sensitivity) * float(scale_sensitivity) + @dataclass class Transducer: id: Annotated[str, OpenLIFUFieldData("Transducer ID", "Unique identifier for transducer")] = "transducer" @@ -55,8 +78,8 @@ class Transducer: The units of this transform are assumed to be the native units of the transducer, the `Transducer.units` field. """ - sensitivity: Annotated[float | None, OpenLIFUFieldData("Sensitivity", "Sensitivity of the element (Pa/V)")] = None - """Sensitivity of the element (Pa/V)""" + sensitivity: Annotated[float | dict[float, float], OpenLIFUFieldData("Sensitivity", "Sensitivity of the element (Pa/V), scalar or {frequency_hz: sensitivity}")] = 1.0 + """Sensitivity of the transducer (Pa/V), scalar or frequency-dependent dictionary.""" crosstalk_frac: Annotated[float, OpenLIFUFieldData("Crosstalk fraction", "Fraction of the signal that leaks into other elements due to crosstalk")] = 0.0 """Fraction of the signal that leaks into other elements due to crosstalk""" @@ -64,12 +87,6 @@ class Transducer: crosstalk_dist: Annotated[float, OpenLIFUFieldData("Crosstalk distance", "Distance within which elements experience crosstalk")] = 0.0 """Distance within which elements experience crosstalk""" - impulse_response: Annotated[np.ndarray | None, OpenLIFUFieldData("Impulse response", "Impulse response of the element")] = None - """Impulse response of the element, can be a single value or an array of values. If an array, `impulse_dt` must be set to the time step of the impulse response. Is convolved with the input signal.""" - - impulse_dt: Annotated[float | None, OpenLIFUFieldData("Impulse response timestep", """Impulse response timestep""")] = None - """Impulse response timestep. If `impulse_response` is an array, this is the time step of the impulse response.""" - module_invert: Annotated[List[bool], OpenLIFUFieldData("Invert polarity", "Whether to invert the polarity of the transducer output, per module")] = field(default_factory=lambda: [False]) """Whether to invert the polarity of the transducer output""" @@ -79,38 +96,36 @@ def __post_init__(self): self.name = self.id for element in self.elements: element.rescale(self.units) - if self.impulse_response is not None: - self.impulse_response = np.array(self.impulse_response, dtype=np.float64) - if self.impulse_response.ndim != 1 or len(self.impulse_response)<2: - raise ValueError("Impulse response must be a 1-dimensional array.") - if self.impulse_dt is None: - raise ValueError("Impulse response timestep must be set if impulse response is set.") - - - def interp_impulse_response(self, dt=None): - if dt is None: - dt = self.impulse_dt - n0 = len(self.impulse_response) - t0 = self.impulse_dt * np.arange(n0) - t1 = np.arange(0, t0[-1] + dt, dt) - impulse_response = np.interp(t1, t0, self.impulse_response) - impulse_t = np.arange(len(impulse_response)) * dt - impulse_t = impulse_t - np.mean(impulse_t) - return impulse_response, impulse_t - - def calc_output(self, input_signal, dt, delays: np.ndarray = None, apod: np.ndarray = None): + if self.sensitivity is None: + self.sensitivity = 1.0 + elif isinstance(self.sensitivity, dict): + if len(self.sensitivity) == 0: + raise ValueError("Sensitivity dictionary must not be empty.") + mapping = {float(k): float(v) for k, v in self.sensitivity.items()} + freqs = np.array(sorted(mapping.keys()), dtype=np.float64) + if np.any(np.diff(freqs) <= 0): + raise ValueError("Sensitivity dictionary frequencies must be strictly increasing.") + self.sensitivity = {float(f): mapping[float(f)] for f in freqs} + else: + self.sensitivity = float(self.sensitivity) + + def get_sensitivity(self, frequency: float) -> float: + return sensitivity_at_frequency(self.sensitivity, frequency) + + def calc_output(self, input_signal, cycles: float, frequency: float, dt: float, delays: np.ndarray = None, apod: np.ndarray = None): if delays is None: delays = np.zeros(self.numelements()) if apod is None: apod = np.ones(self.numelements()) - if self.impulse_response is None: - filtered_input_signal = input_signal * 1 - else: - impulse = self.interp_impulse_response(dt) - filtered_input_signal = np.convolve(input_signal, impulse, mode='full') - if self.sensitivity is not None: - filtered_input_signal = filtered_input_signal * self.sensitivity - outputs = [np.concatenate([np.zeros(int(delay/dt)), a*element.calc_output(filtered_input_signal, dt)],axis=0) for element, delay, a, in zip(self.elements, delays, apod)] + drive_signal = generate_drive_signal(input_signal, cycles=cycles, frequency=frequency, dt=dt) + base_output = drive_signal * self.get_sensitivity(frequency) + outputs = [ + np.concatenate( + [np.zeros(int(delay / dt)), a * element.get_sensitivity(frequency) * base_output], + axis=0, + ) + for element, delay, a, in zip(self.elements, delays, apod) + ] max_len = max([len(o) for o in outputs]) output_signal = np.zeros([self.numelements(), max_len]) for i, o in enumerate(outputs): @@ -235,22 +250,32 @@ def get_standoff_transform_in_units(self, units:str) -> np.ndarray: @staticmethod def merge(list_of_transducers:List[Transducer], offset_pins:bool=False, offset_indices:bool=False, merge_mismatched_sensitivity=True, merged_attrs:dict={}) -> Transducer: array_copies = [arr.copy() for arr in list_of_transducers] - sensitivities = np.array([arr.sensitivity for arr in array_copies if arr.sensitivity is not None]) - if len(sensitivities) > 0 and len(sensitivities) < len(array_copies): - raise ValueError("If one transducer has a sensitivity, all must have a sensitivity.") - if len(set(sensitivities)) > 1: - if not merge_mismatched_sensitivity: - raise ValueError("Transducers have different sensitivities. Use merge_mismatched_sensitivity=True to merge the relative sensitivities into the merged elements") + dict_key_sets = set() + for array in array_copies: + if isinstance(array.sensitivity, dict): + dict_key_sets.add(tuple(array.sensitivity.keys())) + for el in array.elements: + if isinstance(el.sensitivity, dict): + dict_key_sets.add(tuple(el.sensitivity.keys())) + if len(dict_key_sets) > 1: + raise ValueError("Cannot merge sensitivities with different frequency keys.") + + sensitivity_signatures = [] + for array in array_copies: + if isinstance(array.sensitivity, dict): + sensitivity_signatures.append((tuple(array.sensitivity.keys()), tuple(array.sensitivity.values()))) else: - max_sensitivity = sensitivities.max() - relative_sensitivities = sensitivities/max_sensitivity - for array, relative_sensitivity in zip(array_copies, relative_sensitivities): - for el in array.elements: - if el.sensitivity is not None: - el.sensitivity = el.sensitivity * relative_sensitivity - else: - el.sensitivity = relative_sensitivity - array.sensitivity = max_sensitivity + sensitivity_signatures.append(float(array.sensitivity)) + + if not merge_mismatched_sensitivity and len(set(sensitivity_signatures)) > 1: + raise ValueError("Transducers have different sensitivities. Use merge_mismatched_sensitivity=True to merge them into the merged elements") + + for array in array_copies: + transducer_sensitivity = array.sensitivity + for el in array.elements: + el.sensitivity = _combine_sensitivities(el.sensitivity, transducer_sensitivity) + array.sensitivity = 1.0 + merged_array = array_copies[0] for xform_array in array_copies[1:]: if offset_pins: @@ -287,12 +312,6 @@ def sort_by_pin(self): def to_dict(self): d = self.__dict__.copy() d["elements"] = [element.to_dict() for element in d["elements"]] - if self.impulse_response is None: - del d["impulse_response"] - else: - d["impulse_response"] = d["impulse_response"].tolist() - if self.impulse_dt is None: - del d["impulse_dt"] d["standoff_transform"] = d["standoff_transform"].tolist() return d @@ -347,12 +366,11 @@ def from_file(filename): def from_dict(d, **kwargs): d = d.copy() d["elements"] = [Element.from_dict(element) for element in d["elements"]] - if "impulse_response" in d and d["impulse_response"] is not None: - if len(d["impulse_response"]) == 1 and "sensitivity" not in d: - d["sensitivity"] = d["impulse_response"][0] - del d["impulse_response"] - else: - d["impulse_response"] = np.array(d["impulse_response"]) + # Backward compatibility: legacy impulse fields are ignored. + d.pop("impulse_response", None) + d.pop("impulse_dt", None) + if "sensitivity" not in d or d["sensitivity"] is None: + d["sensitivity"] = 1.0 if "standoff_transform" in d and d["standoff_transform"] is not None: d["standoff_transform"] = np.array(d["standoff_transform"]) return Transducer(**d, **kwargs) @@ -385,8 +403,6 @@ def gen_matrix_array(nx=2, ny=2, pitch=1, kerf=0, units="mm", **kwargs): pitch: distance between element centers kerf: distance between element edges units: units of the array dimensions - impulse_response: impulse response of the elements - impulse_dt: time step of the impulse response id: unique identifier name: name of the array attrs: additional attributes diff --git a/src/openlifu/xdc/transducerarray.py b/src/openlifu/xdc/transducerarray.py index 40e61b9c..9eaaf978 100644 --- a/src/openlifu/xdc/transducerarray.py +++ b/src/openlifu/xdc/transducerarray.py @@ -45,8 +45,8 @@ def from_dict(data: dict): if "attrs" in d: if "standoff_transform" in d["attrs"] and d["attrs"]["standoff_transform"] is not None: d["attrs"]["standoff_transform"] = np.array(d["attrs"]["standoff_transform"]) - if "impulse_response" in d["attrs"] and d["attrs"]["impulse_response"] is not None: - d["attrs"]["impulse_response"] = np.array(d["attrs"]["impulse_response"]) + d["attrs"].pop("impulse_response", None) + d["attrs"].pop("impulse_dt", None) return TransducerArray(**d) diff --git a/tests/test_transducer.py b/tests/test_transducer.py index 526dfc58..cc79872f 100644 --- a/tests/test_transducer.py +++ b/tests/test_transducer.py @@ -105,6 +105,95 @@ def test_transducer_array_to_transducer_data_types(transducer_array_id): transducer_array : TransducerArray = load_transducer_array(transducer_array_id) transducer = transducer_array.to_transducer() assert isinstance(transducer.standoff_transform, np.ndarray) - assert isinstance(transducer.impulse_response, np.ndarray) + assert not hasattr(transducer, "impulse_response") + assert not hasattr(transducer, "impulse_dt") if len(transducer.elements) > 0: assert isinstance(transducer.elements[0], Element) + + +def test_transducer_calc_output_interpolates_dictionary_sensitivity(): + transducer = Transducer.gen_matrix_array( + nx=1, + ny=1, + units="mm", + sensitivity={100e3: 1.0, 300e3: 3.0}, + ) + transducer.elements[0].sensitivity = 1.0 + input_signal = np.array([1.0, -1.0, 0.5], dtype=float) + + output_mid = transducer.calc_output(input_signal, cycles=3, frequency=200e3, dt=1e-7) + output_low = transducer.calc_output(input_signal, cycles=3, frequency=100e3, dt=1e-7) + + np.testing.assert_allclose(output_mid[0, :len(input_signal)], 2.0 * input_signal) + np.testing.assert_allclose(output_low[0, :len(input_signal)], 1.0 * input_signal) + + +def test_element_calc_output_generates_signal_from_scalar_input(): + element = Element(sensitivity=2.0) + cycles = 4 + frequency = 100e3 + dt = 1e-7 + n_samples = int(np.round(cycles / (frequency * dt))) + + output = element.calc_output(3.0, cycles=cycles, frequency=frequency, dt=dt) + t = np.arange(n_samples) * dt + expected = 2.0 * 3.0 * np.sin(2 * np.pi * frequency * t) + + np.testing.assert_allclose(output, expected) + + +def test_element_calc_output_enforces_cycles_duration_for_array_input(): + element = Element(sensitivity=1.0) + cycles = 1 + frequency = 200e3 + dt = 1e-6 + n_samples = int(np.round(cycles / (frequency * dt))) + input_signal = np.arange(20, dtype=float) + + output = element.calc_output(input_signal, cycles=cycles, frequency=frequency, dt=dt) + + assert len(output) == n_samples + np.testing.assert_allclose(output, input_signal[:n_samples]) + + +def test_merge_pushes_transducer_sensitivity_into_elements(): + transducer_a = Transducer.gen_matrix_array( + nx=1, + ny=1, + units="mm", + sensitivity={100e3: 2.0, 300e3: 4.0}, + ) + transducer_b = Transducer.gen_matrix_array( + nx=1, + ny=1, + units="mm", + sensitivity={100e3: 3.0, 300e3: 6.0}, + ) + transducer_a.elements[0].sensitivity = 5.0 + transducer_b.elements[0].sensitivity = 7.0 + + merged = Transducer.merge([transducer_a, transducer_b], merge_mismatched_sensitivity=True) + + assert merged.sensitivity == 1.0 + assert merged.elements[0].sensitivity == {100e3: 10.0, 300e3: 20.0} + assert merged.elements[1].sensitivity == {100e3: 21.0, 300e3: 42.0} + + +def test_merge_rejects_mismatched_sensitivity_keys(): + transducer_a = Transducer.gen_matrix_array( + nx=1, + ny=1, + units="mm", + sensitivity={100e3: 2.0, 300e3: 4.0}, + ) + transducer_b = Transducer.gen_matrix_array( + nx=1, + ny=1, + units="mm", + sensitivity={100e3: 3.0, 400e3: 6.0}, + ) + transducer_a.elements[0].sensitivity = {100e3: 5.0, 300e3: 7.0} + transducer_b.elements[0].sensitivity = {100e3: 11.0, 400e3: 13.0} + + with pytest.raises(ValueError, match="different frequency keys"): + Transducer.merge([transducer_a, transducer_b], merge_mismatched_sensitivity=True) From 1f24481fca8e30be62ef639b5efd370bcb69a032 Mon Sep 17 00:00:00 2001 From: Peter Hollender Date: Tue, 31 Mar 2026 11:18:25 -0400 Subject: [PATCH 2/9] add loading functions to xdc.util To aid in using the functionality in #439 --- src/openlifu/xdc/util.py | 274 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 274 insertions(+) diff --git a/src/openlifu/xdc/util.py b/src/openlifu/xdc/util.py index 9c2b8d38..c9970c15 100644 --- a/src/openlifu/xdc/util.py +++ b/src/openlifu/xdc/util.py @@ -1,11 +1,243 @@ from __future__ import annotations import json +import re + +import pandas as pd +import xarray as xa from openlifu.util.types import PathLike from openlifu.xdc.transducer import Transducer from openlifu.xdc.transducerarray import TransducerArray +FOCAL_GAIN_LUT = xa.DataArray.from_dict( + {'dims': ('f0', 'crosstalk'), + 'attrs': {}, + 'data': [[2.807589054107666, + 3.2286391258239746, + 3.649686813354492, + 3.8181092739105225, + 4.07073974609375, + 4.491786956787109, + 4.912837505340576, + 5.333885669708252, + 5.754938125610352, + 6.175983428955078, + 6.597033500671387, + 7.01808500289917], + [2.90433931350708, + 3.3324313163757324, + 3.760524272918701, + 3.931760549545288, + 4.188616752624512, + 4.616710662841797, + 5.044803142547607, + 5.472893714904785, + 5.900986671447754, + 6.3290791511535645, + 6.757172107696533, + 7.185482501983643], + [2.9909276962280273, + 3.428293466567993, + 3.865659713745117, + 4.040605068206787, + 4.3030242919921875, + 4.740390777587891, + 5.1777544021606445, + 5.615119934082031, + 6.052487373352051, + 6.489851951599121, + 6.927217483520508, + 7.364583969116211], + [3.0771772861480713, + 3.5201354026794434, + 3.9643561840057373, + 4.142045497894287, + 4.408576965332031, + 4.852799415588379, + 5.297021865844727, + 5.741242408752441, + 6.185462474822998, + 6.629685878753662, + 7.073910713195801, + 7.518129348754883], + [3.170368194580078, + 3.617199182510376, + 4.064029693603516, + 4.242762565612793, + 4.51104211807251, + 4.961826801300049, + 5.4126152992248535, + 5.863399505615234, + 6.314184665679932, + 6.764969825744629, + 7.215755462646484, + 7.666543960571289], + [3.242729663848877, + 3.6979565620422363, + 4.153182506561279, + 4.335274696350098, + 4.6084089279174805, + 5.064931869506836, + 5.521984100341797, + 5.979033946990967, + 6.4360833168029785, + 6.89313268661499, + 7.350184440612793, + 7.8072357177734375], + [3.327850103378296, + 3.783803701400757, + 4.245124340057373, + 4.429731845855713, + 4.706640720367432, + 5.168155193328857, + 5.6296706199646, + 6.091186046600342, + 6.552703380584717, + 7.014214992523193, + 7.475728988647461, + 7.937244415283203], + [3.415055990219116, + 3.8783867359161377, + 4.341715335845947, + 4.527048110961914, + 4.8050456047058105, + 5.268375396728516, + 5.731706142425537, + 6.195037841796875, + 6.658364295959473, + 7.12183952331543, + 7.587955474853516, + 8.054073333740234], + [5.705652713775635, + 5.966911792755127, + 6.2332234382629395, + 6.343565940856934, + 6.509076118469238, + 6.784928798675537, + 7.060781478881836, + 7.336633682250977, + 7.612486362457275, + 7.888339519500732, + 8.164192199707031, + 8.447998046875], + [5.73893404006958, + 5.998416423797607, + 6.260423183441162, + 6.365228652954102, + 6.522432327270508, + 6.784440994262695, + 7.046449184417725, + 7.310236930847168, + 7.583878517150879, + 7.8575215339660645, + 8.131163597106934, + 8.404805183410645], + [5.780664920806885, + 6.028132915496826, + 6.275601387023926, + 6.3745880126953125, + 6.523068904876709, + 6.777853965759277, + 7.03900146484375, + 7.300150394439697, + 7.561298370361328, + 7.822445869445801, + 8.083595275878906, + 8.350235939025879], + [5.814091205596924, + 6.0464091300964355, + 6.284290313720703, + 6.383488178253174, + 6.532283306121826, + 6.780277252197266, + 7.028269290924072, + 7.27626371383667, + 7.524255752563477, + 7.772250175476074, + 8.040055274963379, + 8.318523406982422], + [5.836524486541748, + 6.067535400390625, + 6.301788806915283, + 6.3954901695251465, + 6.536042213439941, + 6.770293712615967, + 7.00454568862915, + 7.238796710968018, + 7.482921600341797, + 7.740076541900635, + 8.005291938781738, + 8.270505905151367], + [5.86801815032959, + 6.0880255699157715, + 6.308032989501953, + 6.396038055419922, + 6.528041839599609, + 6.749823093414307, + 6.9840264320373535, + 7.218224048614502, + 7.4528584480285645, + 7.70409631729126, + 7.955334663391113, + 8.206572532653809], + [5.892360210418701, + 6.097702980041504, + 6.303044319152832, + 6.390904903411865, + 6.523708343505859, + 6.7450480461120605, + 6.966385841369629, + 7.187726020812988, + 7.417387962341309, + 7.661881923675537, + 7.91318941116333, + 8.164498329162598], + [5.90617036819458, + 6.104805946350098, + 6.312875747680664, + 6.396102428436279, + 6.520944595336914, + 6.729012966156006, + 6.937082290649414, + 7.15734338760376, + 7.39542818069458, + 7.6335129737854, + 7.8715996742248535, + 8.1096830368042]], + 'coords': {'f0': {'dims': ('f0',), + 'attrs': {}, + 'data': [130000.0, + 135000.0, + 140000.0, + 145000.0, + 150000.0, + 155000.0, + 160000.0, + 165000.0, + 375000.0, + 380000.0, + 385000.0, + 390000.0, + 395000.0, + 400000.0, + 405000.0, + 410000.0]}, + 'crosstalk': {'dims': ('crosstalk',), + 'attrs': {}, + 'data': [0.0, + 0.05, + 0.1, + 0.12, + 0.15000000000000002, + 0.2, + 0.25, + 0.30000000000000004, + 0.35000000000000003, + 0.4, + 0.45, + 0.5]}}, + 'name': 'focal_gain'}) def load_transducer_from_file(transducer_filepath : PathLike, convert_array:bool = True) -> Transducer|TransducerArray: """Load a Transducer or TransducerArray from file, depending on the "type" field in the file. @@ -28,3 +260,45 @@ def load_transducer_from_file(transducer_filepath : PathLike, convert_array:bool else: transducer = Transducer.from_file(transducer_filepath) return transducer + +def read_test_report(filename: PathLike) -> pd.DataFrame: + sections = [{"name": "info", "start_row": 3, "nrows": 3}, + {"name": "txm", "start_row": 9, "nrows": 5}, + {"name": "console", "start_row": 17, "nrows": 3}, + {"name": "scans", "start_row": 23, "nrows": 7}, + {"name": "freq", "start_row": 33, "nrows": 9}, + {"name": "voltage", "start_row": 45, "nrows": 7}] + + all_data = [] + for section in sections: + report_df = pd.read_excel(filename, sheet_name="Report", skiprows=section["start_row"], nrows=section["nrows"], index_col=0, usecols="A:C") + report_df["Section"] = section["name"] + all_data.append(report_df) + + report_df = pd.concat(all_data) + return report_df + +def report_to_matrix_dict(report_df: pd.DataFrame, focal_gain_lut=FOCAL_GAIN_LUT) -> dict: + ROW_SN = 'B.1' + ROW_FREQ = 'B.2' + ROW_VOLTAGE = 'E.1' + LIFU_400 = {'id': r'txm_400_{sn}', 'name': r'TXM 400kHz (S/N {sn})', 'nx': 8, 'ny': 8, 'pitch': 5, 'frequency': 400e3, 'kerf': 0.3, 'crosstalk_frac': 0.12, 'crosstalk_dist': 5.05e-3} + LIFU_155 = {'id': r'txm_155_{sn}', 'name': r'TXM 155kHz (S/N {sn})', 'nx': 8, 'ny': 8, 'pitch': 5, 'frequency': 155e3, 'kerf': 0.3, 'crosstalk_frac': 0.12, 'crosstalk_dist': 5.05e-3} + LIFU_MODULES = {400: LIFU_400, 155: LIFU_155} + freq_kHz = report_df.loc[ROW_FREQ]["Value"] + voltage = report_df.loc[ROW_VOLTAGE]["Value"] + sn = report_df.loc[ROW_SN]["Value"] + pattern = r'[^a-zA-Z0-9\-\_]' + replacement = '' + sn = re.sub(pattern, replacement, sn) + matrix_dict = LIFU_MODULES[freq_kHz] + freq_df = report_df[report_df["Section"] == "freq"].copy().drop(columns=["Section"]) + freq_df = freq_df.rename(columns={"Value": "PNP"}) + freq_df = freq_df[freq_df['Item'].str.startswith("PNP")] + freq_df["Frequency"] = freq_df['Item'].apply(lambda x: float(re.search(r"(?<=^PNP \()\d+(?= kHz\)$)", x).group(0))) + freq_df['focal_gain'] = freq_df['Frequency'].apply(lambda f: focal_gain_lut.interp(f0=f*1e3, crosstalk=matrix_dict['crosstalk_frac']).item()) + freq_df['Sensitivity'] = freq_df['PNP'].astype(float)*1e6/freq_df['focal_gain']/voltage + matrix_dict['sensitivity'] = {f*1e3:sens for f, sens, in zip(freq_df["Frequency"], freq_df['Sensitivity'])} + matrix_dict['id'] = matrix_dict['id'].format(sn=sn.lower()) + matrix_dict['name'] = matrix_dict['name'].format(sn=sn) + return matrix_dict From dec962d61acf45bd4afaf19aec13863440518705 Mon Sep 17 00:00:00 2001 From: Peter Hollender Date: Tue, 31 Mar 2026 15:57:54 -0400 Subject: [PATCH 3/9] dynamic reading report improves loading for #440 --- src/openlifu/xdc/util.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/openlifu/xdc/util.py b/src/openlifu/xdc/util.py index c9970c15..e7407fc3 100644 --- a/src/openlifu/xdc/util.py +++ b/src/openlifu/xdc/util.py @@ -262,16 +262,18 @@ def load_transducer_from_file(transducer_filepath : PathLike, convert_array:bool return transducer def read_test_report(filename: PathLike) -> pd.DataFrame: - sections = [{"name": "info", "start_row": 3, "nrows": 3}, - {"name": "txm", "start_row": 9, "nrows": 5}, - {"name": "console", "start_row": 17, "nrows": 3}, - {"name": "scans", "start_row": 23, "nrows": 7}, - {"name": "freq", "start_row": 33, "nrows": 9}, - {"name": "voltage", "start_row": 45, "nrows": 7}] - + sections = [{"name": "info", "start_row": "A"}, + {"name": "txm", "start_row": "B"}, + {"name": "console", "start_row": "C"}, + {"name": "scans", "start_row": "D"}, + {"name": "freq", "start_row": "E"}, + {"name": "voltage", "start_row": "F"}] + raw = pd.read_excel(filename, sheet_name="Report", header=None, usecols="A").rename({0: "Index"}, axis=1) all_data = [] for section in sections: - report_df = pd.read_excel(filename, sheet_name="Report", skiprows=section["start_row"], nrows=section["nrows"], index_col=0, usecols="A:C") + skiprows = raw.loc[raw["Index"] == section["start_row"]].index[0]+1 + nrows = raw['Index'].str.startswith(f'{section["start_row"]}.').sum() + report_df = pd.read_excel(filename, sheet_name="Report", skiprows=skiprows, nrows=nrows, index_col=0, usecols="A:C") report_df["Section"] = section["name"] all_data.append(report_df) From 71b07ac0e5fd0c277d00a4c081c1de39b4aac623 Mon Sep 17 00:00:00 2001 From: Peter Hollender Date: Mon, 6 Apr 2026 16:18:56 -0400 Subject: [PATCH 4/9] fix some frequency-depended issues (#439) --- src/openlifu/plan/solution.py | 9 ++--- src/openlifu/sim/kwave_if.py | 17 +++++---- src/openlifu/xdc/element.py | 70 +++++++++++++++++++++++++++------- src/openlifu/xdc/transducer.py | 4 +- 4 files changed, 73 insertions(+), 27 deletions(-) diff --git a/src/openlifu/plan/solution.py b/src/openlifu/plan/solution.py index a206378e..0ab43b5c 100644 --- a/src/openlifu/plan/solution.py +++ b/src/openlifu/plan/solution.py @@ -256,15 +256,14 @@ def analyze(self, apodization = self.apodizations[focus_index] origin = self.transducer.get_effective_origin(apodizations=apodization, units=options.distance_units) - output_signal_Pa = self.transducer.calc_output( - input_signal_V, - cycles=self.pulse.duration * self.pulse.frequency, + p0_Pa = np.max(self.transducer.calc_output( + cycles=1, frequency=self.pulse.frequency, dt=dt, delays=self.delays[focus_index, :], apod=self.apodizations[focus_index, :], - ) - p0_Pa = np.max(output_signal_Pa, axis=1) + amplitude=self.pulse.amplitude * self.voltage, + ), axis=1) mainlobe_mask = get_mask( pnp_MPa, diff --git a/src/openlifu/sim/kwave_if.py b/src/openlifu/sim/kwave_if.py index 6687ea03..3c87719c 100644 --- a/src/openlifu/sim/kwave_if.py +++ b/src/openlifu/sim/kwave_if.py @@ -139,7 +139,8 @@ def run_simulation(arr: xdc.Transducer, apod = apod if apod is not None else np.ones(arr.numelements()) kgrid = get_kgrid(params.coords, dt=dt, t_end=t_end, cfl=cfl) t = np.arange(0, np.min([cycles / freq, (kgrid.Nt-np.ceil(max(delays)/kgrid.dt))*kgrid.dt]), kgrid.dt) - input_signal = amplitude * np.sin(2 * np.pi * freq * t) + if cycles/freq > t[-1]: + cycles = int(t[-1]*freq) units = [params[dim].attrs['units'] for dim in params.dims] if not all(unit == units[0] for unit in units): raise ValueError("All dimensions must have the same units") @@ -156,14 +157,15 @@ def run_simulation(arr: xdc.Transducer, if _source is not None: source = _source else: - source_mat = arr.calc_output(input_signal, cycles=cycles, frequency=freq, dt=kgrid.dt, delays=delays, apod=apod) + source_mat = arr.calc_output(amplitude=amplitude, cycles=cycles, frequency=freq, dt=kgrid.dt, delays=delays, apod=apod) + source_mat = source_mat[:, :kgrid.Nt] if arr.crosstalk_frac != 0: # Simulate crosstalk by adding additional elements to the array for each element that # is within the crosstalk distance, with the signal scaled by the crosstalk fraction. # This is a simple model of crosstalk and may not capture all the complexities of real # crosstalk, but it can be useful for testing and simulation purposes. - crosstalk_arr = arr.copy() - crosstalk_mat = source_mat + #crosstalk_arr = arr.copy() + crosstalk_mat = source_mat.copy() positions = arr.get_positions(units="m") for src_idx in range(arr.numelements()): for dst_idx in range(arr.numelements()): @@ -173,9 +175,10 @@ def run_simulation(arr: xdc.Transducer, dst_pos = np.array(positions[dst_idx]) dist = np.linalg.norm(src_pos - dst_pos) if dist <= arr.crosstalk_dist: - crosstalk_arr.elements += [arr.elements[dst_idx].copy()] - crosstalk_mat = np.vstack((crosstalk_mat, arr.crosstalk_frac*source_mat[src_idx,:])) - arr = crosstalk_arr + # crosstalk_arr.elements += [arr.elements[dst_idx].copy()] + crosstalk_mat[dst_idx,:] += arr.crosstalk_frac*source_mat[src_idx,:] + #crosstalk_mat = np.vstack((crosstalk_mat, arr.crosstalk_frac*source_mat[src_idx,:])) + #arr = crosstalk_arr source_mat = crosstalk_mat karray = get_karray(arr, translation=array_offset, diff --git a/src/openlifu/xdc/element.py b/src/openlifu/xdc/element.py index b961a3ae..f6b58ca5 100644 --- a/src/openlifu/xdc/element.py +++ b/src/openlifu/xdc/element.py @@ -9,16 +9,66 @@ from openlifu.util.annotations import OpenLIFUFieldData from openlifu.util.units import getunitconversion +SENS_FREQ_KEY = "freq_Hz" +SENS_VALUE_KEY = "values_Pa_per_V" -def sensitivity_at_frequency(sensitivity: float | dict[float, float], frequency: float) -> float: + +def normalize_sensitivity(sensitivity: float | dict) -> float | dict[str, list[float]]: + """Normalize sensitivity to a canonical representation. + + Canonical frequency-dependent representation is: + {"freq_Hz": [...], "values_Pa_per_V": [...]} + + Backward-compatible legacy representation with frequency keys is accepted and + converted to the canonical representation. + """ if isinstance(sensitivity, dict): - freqs = np.array(list(sensitivity.keys()), dtype=np.float64) - values = np.array(list(sensitivity.values()), dtype=np.float64) - return float(np.interp(frequency, freqs, values, left=values[0], right=values[-1])) + if SENS_FREQ_KEY in sensitivity or SENS_VALUE_KEY in sensitivity: + if SENS_FREQ_KEY not in sensitivity or SENS_VALUE_KEY not in sensitivity: + raise ValueError("Sensitivity dictionary must include both 'freq_Hz' and 'values_Pa_per_V'.") + freqs = np.asarray(sensitivity[SENS_FREQ_KEY], dtype=np.float64).reshape(-1) + values = np.asarray(sensitivity[SENS_VALUE_KEY], dtype=np.float64).reshape(-1) + else: + # Legacy format: {frequency_hz: sensitivity} + if len(sensitivity) == 0: + raise ValueError("Sensitivity dictionary must not be empty.") + mapping = {float(k): float(v) for k, v in sensitivity.items()} + freqs = np.array(list(mapping.keys()), dtype=np.float64) + values = np.array(list(mapping.values()), dtype=np.float64) + + if len(freqs) == 0: + raise ValueError("Sensitivity frequency list must not be empty.") + if len(freqs) != len(values): + raise ValueError("Sensitivity frequency and value lists must have the same length.") + + order = np.argsort(freqs) + freqs = freqs[order] + values = values[order] + if np.any(np.diff(freqs) <= 0): + raise ValueError("Sensitivity frequencies must be strictly increasing.") + + return { + SENS_FREQ_KEY: [float(f) for f in freqs], + SENS_VALUE_KEY: [float(v) for v in values], + } + + return float(sensitivity) + + +def sensitivity_at_frequency(sensitivity: float | dict, frequency: float) -> float: + sensitivity = normalize_sensitivity(sensitivity) + if isinstance(sensitivity, dict): + if frequency in sensitivity[SENS_FREQ_KEY]: + idx = sensitivity[SENS_FREQ_KEY].index(frequency) + return float(sensitivity[SENS_VALUE_KEY][idx]) + else: + freqs = np.array(sensitivity[SENS_FREQ_KEY], dtype=np.float64) + values = np.array(sensitivity[SENS_VALUE_KEY], dtype=np.float64) + return float(np.interp(frequency, freqs, values, left=values[0], right=values[-1])) return float(sensitivity) -def generate_drive_signal(input_signal, cycles: float, frequency: float, dt: float) -> np.ndarray: +def generate_drive_signal(cycles: float, frequency: float, dt: float, amplitude: float = 1.0) -> np.ndarray: """Generate a drive signal with duration constrained by cycles/frequency.""" if dt <= 0: raise ValueError("dt must be positive.") @@ -27,14 +77,8 @@ def generate_drive_signal(input_signal, cycles: float, frequency: float, dt: flo if cycles <= 0: raise ValueError("cycles must be positive.") n_samples = max(1, int(np.round(cycles / (frequency * dt)))) - if np.isscalar(input_signal): - t = np.arange(n_samples, dtype=np.float64) * dt - return float(input_signal) * np.sin(2 * np.pi * frequency * t) - base = np.asarray(input_signal, dtype=np.float64).reshape(-1) - drive_signal = np.zeros(n_samples, dtype=np.float64) - n_copy = min(n_samples, len(base)) - drive_signal[:n_copy] = base[:n_copy] - return drive_signal + t = np.arange(n_samples, dtype=np.float64) * dt + return amplitude * np.sin(2 * np.pi * frequency * t) def matrix2xyz(matrix): diff --git a/src/openlifu/xdc/transducer.py b/src/openlifu/xdc/transducer.py index ebeaaaca..311f304e 100644 --- a/src/openlifu/xdc/transducer.py +++ b/src/openlifu/xdc/transducer.py @@ -112,12 +112,12 @@ def __post_init__(self): def get_sensitivity(self, frequency: float) -> float: return sensitivity_at_frequency(self.sensitivity, frequency) - def calc_output(self, input_signal, cycles: float, frequency: float, dt: float, delays: np.ndarray = None, apod: np.ndarray = None): + def calc_output(self, cycles: float, frequency: float, dt: float, delays: np.ndarray = None, apod: np.ndarray = None, amplitude: float = 1.0) -> np.ndarray: if delays is None: delays = np.zeros(self.numelements()) if apod is None: apod = np.ones(self.numelements()) - drive_signal = generate_drive_signal(input_signal, cycles=cycles, frequency=frequency, dt=dt) + drive_signal = generate_drive_signal(cycles=cycles, frequency=frequency, dt=dt, amplitude=amplitude) base_output = drive_signal * self.get_sensitivity(frequency) outputs = [ np.concatenate( From 2956efd6d88aadc44d02155676cc1d0b007e9e83 Mon Sep 17 00:00:00 2001 From: Peter Hollender Date: Fri, 3 Apr 2026 16:22:19 -0400 Subject: [PATCH 5/9] adjust frequency-dependent sensitivity (#439) --- src/openlifu/xdc/element.py | 13 ++------ src/openlifu/xdc/transducer.py | 59 +++++++++++++++++++++------------- src/openlifu/xdc/util.py | 5 ++- tests/test_transducer.py | 28 ++++++++++------ 4 files changed, 61 insertions(+), 44 deletions(-) diff --git a/src/openlifu/xdc/element.py b/src/openlifu/xdc/element.py index f6b58ca5..65b451ce 100644 --- a/src/openlifu/xdc/element.py +++ b/src/openlifu/xdc/element.py @@ -114,7 +114,7 @@ class Element: size: Annotated[np.ndarray, OpenLIFUFieldData("Size", "Size of the element in 2D")] = field(default_factory=lambda: np.array([1., 1.])) """ Size of the element in 2D as a numpy array [width, length].""" - sensitivity: Annotated[float | dict[float, float], OpenLIFUFieldData("Sensitivity", "Sensitivity of the element (Pa/V), scalar or {frequency_hz: sensitivity}")] = 1.0 + sensitivity: Annotated[float | dict, OpenLIFUFieldData("Sensitivity", "Sensitivity of the element (Pa/V), scalar or {'freq_Hz':[...], 'values_Pa_per_V':[...]}")] = 1.0 """Sensitivity of the element (Pa/V)""" pin: Annotated[int, OpenLIFUFieldData("Pin", "Channel pin to which the element is connected")] = -1 @@ -135,16 +135,7 @@ def __post_init__(self): raise ValueError("Size must be a 2-element array.") if self.sensitivity is None: self.sensitivity = 1.0 - elif isinstance(self.sensitivity, dict): - if len(self.sensitivity) == 0: - raise ValueError("Sensitivity dictionary must not be empty.") - mapping = {float(k): float(v) for k, v in self.sensitivity.items()} - freqs = np.array(sorted(mapping.keys()), dtype=np.float64) - if np.any(np.diff(freqs) <= 0): - raise ValueError("Sensitivity dictionary frequencies must be strictly increasing.") - self.sensitivity = {float(f): mapping[float(f)] for f in freqs} - else: - self.sensitivity = float(self.sensitivity) + self.sensitivity = normalize_sensitivity(self.sensitivity) @property def x(self): diff --git a/src/openlifu/xdc/transducer.py b/src/openlifu/xdc/transducer.py index 311f304e..e67f2d9c 100644 --- a/src/openlifu/xdc/transducer.py +++ b/src/openlifu/xdc/transducer.py @@ -14,6 +14,7 @@ from openlifu.xdc.element import ( Element, generate_drive_signal, + normalize_sensitivity, sensitivity_at_frequency, ) @@ -22,21 +23,37 @@ def _combine_sensitivities( - base_sensitivity: float | dict[float, float], - scale_sensitivity: float | dict[float, float], -) -> float | dict[float, float]: + base_sensitivity: float | dict, + scale_sensitivity: float | dict, +) -> float | dict[str, list[float]]: + base_sensitivity = normalize_sensitivity(base_sensitivity) + scale_sensitivity = normalize_sensitivity(scale_sensitivity) + if isinstance(base_sensitivity, dict) and isinstance(scale_sensitivity, dict): - base_keys = tuple(base_sensitivity.keys()) - scale_keys = tuple(scale_sensitivity.keys()) - if base_keys != scale_keys: + base_freqs = np.asarray(base_sensitivity["freq_Hz"], dtype=np.float64) + scale_freqs = np.asarray(scale_sensitivity["freq_Hz"], dtype=np.float64) + if not np.array_equal(base_freqs, scale_freqs): raise ValueError("Cannot combine sensitivity dictionaries with different frequency keys.") - return {frequency: base_sensitivity[frequency] * scale_sensitivity[frequency] for frequency in base_keys} + base_values = np.asarray(base_sensitivity["values_Pa_per_V"], dtype=np.float64) + scale_values = np.asarray(scale_sensitivity["values_Pa_per_V"], dtype=np.float64) + return { + "freq_Hz": [float(f) for f in base_freqs], + "values_Pa_per_V": [float(v) for v in (base_values * scale_values)], + } if isinstance(base_sensitivity, dict): factor = float(scale_sensitivity) - return {frequency: value * factor for frequency, value in base_sensitivity.items()} + values = np.asarray(base_sensitivity["values_Pa_per_V"], dtype=np.float64) + return { + "freq_Hz": [float(f) for f in base_sensitivity["freq_Hz"]], + "values_Pa_per_V": [float(v) for v in (values * factor)], + } if isinstance(scale_sensitivity, dict): factor = float(base_sensitivity) - return {frequency: factor * value for frequency, value in scale_sensitivity.items()} + values = np.asarray(scale_sensitivity["values_Pa_per_V"], dtype=np.float64) + return { + "freq_Hz": [float(f) for f in scale_sensitivity["freq_Hz"]], + "values_Pa_per_V": [float(v) for v in (factor * values)], + } return float(base_sensitivity) * float(scale_sensitivity) @dataclass @@ -78,7 +95,7 @@ class Transducer: The units of this transform are assumed to be the native units of the transducer, the `Transducer.units` field. """ - sensitivity: Annotated[float | dict[float, float], OpenLIFUFieldData("Sensitivity", "Sensitivity of the element (Pa/V), scalar or {frequency_hz: sensitivity}")] = 1.0 + sensitivity: Annotated[float | dict, OpenLIFUFieldData("Sensitivity", "Sensitivity of the transducer (Pa/V), scalar or {'freq_Hz':[...], 'values_Pa_per_V':[...]}")] = 1.0 """Sensitivity of the transducer (Pa/V), scalar or frequency-dependent dictionary.""" crosstalk_frac: Annotated[float, OpenLIFUFieldData("Crosstalk fraction", "Fraction of the signal that leaks into other elements due to crosstalk")] = 0.0 @@ -98,16 +115,7 @@ def __post_init__(self): element.rescale(self.units) if self.sensitivity is None: self.sensitivity = 1.0 - elif isinstance(self.sensitivity, dict): - if len(self.sensitivity) == 0: - raise ValueError("Sensitivity dictionary must not be empty.") - mapping = {float(k): float(v) for k, v in self.sensitivity.items()} - freqs = np.array(sorted(mapping.keys()), dtype=np.float64) - if np.any(np.diff(freqs) <= 0): - raise ValueError("Sensitivity dictionary frequencies must be strictly increasing.") - self.sensitivity = {float(f): mapping[float(f)] for f in freqs} - else: - self.sensitivity = float(self.sensitivity) + self.sensitivity = normalize_sensitivity(self.sensitivity) def get_sensitivity(self, frequency: float) -> float: return sensitivity_at_frequency(self.sensitivity, frequency) @@ -252,18 +260,23 @@ def merge(list_of_transducers:List[Transducer], offset_pins:bool=False, offset_i array_copies = [arr.copy() for arr in list_of_transducers] dict_key_sets = set() for array in array_copies: + array.sensitivity = normalize_sensitivity(array.sensitivity) if isinstance(array.sensitivity, dict): - dict_key_sets.add(tuple(array.sensitivity.keys())) + dict_key_sets.add(tuple(array.sensitivity["freq_Hz"])) for el in array.elements: + el.sensitivity = normalize_sensitivity(el.sensitivity) if isinstance(el.sensitivity, dict): - dict_key_sets.add(tuple(el.sensitivity.keys())) + dict_key_sets.add(tuple(el.sensitivity["freq_Hz"])) if len(dict_key_sets) > 1: raise ValueError("Cannot merge sensitivities with different frequency keys.") sensitivity_signatures = [] for array in array_copies: if isinstance(array.sensitivity, dict): - sensitivity_signatures.append((tuple(array.sensitivity.keys()), tuple(array.sensitivity.values()))) + sensitivity_signatures.append(( + tuple(array.sensitivity["freq_Hz"]), + tuple(array.sensitivity["values_Pa_per_V"]), + )) else: sensitivity_signatures.append(float(array.sensitivity)) diff --git a/src/openlifu/xdc/util.py b/src/openlifu/xdc/util.py index e7407fc3..21f84c84 100644 --- a/src/openlifu/xdc/util.py +++ b/src/openlifu/xdc/util.py @@ -300,7 +300,10 @@ def report_to_matrix_dict(report_df: pd.DataFrame, focal_gain_lut=FOCAL_GAIN_LUT freq_df["Frequency"] = freq_df['Item'].apply(lambda x: float(re.search(r"(?<=^PNP \()\d+(?= kHz\)$)", x).group(0))) freq_df['focal_gain'] = freq_df['Frequency'].apply(lambda f: focal_gain_lut.interp(f0=f*1e3, crosstalk=matrix_dict['crosstalk_frac']).item()) freq_df['Sensitivity'] = freq_df['PNP'].astype(float)*1e6/freq_df['focal_gain']/voltage - matrix_dict['sensitivity'] = {f*1e3:sens for f, sens, in zip(freq_df["Frequency"], freq_df['Sensitivity'])} + matrix_dict['sensitivity'] = { + 'freq_Hz': [float(f * 1e3) for f in freq_df["Frequency"]], + 'values_Pa_per_V': [float(sens) for sens in freq_df['Sensitivity']], + } matrix_dict['id'] = matrix_dict['id'].format(sn=sn.lower()) matrix_dict['name'] = matrix_dict['name'].format(sn=sn) return matrix_dict diff --git a/tests/test_transducer.py b/tests/test_transducer.py index cc79872f..e3cdc0f1 100644 --- a/tests/test_transducer.py +++ b/tests/test_transducer.py @@ -116,7 +116,7 @@ def test_transducer_calc_output_interpolates_dictionary_sensitivity(): nx=1, ny=1, units="mm", - sensitivity={100e3: 1.0, 300e3: 3.0}, + sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [1.0, 3.0]}, ) transducer.elements[0].sensitivity = 1.0 input_signal = np.array([1.0, -1.0, 0.5], dtype=float) @@ -128,6 +128,16 @@ def test_transducer_calc_output_interpolates_dictionary_sensitivity(): np.testing.assert_allclose(output_low[0, :len(input_signal)], 1.0 * input_signal) +def test_legacy_sensitivity_mapping_is_normalized_to_schema(): + transducer = Transducer( + sensitivity={"100000.0": 1.0, "300000.0": 3.0}, + ) + assert transducer.sensitivity == { + "freq_Hz": [100000.0, 300000.0], + "values_Pa_per_V": [1.0, 3.0], + } + + def test_element_calc_output_generates_signal_from_scalar_input(): element = Element(sensitivity=2.0) cycles = 4 @@ -161,13 +171,13 @@ def test_merge_pushes_transducer_sensitivity_into_elements(): nx=1, ny=1, units="mm", - sensitivity={100e3: 2.0, 300e3: 4.0}, + sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [2.0, 4.0]}, ) transducer_b = Transducer.gen_matrix_array( nx=1, ny=1, units="mm", - sensitivity={100e3: 3.0, 300e3: 6.0}, + sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [3.0, 6.0]}, ) transducer_a.elements[0].sensitivity = 5.0 transducer_b.elements[0].sensitivity = 7.0 @@ -175,8 +185,8 @@ def test_merge_pushes_transducer_sensitivity_into_elements(): merged = Transducer.merge([transducer_a, transducer_b], merge_mismatched_sensitivity=True) assert merged.sensitivity == 1.0 - assert merged.elements[0].sensitivity == {100e3: 10.0, 300e3: 20.0} - assert merged.elements[1].sensitivity == {100e3: 21.0, 300e3: 42.0} + assert merged.elements[0].sensitivity == {"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [10.0, 20.0]} + assert merged.elements[1].sensitivity == {"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [21.0, 42.0]} def test_merge_rejects_mismatched_sensitivity_keys(): @@ -184,16 +194,16 @@ def test_merge_rejects_mismatched_sensitivity_keys(): nx=1, ny=1, units="mm", - sensitivity={100e3: 2.0, 300e3: 4.0}, + sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [2.0, 4.0]}, ) transducer_b = Transducer.gen_matrix_array( nx=1, ny=1, units="mm", - sensitivity={100e3: 3.0, 400e3: 6.0}, + sensitivity={"freq_Hz": [100e3, 400e3], "values_Pa_per_V": [3.0, 6.0]}, ) - transducer_a.elements[0].sensitivity = {100e3: 5.0, 300e3: 7.0} - transducer_b.elements[0].sensitivity = {100e3: 11.0, 400e3: 13.0} + transducer_a.elements[0].sensitivity = {"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [5.0, 7.0]} + transducer_b.elements[0].sensitivity = {"freq_Hz": [100e3, 400e3], "values_Pa_per_V": [11.0, 13.0]} with pytest.raises(ValueError, match="different frequency keys"): Transducer.merge([transducer_a, transducer_b], merge_mismatched_sensitivity=True) From e00eaba25e028e5c745edd7c2ccb14cbcbc6f276 Mon Sep 17 00:00:00 2001 From: Peter Hollender Date: Mon, 6 Apr 2026 15:30:29 -0400 Subject: [PATCH 6/9] Clean up k-wave files #435 Closes #435 --- src/openlifu/sim/kwave_if.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/openlifu/sim/kwave_if.py b/src/openlifu/sim/kwave_if.py index 3c87719c..5e95b9ea 100644 --- a/src/openlifu/sim/kwave_if.py +++ b/src/openlifu/sim/kwave_if.py @@ -2,6 +2,7 @@ import logging from copy import deepcopy +from pathlib import Path from typing import List import numpy as np @@ -222,6 +223,11 @@ def run_simulation(arr: xdc.Transducer, coords=[pcoords[dim] for dim in ['t','x','y','z']], attrs={'units':'Pa', 'long_name':'Pressure'}) + # clean up temporary files created by k-wave + for filename in [simulation_options.input_filename, simulation_options.output_filename]: + Path(filename).unlink(missing_ok=True) + + ds = xa.Dataset(ds_dict) if return_kwave_outputs and return_kwave_inputs: return ds, output, inputs From aea329661b7ba11d0d11235c583ac1bf1174005a Mon Sep 17 00:00:00 2001 From: Peter Hollender Date: Tue, 7 Apr 2026 11:44:59 -0400 Subject: [PATCH 7/9] Fix bug in `element.calc_output` #439 --- src/openlifu/xdc/element.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/openlifu/xdc/element.py b/src/openlifu/xdc/element.py index 65b451ce..f410d54e 100644 --- a/src/openlifu/xdc/element.py +++ b/src/openlifu/xdc/element.py @@ -204,9 +204,9 @@ def length(self, value): def get_sensitivity(self, frequency: float) -> float: return sensitivity_at_frequency(self.sensitivity, frequency) - def calc_output(self, input_signal, cycles: float, frequency: float, dt: float): - drive_signal = generate_drive_signal(input_signal, cycles=cycles, frequency=frequency, dt=dt) - return drive_signal * self.get_sensitivity(frequency) + def calc_output(self, cycles: float, frequency: float, dt: float, amplitude: float = 1.0) -> np.ndarray: + drive_signal = generate_drive_signal(cycles=cycles, frequency=frequency, dt=dt) + return drive_signal * self.get_sensitivity(frequency) * amplitude def copy(self): return copy.deepcopy(self) From f1423637e59f42d3330f3c512c63a90bebf0ea6e Mon Sep 17 00:00:00 2001 From: Peter Hollender Date: Tue, 7 Apr 2026 14:15:11 -0400 Subject: [PATCH 8/9] Clean up transducerarray to move elegantly handle differently-sensitive modules #439 --- src/openlifu/sim/kwave_if.py | 2 +- src/openlifu/xdc/transducerarray.py | 49 ++++++-- tests/test_transducer.py | 175 ++++++++++++++++++++++++++-- 3 files changed, 205 insertions(+), 21 deletions(-) diff --git a/src/openlifu/sim/kwave_if.py b/src/openlifu/sim/kwave_if.py index 5e95b9ea..6f2f2e7f 100644 --- a/src/openlifu/sim/kwave_if.py +++ b/src/openlifu/sim/kwave_if.py @@ -141,7 +141,7 @@ def run_simulation(arr: xdc.Transducer, kgrid = get_kgrid(params.coords, dt=dt, t_end=t_end, cfl=cfl) t = np.arange(0, np.min([cycles / freq, (kgrid.Nt-np.ceil(max(delays)/kgrid.dt))*kgrid.dt]), kgrid.dt) if cycles/freq > t[-1]: - cycles = int(t[-1]*freq) + cycles = np.ceil(t[-1]*freq) units = [params[dim].attrs['units'] for dim in params.dims] if not all(unit == units[0] for unit in units): raise ValueError("All dimensions must have the same units") diff --git a/src/openlifu/xdc/transducerarray.py b/src/openlifu/xdc/transducerarray.py index 9eaaf978..e552f20b 100644 --- a/src/openlifu/xdc/transducerarray.py +++ b/src/openlifu/xdc/transducerarray.py @@ -23,6 +23,15 @@ def get_angle_from_gap(width, gap, roc): def get_roc_from_angle(width, gap, dth): return (0.5*gap + (0.5 * width * np.cos(dth))) / np.sin(dth) +def get_gap_from_angle(width, dth, roc): + a = roc + b = width/2 + mag = np.sqrt(a**2 + b**2) + A = a/mag + B = b/mag + gap = 2*mag*np.sin(dth - np.arcsin(B)) + return gap if A >= 0 else -gap + @dataclass class TransducerArray(DictMixin): id: str = "transducer_array" @@ -84,22 +93,44 @@ def to_file(self, file_path: str, compact: bool = False) -> None: f.write(json_string) @staticmethod - def get_concave_cylinder(trans, rows=1, cols=1, width=40, gap=0, dth=None, roc=np.inf, units="mm", id="transducer_array", name="Transducer Array", attrs: dict={}): - scl = getunitconversion(units, trans.units) + def get_concave_cylinder(trans, rows=1, cols=1, width=40, gap=None, dth=None, roc=None, units="mm", id="transducer_array", name="Transducer Array", attrs: dict={}): + modules = [] - if roc == np.inf: + if isinstance(trans, Transducer): + trans_arr = np.array([[trans]*cols for _ in range(rows)]) + else: + trans_arr = np.array(trans).reshape(rows, cols) + scl = getunitconversion(units, trans_arr[0,0].units) + if gap is None: + if dth is not None and roc is not None: + gap = get_gap_from_angle(width, dth, roc) + else: + gap = 0 + elif dth is not None and roc is not None: + raise ValueError("Invalid combination of parameters: cannot specify all of gap, dth, and roc.") + + if dth is None: + if roc is not None: + dth = get_angle_from_gap(width, gap, roc) + else: + dth = 0 + roc = np.inf + + if roc is None: + if np.isclose(dth, 0.0): + roc = np.inf + else: + roc = get_roc_from_angle(width, gap, dth) + + if dth == 0: for i in range(rows): y = (width+gap)*(i-(rows-1)/2)*scl for j in range(cols): dx = (width+gap)*(j-(cols-1)/2)*scl M = np.array([[1,0,0,dx], [0,1,0,y], [0,0,1,0], [0,0,0,1]]) - trans_new = TransformedTransducer.from_transducer(trans, transform=np.linalg.inv(M)) + trans_new = TransformedTransducer.from_transducer(trans_arr[i,j], transform=np.linalg.inv(M)) modules.append(trans_new) else: - if dth is None: - dth = get_angle_from_gap(width, gap, roc) - elif roc is None and dth is not None and gap is not None: - roc = get_roc_from_angle(width, gap, dth) for i in range(rows): y = (width+gap)*(i-(rows-1)/2)*scl for j in range(cols): @@ -110,7 +141,7 @@ def get_concave_cylinder(trans, rows=1, cols=1, width=40, gap=0, dth=None, roc=n [0,1,0,y], [np.sin(th),0,np.cos(th),z], [0,0,0,1]]) - trans_new = TransformedTransducer.from_transducer(trans, transform=np.linalg.inv(M)) + trans_new = TransformedTransducer.from_transducer(trans_arr[i,j], transform=np.linalg.inv(M)) modules.append(trans_new) return TransducerArray(modules=modules, id=id, name=name, attrs=attrs) diff --git a/tests/test_transducer.py b/tests/test_transducer.py index e3cdc0f1..417b4ad2 100644 --- a/tests/test_transducer.py +++ b/tests/test_transducer.py @@ -7,6 +7,11 @@ from helpers import dataclasses_are_equal from openlifu.xdc import Element, Transducer, TransducerArray +from openlifu.xdc.transducerarray import ( + get_angle_from_gap, + get_gap_from_angle, + get_roc_from_angle, +) @pytest.fixture() @@ -119,13 +124,21 @@ def test_transducer_calc_output_interpolates_dictionary_sensitivity(): sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [1.0, 3.0]}, ) transducer.elements[0].sensitivity = 1.0 - input_signal = np.array([1.0, -1.0, 0.5], dtype=float) + cycles = 3 + dt = 1e-7 - output_mid = transducer.calc_output(input_signal, cycles=3, frequency=200e3, dt=1e-7) - output_low = transducer.calc_output(input_signal, cycles=3, frequency=100e3, dt=1e-7) + output_mid = transducer.calc_output(cycles=cycles, frequency=200e3, dt=dt) + output_low = transducer.calc_output(cycles=cycles, frequency=100e3, dt=dt) - np.testing.assert_allclose(output_mid[0, :len(input_signal)], 2.0 * input_signal) - np.testing.assert_allclose(output_low[0, :len(input_signal)], 1.0 * input_signal) + n_samples_mid = int(np.round(cycles / (200e3 * dt))) + n_samples_low = int(np.round(cycles / (100e3 * dt))) + t_mid = np.arange(n_samples_mid) * dt + t_low = np.arange(n_samples_low) * dt + expected_mid = 2.0 * np.sin(2 * np.pi * 200e3 * t_mid) + expected_low = 1.0 * np.sin(2 * np.pi * 100e3 * t_low) + + np.testing.assert_allclose(output_mid[0], expected_mid) + np.testing.assert_allclose(output_low[0], expected_low) def test_legacy_sensitivity_mapping_is_normalized_to_schema(): @@ -145,25 +158,25 @@ def test_element_calc_output_generates_signal_from_scalar_input(): dt = 1e-7 n_samples = int(np.round(cycles / (frequency * dt))) - output = element.calc_output(3.0, cycles=cycles, frequency=frequency, dt=dt) + output = element.calc_output(cycles=cycles, frequency=frequency, dt=dt, amplitude=3.0) t = np.arange(n_samples) * dt expected = 2.0 * 3.0 * np.sin(2 * np.pi * frequency * t) np.testing.assert_allclose(output, expected) -def test_element_calc_output_enforces_cycles_duration_for_array_input(): +def test_element_calc_output_enforces_cycles_duration_for_generated_signal(): element = Element(sensitivity=1.0) cycles = 1 frequency = 200e3 dt = 1e-6 n_samples = int(np.round(cycles / (frequency * dt))) - input_signal = np.arange(20, dtype=float) - - output = element.calc_output(input_signal, cycles=cycles, frequency=frequency, dt=dt) + output = element.calc_output(cycles=cycles, frequency=frequency, dt=dt) + t = np.arange(n_samples) * dt + expected = np.sin(2 * np.pi * frequency * t) assert len(output) == n_samples - np.testing.assert_allclose(output, input_signal[:n_samples]) + np.testing.assert_allclose(output, expected) def test_merge_pushes_transducer_sensitivity_into_elements(): @@ -207,3 +220,143 @@ def test_merge_rejects_mismatched_sensitivity_keys(): with pytest.raises(ValueError, match="different frequency keys"): Transducer.merge([transducer_a, transducer_b], merge_mismatched_sensitivity=True) + + +@pytest.mark.parametrize( + ("width", "dth", "roc"), + [ + (8.0, 0.08, 25.0), + (10.0, 0.12, 30.0), + (12.0, 0.18, 45.0), + ], +) +def test_concave_geometry_helpers_are_mutual_inverses(width: float, dth: float, roc: float): + gap = get_gap_from_angle(width, dth, roc) + recovered_roc = get_roc_from_angle(width, gap, dth) + recovered_dth = get_angle_from_gap(width, gap, roc) + recovered_gap = get_gap_from_angle(width, recovered_dth, roc) + + assert np.isclose(recovered_roc, roc) + assert np.isclose(recovered_dth, dth) + assert np.isclose(recovered_gap, gap) + + +def test_get_concave_cylinder_computes_gap_from_dth_and_roc_layout_spacing(): + base = Transducer.gen_matrix_array(nx=1, ny=1, units="mm") + width = 8.0 + dth = 0.12 + roc = 25.0 + array = TransducerArray.get_concave_cylinder( + base, + rows=2, + cols=1, + width=width, + dth=dth, + roc=roc, + units="mm", + ) + merged = array.to_transducer() + positions = merged.get_positions(units="mm") + + expected_gap = get_gap_from_angle(width, dth, roc) + y_spacing = np.abs(positions[1, 1] - positions[0, 1]) + + assert np.isclose(y_spacing, width + expected_gap) + + +def test_get_concave_cylinder_handles_zero_dth_without_roc(): + base = Transducer.gen_matrix_array(nx=1, ny=1, units="mm") + width = 10.0 + gap = 2.0 + array = TransducerArray.get_concave_cylinder( + base, + rows=1, + cols=2, + width=width, + gap=gap, + dth=0.0, + units="mm", + ) + merged = array.to_transducer() + positions = merged.get_positions(units="mm") + + x_spacing = np.abs(positions[1, 0] - positions[0, 0]) + z_values = positions[:, 2] + + assert np.isclose(x_spacing, width + gap) + np.testing.assert_allclose(z_values, np.zeros_like(z_values)) + + +def test_get_concave_cylinder_rejects_gap_dth_roc_together(): + base = Transducer.gen_matrix_array(nx=1, ny=1, units="mm") + with pytest.raises(ValueError, match="cannot specify all of gap, dth, and roc"): + TransducerArray.get_concave_cylinder( + base, + rows=1, + cols=2, + width=10.0, + gap=1.0, + dth=0.2, + roc=20.0, + units="mm", + ) + + +def test_transducer_calc_output_combines_frequency_dependent_sensitivities(): + transducer = Transducer.gen_matrix_array( + nx=1, + ny=1, + units="mm", + sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [2.0, 4.0]}, + ) + transducer.elements[0].sensitivity = {"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [5.0, 9.0]} + + frequency = 200e3 + dt = 1e-7 + cycles = 3 + n_samples = int(np.round(cycles / (frequency * dt))) + t = np.arange(n_samples) * dt + expected_drive = np.sin(2 * np.pi * frequency * t) + + output = transducer.calc_output(cycles=cycles, frequency=frequency, dt=dt) + + np.testing.assert_allclose(output[0], 21.0 * expected_drive) + + +def test_transducer_array_to_transducer_preserves_frequency_dependent_sensitivities(): + transducer_a = Transducer.gen_matrix_array( + nx=1, + ny=1, + units="mm", + sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [2.0, 4.0]}, + ) + transducer_b = Transducer.gen_matrix_array( + nx=1, + ny=1, + units="mm", + sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [1.0, 3.0]}, + ) + transducer_a.elements[0].sensitivity = 5.0 + transducer_b.elements[0].sensitivity = 7.0 + + array = TransducerArray.get_concave_cylinder( + [transducer_a, transducer_b], + rows=1, + cols=2, + width=10.0, + gap=0.0, + units="mm", + ) + merged = array.to_transducer() + + frequency = 200e3 + dt = 1e-7 + cycles = 2 + n_samples = int(np.round(cycles / (frequency * dt))) + t = np.arange(n_samples) * dt + expected_drive = np.sin(2 * np.pi * frequency * t) + + output = merged.calc_output(cycles=cycles, frequency=frequency, dt=dt) + + np.testing.assert_allclose(output[0], 15.0 * expected_drive) + np.testing.assert_allclose(output[1], 14.0 * expected_drive) From a9772307a297407f66125207939f8af57430d172 Mon Sep 17 00:00:00 2001 From: Peter Hollender Date: Thu, 9 Apr 2026 13:47:43 -0400 Subject: [PATCH 9/9] Update to list of tuples form #439 --- src/openlifu/xdc/element.py | 68 +++++----------------------- src/openlifu/xdc/transducer.py | 77 +++++++++++++------------------- src/openlifu/xdc/util.py | 5 +-- tests/test_transducer.py | 81 ++++++++++++++++++++++++---------- 4 files changed, 102 insertions(+), 129 deletions(-) diff --git a/src/openlifu/xdc/element.py b/src/openlifu/xdc/element.py index f410d54e..1d5bc973 100644 --- a/src/openlifu/xdc/element.py +++ b/src/openlifu/xdc/element.py @@ -2,68 +2,23 @@ import copy from dataclasses import dataclass, field -from typing import Annotated +from typing import Annotated, List import numpy as np from openlifu.util.annotations import OpenLIFUFieldData from openlifu.util.units import getunitconversion -SENS_FREQ_KEY = "freq_Hz" -SENS_VALUE_KEY = "values_Pa_per_V" - -def normalize_sensitivity(sensitivity: float | dict) -> float | dict[str, list[float]]: - """Normalize sensitivity to a canonical representation. - - Canonical frequency-dependent representation is: - {"freq_Hz": [...], "values_Pa_per_V": [...]} - - Backward-compatible legacy representation with frequency keys is accepted and - converted to the canonical representation. - """ - if isinstance(sensitivity, dict): - if SENS_FREQ_KEY in sensitivity or SENS_VALUE_KEY in sensitivity: - if SENS_FREQ_KEY not in sensitivity or SENS_VALUE_KEY not in sensitivity: - raise ValueError("Sensitivity dictionary must include both 'freq_Hz' and 'values_Pa_per_V'.") - freqs = np.asarray(sensitivity[SENS_FREQ_KEY], dtype=np.float64).reshape(-1) - values = np.asarray(sensitivity[SENS_VALUE_KEY], dtype=np.float64).reshape(-1) - else: - # Legacy format: {frequency_hz: sensitivity} - if len(sensitivity) == 0: - raise ValueError("Sensitivity dictionary must not be empty.") - mapping = {float(k): float(v) for k, v in sensitivity.items()} - freqs = np.array(list(mapping.keys()), dtype=np.float64) - values = np.array(list(mapping.values()), dtype=np.float64) - - if len(freqs) == 0: - raise ValueError("Sensitivity frequency list must not be empty.") - if len(freqs) != len(values): - raise ValueError("Sensitivity frequency and value lists must have the same length.") - - order = np.argsort(freqs) - freqs = freqs[order] - values = values[order] - if np.any(np.diff(freqs) <= 0): - raise ValueError("Sensitivity frequencies must be strictly increasing.") - - return { - SENS_FREQ_KEY: [float(f) for f in freqs], - SENS_VALUE_KEY: [float(v) for v in values], - } - - return float(sensitivity) - - -def sensitivity_at_frequency(sensitivity: float | dict, frequency: float) -> float: - sensitivity = normalize_sensitivity(sensitivity) - if isinstance(sensitivity, dict): - if frequency in sensitivity[SENS_FREQ_KEY]: - idx = sensitivity[SENS_FREQ_KEY].index(frequency) - return float(sensitivity[SENS_VALUE_KEY][idx]) +def sensitivity_at_frequency(sensitivity: float | List[tuple[float, float]], frequency: float) -> float: + if isinstance(sensitivity, list): + freqs, values = zip(*sensitivity) + freqs = np.array(freqs, dtype=np.float64) + values = np.array(values, dtype=np.float64) + if frequency in freqs: + idx = np.where(freqs == frequency)[0][0] + return float(values[idx]) else: - freqs = np.array(sensitivity[SENS_FREQ_KEY], dtype=np.float64) - values = np.array(sensitivity[SENS_VALUE_KEY], dtype=np.float64) return float(np.interp(frequency, freqs, values, left=values[0], right=values[-1])) return float(sensitivity) @@ -114,7 +69,7 @@ class Element: size: Annotated[np.ndarray, OpenLIFUFieldData("Size", "Size of the element in 2D")] = field(default_factory=lambda: np.array([1., 1.])) """ Size of the element in 2D as a numpy array [width, length].""" - sensitivity: Annotated[float | dict, OpenLIFUFieldData("Sensitivity", "Sensitivity of the element (Pa/V), scalar or {'freq_Hz':[...], 'values_Pa_per_V':[...]}")] = 1.0 + sensitivity: Annotated[float | List[tuple[float, float]], OpenLIFUFieldData("Sensitivity", "Sensitivity of the element (Pa/V), scalar or list of (frequency, value) tuples")] = 1.0 """Sensitivity of the element (Pa/V)""" pin: Annotated[int, OpenLIFUFieldData("Pin", "Channel pin to which the element is connected")] = -1 @@ -135,7 +90,8 @@ def __post_init__(self): raise ValueError("Size must be a 2-element array.") if self.sensitivity is None: self.sensitivity = 1.0 - self.sensitivity = normalize_sensitivity(self.sensitivity) + elif isinstance(self.sensitivity, list): + self.sensitivity = [(float(f), float(v)) for f, v in self.sensitivity] @property def x(self): diff --git a/src/openlifu/xdc/transducer.py b/src/openlifu/xdc/transducer.py index e67f2d9c..f447d025 100644 --- a/src/openlifu/xdc/transducer.py +++ b/src/openlifu/xdc/transducer.py @@ -14,7 +14,6 @@ from openlifu.xdc.element import ( Element, generate_drive_signal, - normalize_sensitivity, sensitivity_at_frequency, ) @@ -23,38 +22,28 @@ def _combine_sensitivities( - base_sensitivity: float | dict, - scale_sensitivity: float | dict, -) -> float | dict[str, list[float]]: - base_sensitivity = normalize_sensitivity(base_sensitivity) - scale_sensitivity = normalize_sensitivity(scale_sensitivity) - - if isinstance(base_sensitivity, dict) and isinstance(scale_sensitivity, dict): - base_freqs = np.asarray(base_sensitivity["freq_Hz"], dtype=np.float64) - scale_freqs = np.asarray(scale_sensitivity["freq_Hz"], dtype=np.float64) + base_sensitivity: float | List[tuple[float, float]], + scale_sensitivity: float | List[tuple[float, float]], +) -> float | List[tuple[float, float]]: + + if isinstance(base_sensitivity, list) and isinstance(scale_sensitivity, list): + base_freqs = np.asarray([f for f, _ in base_sensitivity], dtype=np.float64) + scale_freqs = np.asarray([f for f, _ in scale_sensitivity], dtype=np.float64) if not np.array_equal(base_freqs, scale_freqs): raise ValueError("Cannot combine sensitivity dictionaries with different frequency keys.") - base_values = np.asarray(base_sensitivity["values_Pa_per_V"], dtype=np.float64) - scale_values = np.asarray(scale_sensitivity["values_Pa_per_V"], dtype=np.float64) - return { - "freq_Hz": [float(f) for f in base_freqs], - "values_Pa_per_V": [float(v) for v in (base_values * scale_values)], - } - if isinstance(base_sensitivity, dict): + base_values = np.asarray([v for _, v in base_sensitivity], dtype=np.float64) + scale_values = np.asarray([v for _, v in scale_sensitivity], dtype=np.float64) + return [(float(f), float(v)) for f, v in zip(base_freqs, base_values * scale_values)] + elif isinstance(base_sensitivity, list): factor = float(scale_sensitivity) - values = np.asarray(base_sensitivity["values_Pa_per_V"], dtype=np.float64) - return { - "freq_Hz": [float(f) for f in base_sensitivity["freq_Hz"]], - "values_Pa_per_V": [float(v) for v in (values * factor)], - } - if isinstance(scale_sensitivity, dict): + values = np.asarray([v for _, v in base_sensitivity], dtype=np.float64) + return [(float(f), float(v)) for (f, _), v in zip(base_sensitivity, values * factor)] + elif isinstance(scale_sensitivity, list): factor = float(base_sensitivity) - values = np.asarray(scale_sensitivity["values_Pa_per_V"], dtype=np.float64) - return { - "freq_Hz": [float(f) for f in scale_sensitivity["freq_Hz"]], - "values_Pa_per_V": [float(v) for v in (factor * values)], - } - return float(base_sensitivity) * float(scale_sensitivity) + values = np.asarray([v for _, v in scale_sensitivity], dtype=np.float64) + return [(float(f), float(v)) for (f, _), v in zip(scale_sensitivity, factor * values)] + else: + return float(base_sensitivity) * float(scale_sensitivity) @dataclass class Transducer: @@ -95,8 +84,8 @@ class Transducer: The units of this transform are assumed to be the native units of the transducer, the `Transducer.units` field. """ - sensitivity: Annotated[float | dict, OpenLIFUFieldData("Sensitivity", "Sensitivity of the transducer (Pa/V), scalar or {'freq_Hz':[...], 'values_Pa_per_V':[...]}")] = 1.0 - """Sensitivity of the transducer (Pa/V), scalar or frequency-dependent dictionary.""" + sensitivity: Annotated[float | List[tuple[float, float]], OpenLIFUFieldData("Sensitivity", "Sensitivity of the transducer (Pa/V), scalar or list of (frequency, value) tuples")] = 1.0 + """Sensitivity of the transducer (Pa/V), scalar or frequency-dependent list of tuples.""" crosstalk_frac: Annotated[float, OpenLIFUFieldData("Crosstalk fraction", "Fraction of the signal that leaks into other elements due to crosstalk")] = 0.0 """Fraction of the signal that leaks into other elements due to crosstalk""" @@ -115,10 +104,8 @@ def __post_init__(self): element.rescale(self.units) if self.sensitivity is None: self.sensitivity = 1.0 - self.sensitivity = normalize_sensitivity(self.sensitivity) - - def get_sensitivity(self, frequency: float) -> float: - return sensitivity_at_frequency(self.sensitivity, frequency) + elif isinstance(self.sensitivity, list): + self.sensitivity = [(float(f), float(v)) for f, v in self.sensitivity] def calc_output(self, cycles: float, frequency: float, dt: float, delays: np.ndarray = None, apod: np.ndarray = None, amplitude: float = 1.0) -> np.ndarray: if delays is None: @@ -126,10 +113,10 @@ def calc_output(self, cycles: float, frequency: float, dt: float, delays: np.nda if apod is None: apod = np.ones(self.numelements()) drive_signal = generate_drive_signal(cycles=cycles, frequency=frequency, dt=dt, amplitude=amplitude) - base_output = drive_signal * self.get_sensitivity(frequency) + base_output = drive_signal * sensitivity_at_frequency(self.sensitivity, frequency) outputs = [ np.concatenate( - [np.zeros(int(delay / dt)), a * element.get_sensitivity(frequency) * base_output], + [np.zeros(int(delay / dt)), a * sensitivity_at_frequency(element.sensitivity, frequency) * base_output], axis=0, ) for element, delay, a, in zip(self.elements, delays, apod) @@ -260,22 +247,20 @@ def merge(list_of_transducers:List[Transducer], offset_pins:bool=False, offset_i array_copies = [arr.copy() for arr in list_of_transducers] dict_key_sets = set() for array in array_copies: - array.sensitivity = normalize_sensitivity(array.sensitivity) - if isinstance(array.sensitivity, dict): - dict_key_sets.add(tuple(array.sensitivity["freq_Hz"])) + if isinstance(array.sensitivity, list): + dict_key_sets.add(tuple(f for f, _ in array.sensitivity)) for el in array.elements: - el.sensitivity = normalize_sensitivity(el.sensitivity) - if isinstance(el.sensitivity, dict): - dict_key_sets.add(tuple(el.sensitivity["freq_Hz"])) + if isinstance(el.sensitivity, list): + dict_key_sets.add(tuple(f for f, _ in el.sensitivity)) if len(dict_key_sets) > 1: raise ValueError("Cannot merge sensitivities with different frequency keys.") sensitivity_signatures = [] for array in array_copies: - if isinstance(array.sensitivity, dict): + if isinstance(array.sensitivity, list): sensitivity_signatures.append(( - tuple(array.sensitivity["freq_Hz"]), - tuple(array.sensitivity["values_Pa_per_V"]), + tuple(f for f, _ in array.sensitivity), + tuple(v for _, v in array.sensitivity), )) else: sensitivity_signatures.append(float(array.sensitivity)) diff --git a/src/openlifu/xdc/util.py b/src/openlifu/xdc/util.py index 21f84c84..e93d9e32 100644 --- a/src/openlifu/xdc/util.py +++ b/src/openlifu/xdc/util.py @@ -300,10 +300,7 @@ def report_to_matrix_dict(report_df: pd.DataFrame, focal_gain_lut=FOCAL_GAIN_LUT freq_df["Frequency"] = freq_df['Item'].apply(lambda x: float(re.search(r"(?<=^PNP \()\d+(?= kHz\)$)", x).group(0))) freq_df['focal_gain'] = freq_df['Frequency'].apply(lambda f: focal_gain_lut.interp(f0=f*1e3, crosstalk=matrix_dict['crosstalk_frac']).item()) freq_df['Sensitivity'] = freq_df['PNP'].astype(float)*1e6/freq_df['focal_gain']/voltage - matrix_dict['sensitivity'] = { - 'freq_Hz': [float(f * 1e3) for f in freq_df["Frequency"]], - 'values_Pa_per_V': [float(sens) for sens in freq_df['Sensitivity']], - } + matrix_dict['sensitivity'] = [(f*1e3, sens) for f, sens in zip(freq_df["Frequency"], freq_df['Sensitivity'])] matrix_dict['id'] = matrix_dict['id'].format(sn=sn.lower()) matrix_dict['name'] = matrix_dict['name'].format(sn=sn) return matrix_dict diff --git a/tests/test_transducer.py b/tests/test_transducer.py index 417b4ad2..5b34382f 100644 --- a/tests/test_transducer.py +++ b/tests/test_transducer.py @@ -121,7 +121,7 @@ def test_transducer_calc_output_interpolates_dictionary_sensitivity(): nx=1, ny=1, units="mm", - sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [1.0, 3.0]}, + sensitivity=[(100e3, 1.0), (300e3, 3.0)], ) transducer.elements[0].sensitivity = 1.0 cycles = 3 @@ -141,16 +141,6 @@ def test_transducer_calc_output_interpolates_dictionary_sensitivity(): np.testing.assert_allclose(output_low[0], expected_low) -def test_legacy_sensitivity_mapping_is_normalized_to_schema(): - transducer = Transducer( - sensitivity={"100000.0": 1.0, "300000.0": 3.0}, - ) - assert transducer.sensitivity == { - "freq_Hz": [100000.0, 300000.0], - "values_Pa_per_V": [1.0, 3.0], - } - - def test_element_calc_output_generates_signal_from_scalar_input(): element = Element(sensitivity=2.0) cycles = 4 @@ -184,13 +174,13 @@ def test_merge_pushes_transducer_sensitivity_into_elements(): nx=1, ny=1, units="mm", - sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [2.0, 4.0]}, + sensitivity=[(100e3, 2.0), (300e3, 4.0)], ) transducer_b = Transducer.gen_matrix_array( nx=1, ny=1, units="mm", - sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [3.0, 6.0]}, + sensitivity=[(100e3, 3.0), (300e3, 6.0)], ) transducer_a.elements[0].sensitivity = 5.0 transducer_b.elements[0].sensitivity = 7.0 @@ -198,8 +188,8 @@ def test_merge_pushes_transducer_sensitivity_into_elements(): merged = Transducer.merge([transducer_a, transducer_b], merge_mismatched_sensitivity=True) assert merged.sensitivity == 1.0 - assert merged.elements[0].sensitivity == {"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [10.0, 20.0]} - assert merged.elements[1].sensitivity == {"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [21.0, 42.0]} + assert merged.elements[0].sensitivity == [(100e3, 10.0),(300e3, 20.0)] + assert merged.elements[1].sensitivity == [(100e3, 21.0),(300e3, 42.0)] def test_merge_rejects_mismatched_sensitivity_keys(): @@ -207,16 +197,16 @@ def test_merge_rejects_mismatched_sensitivity_keys(): nx=1, ny=1, units="mm", - sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [2.0, 4.0]}, + sensitivity=[(100e3, 2.0), (300e3, 4.0)], ) transducer_b = Transducer.gen_matrix_array( nx=1, ny=1, units="mm", - sensitivity={"freq_Hz": [100e3, 400e3], "values_Pa_per_V": [3.0, 6.0]}, + sensitivity=[(100e3, 2.0), (300e3, 4.0)], ) - transducer_a.elements[0].sensitivity = {"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [5.0, 7.0]} - transducer_b.elements[0].sensitivity = {"freq_Hz": [100e3, 400e3], "values_Pa_per_V": [11.0, 13.0]} + transducer_a.elements[0].sensitivity = [(100e3, 5.0), (300e3, 7.0)] + transducer_b.elements[0].sensitivity = [(100e3, 11.0), (400e3, 13.0)] with pytest.raises(ValueError, match="different frequency keys"): Transducer.merge([transducer_a, transducer_b], merge_mismatched_sensitivity=True) @@ -307,9 +297,9 @@ def test_transducer_calc_output_combines_frequency_dependent_sensitivities(): nx=1, ny=1, units="mm", - sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [2.0, 4.0]}, + sensitivity=[(100e3, 2.0), (300e3, 4.0)], ) - transducer.elements[0].sensitivity = {"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [5.0, 9.0]} + transducer.elements[0].sensitivity = [(100e3, 5.0), (300e3, 9.0)] frequency = 200e3 dt = 1e-7 @@ -328,13 +318,13 @@ def test_transducer_array_to_transducer_preserves_frequency_dependent_sensitivit nx=1, ny=1, units="mm", - sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [2.0, 4.0]}, + sensitivity=[(100e3, 2.0), (300e3, 4.0)], ) transducer_b = Transducer.gen_matrix_array( nx=1, ny=1, units="mm", - sensitivity={"freq_Hz": [100e3, 300e3], "values_Pa_per_V": [1.0, 3.0]}, + sensitivity=[(100e3, 1.0), (300e3, 3.0)], ) transducer_a.elements[0].sensitivity = 5.0 transducer_b.elements[0].sensitivity = 7.0 @@ -360,3 +350,48 @@ def test_transducer_array_to_transducer_preserves_frequency_dependent_sensitivit np.testing.assert_allclose(output[0], 15.0 * expected_drive) np.testing.assert_allclose(output[1], 14.0 * expected_drive) + + +def test_element_sensitivity_from_json_is_list_of_tuples(): + """Sensitivity read from a JSON dict (list-of-lists) is converted to List[tuple[float, float]].""" + d = { + "index": 1, + "position": [0.0, 0.0, 0.0], + "orientation": [0.0, 0.0, 0.0], + "size": [1.0, 1.0], + "pin": 1, + "units": "mm", + "sensitivity": [[100e3, 1.0], [300e3, 3.0]], # JSON encodes tuples as lists + } + element = Element.from_dict(d) + assert isinstance(element.sensitivity, list) + assert all(isinstance(pair, tuple) for pair in element.sensitivity) + assert all(isinstance(f, float) and isinstance(v, float) for f, v in element.sensitivity) + assert element.sensitivity == [(100e3, 1.0), (300e3, 3.0)] + + +def test_transducer_sensitivity_from_json_is_list_of_tuples(): + """Transducer-level sensitivity survives a to_json/from_json round-trip as List[tuple[float, float]].""" + transducer = Transducer.gen_matrix_array( + nx=1, + ny=1, + units="mm", + sensitivity=[(100e3, 2.0), (300e3, 4.0)], + ) + reconstructed = Transducer.from_json(transducer.to_json()) + assert isinstance(reconstructed.sensitivity, list) + assert all(isinstance(pair, tuple) for pair in reconstructed.sensitivity) + assert all(isinstance(f, float) and isinstance(v, float) for f, v in reconstructed.sensitivity) + assert reconstructed.sensitivity == [(100e3, 2.0), (300e3, 4.0)] + + +def test_element_in_transducer_sensitivity_from_json_is_list_of_tuples(): + """Element-level sensitivity inside a Transducer survives a to_json/from_json round-trip as List[tuple[float, float]].""" + transducer = Transducer.gen_matrix_array(nx=1, ny=1, units="mm") + transducer.elements[0].sensitivity = [(100e3, 5.0), (300e3, 9.0)] + reconstructed = Transducer.from_json(transducer.to_json()) + el_sensitivity = reconstructed.elements[0].sensitivity + assert isinstance(el_sensitivity, list) + assert all(isinstance(pair, tuple) for pair in el_sensitivity) + assert all(isinstance(f, float) and isinstance(v, float) for f, v in el_sensitivity) + assert el_sensitivity == [(100e3, 5.0), (300e3, 9.0)]