refactor(pm): add channel variability and NoNoise models

Project.toml

@@ -1,7 +1,7 @@
 name = "UnderwaterAcoustics"
 uuid = "0efb1f7a-1ce7-46d2-9f48-546a4c8fbb99"
 authors = ["Mandar Chitre <[email protected]>"]
-version = "0.3.3"
+version = "0.4.0"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

src/pm_basic.jl

@@ -7,9 +7,9 @@ export ReflectionCoef, FlatSurface, RayleighReflectionCoef, SurfaceLoss
 export Rock, Pebbles, SandyGravel, CoarseSand, MediumSand, FineSand, VeryFineSand
 export ClayeySand, CoarseSilt, SandySilt, Silt, FineSilt, SandyClay, SiltyClay, Clay
 export Vacuum, SeaState0, SeaState1, SeaState2, SeaState3, SeaState4
-export SeaState5, SeaState6, SeaState7, SeaState8, SeaState9
+export SeaState5, SeaState6, SeaState7, SeaState8, SeaState9, NoVariability
 export AcousticReceiverGrid2D, AcousticReceiverGrid3D, NarrowbandAcousticSource
-export WhiteGaussianNoise, RedGaussianNoise, Pinger, SampledAcousticSource
+export WhiteGaussianNoise, RedGaussianNoise, Pinger, SampledAcousticSource, NoNoise
 
 ### sound speed profiles
 
@@ -314,9 +314,23 @@ const SeaState7 = SurfaceLoss(26.5)
 const SeaState8 = SurfaceLoss(30.6)
 const SeaState9 = SurfaceLoss(32.9)
 
+### time-invariant impulse response model
+
+"""
+$(TYPEDEF)
+No time variability.
+"""
+struct NoVariability <: VariabilityModel end
+
+function impulseresponse(::NoVariability, arrivals::Vector{<:Arrival}, fs; abstime=true, ntaps=0)
+  ir = impulseresponse(model, arrivals, fs; abstime, ntaps)
+  n = floor(Int, T / Δt) + 1
+  repeat(ir; outer=(1,n))
+end
+
 ### basic environmental model
 
-Base.@kwdef struct BasicUnderwaterEnvironment{T1<:Altimetry, T2<:Bathymetry, T3<:SoundSpeedProfile, T4<:Number, T5<:ReflectionModel, T6<:ReflectionModel, T7} <: UnderwaterEnvironment
+Base.@kwdef struct BasicUnderwaterEnvironment{T1<:Altimetry, T2<:Bathymetry, T3<:SoundSpeedProfile, T4<:Number, T5<:ReflectionModel, T6<:ReflectionModel, T7<:NoiseModel, T8<:VariabilityModel} <: UnderwaterEnvironment
   altimetry::T1 = FlatSurface()
   bathymetry::T2 = ConstantDepth(20.0)
   ssp::T3 = IsoSSP(soundspeed())
