P/frequency dependent sensitivity

src/openlifu/plan/solution.py

@@ -256,8 +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, :])
-            p0_Pa = np.max(output_signal_Pa, axis=1)
+            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, :],
+                amplitude=self.pulse.amplitude * self.voltage,
+            ), axis=1)
 
             mainlobe_mask = get_mask(
                 pnp_MPa,

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
@@ -139,7 +140,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 = 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")
@@ -156,14 +158,14 @@ 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(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
+        # Simulate crosstalk by accumulating fraction of the source signal from each element into
+        # its neighbors within a certain distance. 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_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,7 @@ 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_mat[dst_idx,:] += arr.crosstalk_frac*source_mat[src_idx,:]
         source_mat = crosstalk_mat
     karray = get_karray(arr,
                         translation=array_offset,
@@ -219,6 +219,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

src/openlifu/xdc/element.py

@@ -2,14 +2,35 @@
 
 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
 
 
+def sensitivity_at_frequency(sensitivity: float | List[tuple[float, float]], frequency: float) -> float:
+    """Return sensitivity at `frequency`; list form assumes frequencies are sorted ascending."""
+    if isinstance(sensitivity, list):
+        freqs, values = zip(*sensitivity)
+        return float(np.interp(frequency, freqs, values))
+    return float(sensitivity)
+
+
+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.")
+    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))))
+    t = np.arange(n_samples, dtype=np.float64) * dt
+    return amplitude * np.sin(2 * np.pi * frequency * t)
+
+
 def matrix2xyz(matrix):
     x = matrix[0, 3]
     y = matrix[1, 3]
@@ -43,15 +64,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 | 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)"""
 
-    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 +83,10 @@ 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, list):
+            self.sensitivity = sorted(((float(f), float(v)) for f, v in self.sensitivity), key=lambda t: t[0])
 
     @property
     def x(self):
@@ -141,17 +152,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, 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)
@@ -225,17 +231,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 +268,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 +278,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)