Propagation & channel modeling

Quickstart guide

Propagation modeling

Let’s get started:

using UnderwaterAcoustics
using Plots

Creating an environmental description

We typically start with an environmental description:

env = UnderwaterEnvironment()

UnderwaterEnvironment(
  bathymetry = 100.0, 
  altimetry = 0.0, 
  temperature = 27.0, 
  salinity = 35.0, 
  pH = 8.1, 
  soundspeed = 1538.9235842, 
  density = 1022.7198310217424, 
  seabed = RigidBoundary, 
  surface = PressureReleaseBoundary, 
)

If the defaults don’t suit our needs, we can customize the environment. For example, if we wanted an environment with 20 m water depth, sandy-clay seabed, and a smooth sound speed profile with 1500 m/s near the surface, 1490 m/s at 10 m depth, and 1520 m/s near the seabed, we could define:

env = UnderwaterEnvironment(
  bathymetry = 20.0,
  seabed = SandyClay,
  soundspeed = SampledField([1500, 1490, 1520]; z=0:-10:-20, interp=CubicSpline())
)

UnderwaterEnvironment(
  bathymetry = 20.0, 
  altimetry = 0.0, 
  temperature = 27.0, 
  salinity = 35.0, 
  pH = 8.1, 
  soundspeed = SampledField(z-varying, 3 samples), 
  density = 1022.7198310217424, 
  seabed = FluidBoundary(ρ=1173.38, c=1504.93, δ=0.00242), 
  surface = PressureReleaseBoundary, 
)

If we have Plots.jl installed, we can use plot recipes to plot the environment or the sound speed profile. For example:

plot(env.soundspeed)

Let us construct a range-independent iso-velocity environment with water depth of 20 m:

env = UnderwaterEnvironment(bathymetry=20, seabed=SandyClay)

UnderwaterEnvironment(
  bathymetry = 20, 
  altimetry = 0.0, 
  temperature = 27.0, 
  salinity = 35.0, 
  pH = 8.1, 
  soundspeed = 1538.9235842, 
  density = 1022.7198310217424, 
  seabed = FluidBoundary(ρ=1173.38, c=1504.93, δ=0.00242), 
  surface = PressureReleaseBoundary, 
)

Selecting a propagation model

Once we have an environment, we need to select a propagation model. Since we have a range-independent iso-velocity environment, we can use the PekerisRayTracer:

pm = PekerisRayTracer(env)

PekerisRayTracer(h=20, max_bounces=3)

Tip

We can get a list of all available propagation models by calling models().

Setting up transmitters and receivers

Next, we need a source and a receiver:

tx = AcousticSource(0.0, -5.0, 1000.0)

TX[1000.0 Hz, 0 dB](0.0, 0.0, -5.0)

rx = AcousticReceiver(100.0, -10.0)

RX(100.0, 0.0, -10.0)

For improved readability, positions can be specified as tuples or named tuples, and the API fully supports Unitful.jl:

tx = AcousticSource((x=0u"m", z=-5u"m"), 1u"kHz")

TX[1000 Hz, 0 dB](0.0, 0.0, -5.0)

rx = AcousticReceiver(100u"m", -10u"m")

RX(100.0, 0.0, -10.0)

Tip

2-tuples are interpreted as (x, z) and 3-tuples as (x, y, z). The coordinate system has x and y axis in the horizontal plane, and z axis pointing upwards, with the nominal water surface being at 0 m. This means that all z coordinates in water are negative. If units are not specified, they are assumed to be S.I. units (meters for distances, Hz for frequency, etc).

We just defined an omnidirectional 1 kHz transmitter tx at a depth of 5 m at the origin, and an omnidirectional receiver rx at a range of 100 m and a depth of 10 m.

Running the model

Now that we have an environment, a propagation model, a transmitter and a receiver, we can modeling. First, we ask for all ray arrivals (eigenrays) between the transmitter and receiver:

rays = arrivals(pm, tx, rx)

7 element Vector{RayArrival}:
∠ -2.9°  0⤒  0⤓ ∠  2.9° |  65.06 ms | -40.0 dB ϕ   0.0°↝
∠  8.5°  1⤒  0⤓ ∠  8.5° |  65.71 ms | -40.1 dB ϕ 180.0°↝
∠-14.0°  0⤒  1⤓ ∠-14.0° |  66.98 ms | -62.6 dB ϕ 170.9°↝
∠ 19.3°  1⤒  1⤓ ∠-19.3° |  68.85 ms | -74.3 dB ϕ -23.4°↝
∠-24.2°  1⤒  1⤓ ∠ 24.2° |  71.26 ms | -80.4 dB ϕ-145.6°↝
∠ 28.8°  2⤒  1⤓ ∠ 28.8° |  74.16 ms | -73.5 dB ϕ  10.8°↝
∠-33.0°  1⤒  2⤓ ∠-33.0° |  77.50 ms | -100.7 dB ϕ-167.2°↝

For each eigenray, this shows us the launch angle, number of surface bounces, number of bottom bounces, arrival angle, travel time, transmission loss along that ray, and phase change. The last “↝” symbol indicates that the complete ray path is also available. We can plot the ray paths:

