Pulse shapes¶

A pulse is the acoustic pressure \(p_\text{ac}(t)\) that drives the bubble. Every pulse is a callable Equinox module: pulse(t) returns the pressure in pascals at time t in seconds. Because pulses are JAX PyTrees, you can compile, batch, and differentiate through their parameters, such as the amplitude, the frequency, or the weights of a neural network.

Every code block on this page runs as written, in order.

Tone bursts¶

A ToneBurst is a periodic carrier shape, multiplied by a peak pressure and gated by an envelope:

where \(P\) is pressure, \(s\) is shape, \(f\) is freq, \(\phi\) is phase, \(t_0\) is initial_time, \(N\) is cycle_num, and \(w\) is envelope.

import jax
import jax.numpy as jnp

from jbubble.pulse import HannEnvelope, ToneBurst
from jbubble.pulse.shapes import Sine

pulse = ToneBurst(freq=1e6, pressure=100e3, shape=Sine(), cycle_num=5)
hann = ToneBurst(
    freq=1e6, pressure=100e3, shape=Sine(), cycle_num=5, envelope=HannEnvelope()
)

ts = jnp.linspace(0.0, 6e-6, 601)  # [s]
p = jax.vmap(pulse)(ts)  # [Pa]
print(
    f"peak {p.max() / 1e3:.1f} kPa, Hann peak {jax.vmap(hann)(ts).max() / 1e3:.1f} kPa"
)

Carrier shapes¶

Shapes from jbubble.pulse.shapes have a nominal amplitude of 1.

from jbubble.pulse.shapes import Rectangular, TimeDomainSquare

square = ToneBurst(freq=1e6, pressure=100e3, shape=TimeDomainSquare(), cycle_num=5)
# Negative half-cycle pulses: -1 for 10 % of each period, 0 otherwise.
monopolar = ToneBurst(
    freq=1e6,
    pressure=100e3,
    shape=Rectangular(duty=0.1, high_level=-1.0, low_level=0.0),
    cycle_num=5,
)

Envelopes¶

The envelope gates the carrier to the active window. A single pulse, such as a tone burst, defaults to SoftRectangularEnvelope, and a sum of pulses to NoEnvelope.

To change the envelope of an existing pulse, call windowed:

from jbubble.pulse import TukeyEnvelope

tapered = pulse.windowed(TukeyEnvelope(alpha=0.2))

Timing: start, stop, and end¶

Each pulse has an active window, from t_start to t_stop. To delay a pulse, set its keyword-only initial_time. t_end, the stop time of a simulation that doesn't set t_max, is the start time plus twice the duration, so the bubble has time to ring down:

delayed = ToneBurst(
    freq=1e6, pressure=100e3, shape=Sine(), cycle_num=5, initial_time=2e-6
)
print(f"start {delayed.t_start * 1e6:.1f} µs, stop {delayed.t_stop * 1e6:.1f} µs, "
      f"end {delayed.t_end * 1e6:.1f} µs")

The solver steps to every edge in window_edges, so it can't step over a pulse that starts late.

Chirps¶

A ChirpPulse sweeps its frequency from freq_start to freq_end over sweep_duration. The sweep law is linear by default; ExponentialSweep sweeps geometrically, with equal time per octave:

from jbubble.pulse import ChirpPulse
from jbubble.pulse.chirp import ExponentialSweep

chirp = ChirpPulse(
    freq_start=0.5e6,  # [Hz]
    freq_end=3e6,  # [Hz]
    pressure=50e3,  # [Pa]
    sweep_duration=10e-6,  # [s]
    sweep=ExponentialSweep(),
    envelope=TukeyEnvelope(alpha=0.2),
)

Measured waveforms¶

A SampledPulse interpolates a sampled pressure trace, such as a hydrophone recording, linearly between the samples. The sample times are absolute: the active window runs from the first sample to the last. The following code builds a stand-in for a measured trace, sampled at 100 MHz:

from jbubble.pulse import RectangularEnvelope, SampledPulse

dt = 10e-9  # 100 MHz sampling [s]
t_rec = jnp.arange(800) * dt
measured = 80e3 * jnp.sin(2 * jnp.pi * 1e6 * t_rec) * jnp.exp(
    -(((t_rec - 4e-6) / 1.5e-6) ** 2)
)  # stand-in for a hydrophone recording [Pa]

recorded = SampledPulse.from_uniform(measured, dt=dt)
print(f"window: {recorded.t_start * 1e6:.2f} to {recorded.t_stop * 1e6:.2f} µs")

The default soft envelope halves the first and last samples. If the trace already starts and ends at zero, as here, that changes nothing. Otherwise, pass envelope=RectangularEnvelope() to keep them, or envelope=HannEnvelope() to taper a trace that's cut off mid-signal.

Neural pulses¶

A NeuralPulse lets a neural network define the waveform: the network maps the normalised time \((t - t_0)/T\), from 0 to 1 across the window, to the pressure in units of pressure_scale. Its weights are parameters that you can optimise, for example to design a drive that maximises a bubble response:

import equinox as eqx

from jbubble.pulse import NeuralPulse

net = eqx.nn.MLP(in_size=1, out_size=1, width_size=32, depth=2, key=jax.random.key(0))
learned = NeuralPulse(net=net, pulse_duration=5e-6, pressure_scale=100e3)

pulse_duration and pressure_scale are fixed configuration, not trainable parameters. To scale a neural pulse by a traced value, multiply it, as in the next section.

Combine pulses¶

Pulses support + and *:

low = ToneBurst(freq=1e6, pressure=80e3, shape=Sine(), cycle_num=10)
high = ToneBurst(freq=2e6, pressure=40e3, shape=Sine(), cycle_num=20)
dual = low + high  # Summed((low, high))
half = 0.5 * dual  # Scaled
smooth_dual = dual.windowed(HannEnvelope())  # windows the sum
print(type(dual).__name__, type(half).__name__, f"{float(dual.duration) * 1e6:.0f} µs")

Each part of a sum keeps its own envelope and timing, so a delayed part starts on time. A sum has no window of its own until you call windowed.

Simulate with a pulse and differentiate it¶

Any pulse drives any equation of motion. Because the pulse is part of the model, jax.grad differentiates a simulated quantity with respect to its parameters. The following code computes how the peak radius of a 2 µm lipid-coated bubble responds to the drive amplitude:

from jbubble import run_simulation
from jbubble.utils.presets import lipid_bubble

eom, _ = lipid_bubble()


def peak_ratio(k):
    result = run_simulation(eom, k * pulse)
    return result.radius.max() / eom.R0


value, slope = jax.value_and_grad(peak_ratio)(1.0)
print(f"peak R/R0 = {value:.3f}, d(peak R/R0)/dk = {slope:.3f}")

For the pulse shapes, the algebra, and simulations under each kind of pulse, with plots, see the example Driving pulses.