@@ -325,6 +339,7 @@ Base.@kwdef struct BasicUnderwaterEnvironment{T1<:Altimetry, T2<:Bathymetry, T3<
   seasurface::T5 = Vacuum
   seabed::T6 = SandySilt
   noise::T7 = RedGaussianNoise(db2amp(120.0))
+  variability::T8 = NoVariability()
 end
 
 """
@@ -342,9 +357,27 @@ waterdensity(env::BasicUnderwaterEnvironment) = env.waterdensity
 seasurface(env::BasicUnderwaterEnvironment) = env.seasurface
 seabed(env::BasicUnderwaterEnvironment) = env.seabed
 noise(env::BasicUnderwaterEnvironment) = env.noise
+variability(env::BasicUnderwaterEnvironment) = env.variability
 
 ### noise models
 
+"""
+$(TYPEDEF)
+Ambient noise model representing no noise.
+
+---
+
+    NoNoise(σ)
+
+Create an ambient noise model representing no noise.
+"""
+struct NoNoise <: NoiseModel end
+
+function record(::NoNoise, duration, fs; start=0.0)
+  n = round(Int, duration*fs)
+  signal(zeros(ComplexF64, n), fs)
+end
+
 """
 $(TYPEDEF)
 Ambient noise model with Gaussian noise with a flat power spectral density.

src/pm_core.jl

@@ -10,7 +10,7 @@ export ReflectionModel, reflectioncoef
 export UnderwaterEnvironment, altimetry, bathymetry, ssp, salinity, seasurface, seabed, noise
 export AcousticSource, AcousticReceiver, location, nominalfrequency, phasor, record, recorder
 export PropagationModel, arrivals, transfercoef, transmissionloss, eigenrays, rays
-export NoiseModel
+export NoiseModel, variability
 export impulseresponse, channelmatrix
 
 ### interface: SoundSpeedProfile
@@ -130,6 +130,13 @@ Get the noise model for the underwater environment.
 """
 function noise end
 
+"""
+    variability(env::UnderwaterEnvironment)::VariabilityModel
+
+Get the variability model for the underwater environment.
+"""
+function variability end
+
 ### interface: AcousticSource
 
 abstract type AcousticSource end
@@ -161,7 +168,7 @@ function phasor end
     record(src::AcousticSource, duration, fs; start=0.0)
     record(noise::NoiseModel, duration, fs; start=0.0)
     record(model::PropagationModel, tx, rx, duration, fs; start=0.0)
-    record(model::PropagationModel, tx, rx, sig; reltime=true)
+    record(model::PropagationModel, tx, rx, sig; abstime=false)
 
 Make a recording of an acoustic source or ambient noise. The `start` time and
 `duration` are specified in seconds, and the recording is made at a sampling
@@ -177,7 +184,7 @@ irrespective of whether the source is real or complex.
 
 When a `sig` is specified, the sources are assumed to transmit the sampled signal in `sig`.
 The number of channels in `sig` must match the number of sources. The returned signal
-is the same type as the input signal (real or complex). If `reltime` is `true`, the
+is the same type as the input signal (real or complex). If `abstime` is `false`, the
 recorded signal starts at the first arrival, otherwise it starts at the beginning of the
 transmission.
 """
@@ -193,6 +200,10 @@ function location end
 abstract type NoiseModel end
 function record end
 
+### interface: VariabilityModel
+
+abstract type VariabilityModel end
+
 ### interface: PropagationModel
 
 abstract type PropagationModel{T<:UnderwaterEnvironment} end
@@ -362,16 +373,14 @@ function (rec::Recorder)(duration, fs; start=0.0)
       end
     end
   end
-  if rec.noisemodel !== missing && rec.noisemodel !== nothing
-    for k = 1:size(rec.arr, 2)
-      x[:,k] .+= record(rec.noisemodel, duration, fs; start)
-    end
+  for k = 1:size(rec.arr, 2)
+    x[:,k] .+= record(rec.noisemodel, duration, fs; start)
   end
   x̄ = rec.rx isa AbstractArray ? x : dropdims(x; dims=2)
   signal(x̄, fs)
 end
 
-function (rec::Recorder)(sig; fs=framerate(sig), reltime=true)
+function (rec::Recorder)(sig; fs=framerate(sig), abstime=false)
   nchannels(sig) == length(rec.tx) || throw(ArgumentError("Input signal must have $(length(rec.tx)) channel(s)"))
   mindelay, maxdelay = extrema(Iterators.flatten([[a1.time for a1 ∈ a] for a ∈ rec.arr]))
   n1 = round(Int, maxdelay * fs)
@@ -392,7 +401,7 @@ function (rec::Recorder)(sig; fs=framerate(sig), reltime=true)
       x[:,k] .+= record(rec.noisemodel, nsamples/fs, fs)
     end
   end
-  n3 = reltime ? round(Int, mindelay * fs) : 1
+  n3 = abstime ? 1 : round(Int, mindelay * fs)
   x̄ = rec.rx isa AbstractArray ? @view(x[n3:end,:]) : dropdims(@view x[n3:end,:]; dims=2)
   isanalytic(sig) ? signal(x̄, fs) : signal(convert.(eltype(sig), real.(x̄) .* √2), fs)
 end
@@ -408,12 +417,12 @@ If `approx` is `true`, a fast algorithm is used to generate a sparse channel mat
 that assigns an arrival to the nearest sampling time.
 """
 function channelmatrix(rec::Recorder, fs, ntaps=0; tx=1, rx=1, approx=false)
-  ir = impulseresponse(rec.arr[tx,rx], fs, ntaps; reltime=true, approx)
+  ir = impulseresponse(rec.arr[tx,rx], fs; approx, ntaps)
   TriangularToeplitz(ir, :L)
 end
 
 function channelmatrix(arrivals::Vector{<:Arrival}, fs, ntaps=0; approx=false)
-  ir = impulseresponse(arrivals, fs, ntaps; reltime=true, approx)
+  ir = impulseresponse(arrivals, fs; approx, ntaps)
   TriangularToeplitz(ir, :L)
 end
 
@@ -439,23 +448,25 @@ end
 
 # delegate recording task to an ephemeral recorder
 record(model::PropagationModel, tx, rx, duration, fs; start=0.0) = recorder(model, tx, rx)(duration, fs; start)
-record(model::PropagationModel, tx, rx, sig; reltime=true) = recorder(model, tx, rx)(sig; reltime)
+record(model::PropagationModel, tx, rx, sig; abstime=false) = recorder(model, tx, rx)(sig; abstime)
 
 """
-$(SIGNATURES)
+    impulseresponse(model::PropagationModel, arrivals::Vector{<:Arrival}, fs; abstime=false, approx=false, ntaps=0)
+    impulseresponse(arrivals::Vector{<:Arrival}, fs; abstime=false, approx=false, ntaps=0)
+
 Convert a vector of arrivals to a sampled impulse response time series at a
-sampling rate of `fs` Hz. If `ntaps` is zero, the number of taps of the impulse
-response are chosen automatically.
+sampling rate of `fs` Hz.
 
-If `reltime` is `true`, the impulse response start time is relative to the
+If `abstime` is `false`, the impulse response start time is relative to the
 first arrival, otherwise it is relative to the absolute time. If `approx`
 is `true`, a fast algorithm is used to generate a sparse impulse response
-that assigns an arrival to the nearest sampling time.
+that assigns an arrival to the nearest sampling time. If `ntaps` is zero,
+the number of taps of the impulse response are chosen automatically.
 """
-function impulseresponse(arrivals::Vector{<:Arrival}, fs, ntaps=0; reltime=false, approx=false)
+function impulseresponse(arrivals::Vector{<:Arrival}, fs; abstime=false, approx=false, ntaps=0)
   length(arrivals) == 0 && throw(ArgumentError("No arrivals"))
   mintime, maxtime = extrema(a.time for a ∈ arrivals)
-  reltime || (mintime = zero(typeof(arrivals[1].time)))
+  abstime && (mintime = zero(typeof(arrivals[1].time)))
   mintaps = ceil(Int, (maxtime-mintime) * fs) + 1
   if approx
     ntaps == 0 && (ntaps = mintaps)
@@ -488,6 +499,26 @@ function impulseresponse(arrivals::Vector{<:Arrival}, fs, ntaps=0; reltime=false
   ir
 end
 
+function impulseresponse(::PropagationModel, arrivals::Vector{<:Arrival}, fs; abstime=false, approx=false, ntaps=0)
+  impulseresponse(arrivals, fs; abstime, approx, ntaps)
+end
+
+"""
+    impulseresponse(model::PropagationModel, arrivals::Vector{<:Arrival}, fs, Δt, T; abstime=false, ntaps=0)
+
+Convert a vector of arrivals to a series of sampled time-varying impulse responses,
+each of which is sampled at a rate of `fs` Hz, with an impulse response generated
+every `Δt` seconds for `T` seconds. Returns a matrix of impulse responses, where
+each column corresponds to an impulse response at a specific time.
+
+If `abstime` is `false`, the impulse response start time is relative to the
+first arrival, otherwise it is relative to the absolute time. If `ntaps` is zero,
+the number of taps of the impulse response are chosen automatically.
+"""
+function impulseresponse(model::PropagationModel, arrivals::Vector{<:Arrival}, fs, Δt, T; abstime=false, ntaps=0)
+  impulseresponse(variability(model), arrivals, fs, Δt, T; abstime, ntaps)
+end
+
 # fast threaded map, assuming all entries have the same result type
 function tmap(f, x)
   x1 = first(x)

test/runtests.jl

@@ -290,8 +290,8 @@ end
   @test all([r[j].raypath[end][1] ≥ 100.0 for j ∈ 1:9])
   @test all([length(r[j].raypath) == r[j].surface + r[j].bottom + 2 for j ∈ 1:9])
 
-  ir1 = impulseresponse(arr, 10000.0; reltime=true, approx=true)
-  ir2 = impulseresponse(arr, 10000.0; reltime=false, approx=true)
+  ir1 = impulseresponse(arr, 10000.0; abstime=false, approx=true)
+  ir2 = impulseresponse(arr, 10000.0; abstime=true, approx=true)
   @test length(ir2) ≈ length(ir1) + round(Int, 10000.0 * arr[1].time) atol=1
   @test length(ir2) == round(Int, 10000.0 * arr[end].time) + 1
   @test sum(ir1 .!= 0.0) == 7
@@ -301,21 +301,21 @@ end
   ndx = findall(abs.(ir2) .> 0)
   @test (ndx .- 1) ./ 10000.0 ≈ [arr[j].time for j ∈ 1:7] atol=1e-4
 
-  ir1a = impulseresponse(arr, 10000.0; reltime=true)
-  ir2a = impulseresponse(arr, 10000.0; reltime=false)
+  ir1a = impulseresponse(arr, 10000.0; abstime=false)
+  ir2a = impulseresponse(arr, 10000.0; abstime=true)
   @test length(ir2a) ≈ length(ir1a) + round(Int, 10000.0 * arr[1].time) atol=1
   @test length(ir2a) ≥ length(ir2)
   @test sum(abs2.(ir1a))/sum(abs2.(ir1)) ≈ 1.0 atol=0.05
   @test sum(abs2.(ir2a))/sum(abs2.(ir2)) ≈ 1.0 atol=0.05
 
-  @test length(impulseresponse(arr, 10000.0, 256; reltime=true, approx=true)) == 256
-  @test length(impulseresponse(arr, 10000.0, 64; reltime=true, approx=true)) == 64
-  @test length(impulseresponse(arr, 10000.0, 256; reltime=true, approx=false)) == 256
-  @test length(impulseresponse(arr, 10000.0, 64; reltime=true, approx=false)) == 64
-  @test length(impulseresponse(arr, 10000.0, 1024; reltime=false, approx=true)) == 1024
-  @test length(impulseresponse(arr, 10000.0, 700; reltime=false, approx=true)) == 700
-  @test length(impulseresponse(arr, 10000.0, 1024; reltime=false, approx=false)) == 1024
-  @test length(impulseresponse(arr, 10000.0, 700; reltime=false, approx=false)) == 700
+  @test length(impulseresponse(arr, 10000.0; ntaps=256, abstime=false, approx=true)) == 256
+  @test length(impulseresponse(arr, 10000.0; ntaps=64, abstime=false, approx=true)) == 64
+  @test length(impulseresponse(arr, 10000.0; ntaps=256, abstime=false, approx=false)) == 256
+  @test length(impulseresponse(arr, 10000.0; ntaps=64, abstime=false, approx=false)) == 64
+  @test length(impulseresponse(arr, 10000.0; ntaps=1024, abstime=true, approx=true)) == 1024
+  @test length(impulseresponse(arr, 10000.0; ntaps=700, abstime=true, approx=true)) == 700
+  @test length(impulseresponse(arr, 10000.0; ntaps=1024, abstime=true, approx=false)) == 1024
+  @test length(impulseresponse(arr, 10000.0; ntaps=700, abstime=true, approx=false)) == 700
 
   env = UnderwaterEnvironment(ssp=IsoSSP(1500.0))
   pm = PekerisRayModel(env, 2)
@@ -408,7 +408,7 @@ end
   sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
   @test size(sig) == (44100,2)
 
-  env = UnderwaterEnvironment(noise=missing)
+  env = UnderwaterEnvironment(noise=NoNoise())
   pm = PekerisRayModel(env, 7)
   tx = Pinger(0.0, -5.0, 1000.0; interval=0.3)
   rx = AcousticReceiver(100.0, -10.0)
@@ -459,7 +459,7 @@ end
   @test sig isa AbstractVector
   @test eltype(sig) === Float64
   @test 44100 < length(sig) < 44700
-  sig = rec(zeros(44100); fs=44100.0, reltime=false)
+  sig = rec(zeros(44100); fs=44100.0, abstime=true)
   @test sig isa AbstractVector
   @test eltype(sig) === Float64
   @test 47000 < length(sig) < 47600
@@ -484,12 +484,12 @@ end
   @test sig isa AbstractMatrix
   @test size(sig, 2) == 2
 
-  env = UnderwaterEnvironment(noise=nothing)
+  env = UnderwaterEnvironment(noise=NoNoise())
   pm = PekerisRayModel(env)
   tx = AcousticSource(0.0, -5.0, 10000.0)
   rx = AcousticReceiver(100.0, -10.0)
   rec = recorder(pm, tx, rx)
-  ir = impulseresponse(arrivals(pm, tx, rx), 44100.0, 500; approx=true, reltime=true)
+  ir = impulseresponse(arrivals(pm, tx, rx), 44100.0; ntaps=500, approx=true, abstime=false)
   X1 = channelmatrix(rec, 44100.0, 500; approx=true)
   X2 = channelmatrix(arrivals(pm, tx, rx), 44100.0, 500; approx=true)
   @test X1 == X2