plot(env; xlims=(-10,110))
plot!(tx)
plot!(rx)
plot!(rays)

The red star is the transmitter and the blue circle is the receiver. The stronger eigenrays are shown in blue, while the weaker ones are shown in red.

Often, we are interested in the arrival structure at a receiver. We generate an impulse response sampled at 48 kSa/s and plot it:

plot(impulse_response(pm, tx, rx, 48000))

The first arrival is at about 1 ms, as the impulse response is generated with timing relative to the first arrival with a guard time of 1 ms. If we wanted timings relative to the transmission time, we can set keyword argument abstime to true:

plot(impulse_response(pm, tx, rx, 48000; abstime=true))

We can also get the complex acoustic field or the transmission loss in dB:

acoustic_field(pm, tx, rx)

0.011324695913076596 + 0.013441491695419415im

transmission_loss(pm, tx, rx)

35.10150548468688

We can also pass in arrays of receivers, if we want to compute transmission loss at many locations simultaneously. Some models are able to compute transmission loss on a Cartesian grid very efficiently. This is useful to plot transmission loss as a function of space.

To define a 1000×200 Cartesian grid with 0.1 m spacing and compute the transmission loss over the grid:

rxs = AcousticReceiverGrid2D(1.0:0.1:100, -20:0.1:0)
x = transmission_loss(pm, tx, rxs)

991×201 Matrix{Float64}:
 19.1283  19.2328  19.5428  19.8492  …   8.82646  11.1099  16.4209  87.1226
 19.1329  19.2373  19.5471  19.8536      8.87496  11.1689  16.4855  87.1228
 19.1381  19.2424  19.552   19.8586      8.92702  11.2323  16.5549  87.1232
 19.1441  19.2483  19.5577  19.8643      8.9824   11.2997  16.6288  87.1235
 19.151   19.255   19.5641  19.8709      9.04085  11.3708  16.7068  87.1239
 19.1587  19.2626  19.5715  19.8784  …   9.10213  11.4454  16.7886  87.1243
 19.1674  19.2712  19.5799  19.8869      9.16605  11.5232  16.8739  87.1247
 19.1772  19.2809  19.5893  19.8965      9.23241  11.604   16.9625  87.1251
 19.1882  19.2918  19.6     19.9073      9.30109  11.6875  17.0541  87.1256
 19.2005  19.3039  19.6119  19.9194      9.37202  11.7737  17.1485  87.1261
  ⋮                                  ⋱                               ⋮
 37.3743  37.4386  37.4768  37.492      57.8842   61.3577  67.2494  98.1616
 37.4077  37.474   37.5142  37.5316     57.9055   61.3787  67.2687  98.1828
 37.442   37.5103  37.5526  37.5721     57.9268   61.4     67.2892  98.2041
 37.4772  37.5475  37.5919  37.6135  …  57.948    61.4216  67.3108  98.2254
 37.5133  37.5856  37.632   37.6557     57.9692   61.4434  67.3335  98.2468
 37.5502  37.6246  37.673   37.6988     57.9903   61.4655  67.3572  98.2682
 37.5881  37.6644  37.7149  37.7427     58.0112   61.4876  67.3819  98.2895
 37.6268  37.7052  37.7577  37.7875     58.0319   61.5099  67.4075  98.311
 37.6664  37.7468  37.8013  37.8331  …  58.0524   61.5322  67.4338  98.3324

plot(env; xlims=(0,100))
plot!(rxs, x)

Channel modeling

Simulating a transmission

Once we have a propagation model, we can setup a channel simulation. Let’s say we had a 10 kHz source transmitting at a source level of 170 dB re 1 µPa @ 1m and we wanted to run our simulations at 192 kSa/s:

tx = AcousticSource((x=0u"m", z=-5u"m"), 10u"kHz"; spl=170)
ch = channel(pm, tx, rx, 192000)

SampledPassbandChannel(1×1, 192000.0 Hz)

This channel accepts one input acoustic channel and yields one output acoustic channel. To transmit a signal through the channel, we create a windowed pulse using SignalAnalysis.jl and transmit() it through the channel ch:

using SignalAnalysis

x = cw(10000, 0.001, 192000; window=(tukey, 0.5)) |> real
y = transmit(ch, x)

plot(
  plot(x; xlims=(0,5)),
  plot(y; xlims=(0,5));
  layout=(2,1)
)

The received signal is scaled to be in µPa, assuming the source level of the transmitter was specified in dB re 1 µPa @ 1 m. We can clearly see the multipath arrivals in the received signal. Had we specified keyword argument abstime = true, we would have also seen the signal delayed by 65 ms.

Adding noise

We can also simulate channels with noise. For example, if we wanted red Gaussian noise with standard deviation σ = 0.5 Pa, we can specify that when generating the channel:

ch = channel(pm, tx, rx, 192000; noise=RedGaussianNoise(0.5e6))
y = transmit(ch, x)

plot(
  plot(x; xlims=(0,5)),
  plot(y; xlims=(0,5));
  layout=(2,1)
)

We can see the noise in the timeseries. We can see the 1/f² variation (red Gaussian noise) in the power spectral density superimposed on the 10 kHz peak from the transmitted pulse:

psd(y)