From 6d93c1bcf37794bab2329867d0e47edbda9447cf Mon Sep 17 00:00:00 2001 From: Aaron Kaw Date: Thu, 15 Aug 2024 01:09:03 +1000 Subject: [PATCH] Really happy with sonar equation implementation prototype --- book/OceanSonar.typ | 234 ++++++++++++++++++++++--- prototyping/convenience_structure.md | 31 ---- prototyping/sonar_equations.jl | 80 +++++++++ prototyping/sonar_scenario.jl | 207 ---------------------- test/99_postliminary/abstract_trees.jl | 5 +- 5 files changed, 289 insertions(+), 268 deletions(-) delete mode 100644 prototyping/convenience_structure.md create mode 100644 prototyping/sonar_equations.jl delete mode 100644 prototyping/sonar_scenario.jl diff --git a/book/OceanSonar.typ b/book/OceanSonar.typ index 3d82e42..9255d72 100644 --- a/book/OceanSonar.typ +++ b/book/OceanSonar.typ @@ -1,85 +1,180 @@ +#import "@preview/unify:0.6.0": unit #import "@preview/jlyfish:0.1.0": * #set heading(numbering: "1.") #set math.equation(numbering: "(1)") +#show heading.where(level: 1): it => {pagebreak(weak: false);it} +#show heading.where(level: 2): it => {pagebreak(weak: false);it} + #read-julia-output(json("OceanSonar-jlyfish.json")) #jl-pkg( + "AbstractTrees", "CairoMakie", "Distributions", + "InteractiveUtils", "Typstry" ) -= Ocean Sonar -== Sonar Oceanography -== Ocean Acoustics -== Sonar Signal Processing +#align(center, text(17pt)[Ocean Sonar]) + +#outline(depth: 2, indent: auto) + +#pagebreak() + +== Pillars of Ocean Sonar + +== Mathematical Preliminaries + +=== Special Functions + +==== Error Function and Variations + +==== Standard Normal Cumulative Distribution Function + +== Programming Preliminaries + += Sonar Oceanography + += Ocean Acoustics -== Statistical Detection Theory += Sonar Signal Processing -=== Probability Metrics +== Sonar Scenarios -The vector of received acoustic data is modelled as $arrow(x)$. The two possibilities to be discerned upon receiving acoustic data are: +#jl(recompute: false, code: true, +``` +using AbstractTrees: print_tree +using InteractiveUtils: subtypes + +import AbstractTrees: children + +abstract type SonarType end + +abstract type Passive <: SonarType end +abstract type Active <: SonarType end + +abstract type Narrowband <: Passive end +abstract type Broadband <: Passive end +abstract type Intercept <: Passive end + +children(::Type{T}) where T <: SonarType = subtypes(T) + +print_tree(SonarType) +``` +) -- Noise only -- Signal present in the noise += Statistical Detection Theory -The hypothesis that the acoustic data contains only noise is labelled $H_0$, and the hypothesis that the acoustic data contains a signal with the noise is labelled $H_1$. +== Probability Metrics -The distributions of the acoustic data under the noise-only hypothesis $X_0 = T(arrow(n))$ and the signal-present hypothesis $X_1 = T(arrow(s) + arrow(n))$ depend on your noise $arrow(n)$ and signal $arrow(s)$ distributions and the choice of data processor $T$ for enhancing the signal-to-noise ratio. +The collection of received acoustic data is modelled in a vector as $arrow(x)$. +The units of such data may be in #unit("uPa") for amplitude data, +or #unit("uPa^2") for energy data. +(TODO: verify this.) -The probability of false alarm $P_f$ is the metric for an incorrect decision that a signal is present when the true scenario is that the signal is absent. It is thus defined as +The objective is to determine if the acoustic data contains noise only, or has a signal present. +We define hypotheses: + +- $H_0$: Noise only +- $H_1$: Signal present + +The noise data $arrow(n)$ and signal data $arrow(s)$ are assumed to combine additively to form the acoustic data. +That is, under each hypothesis: $ -P_f = upright(Pr){X_0 >= h | H_0} +H_0 &: arrow(x)_0 = arrow(n) \ +H_1 &: arrow(x)_1 = arrow(n) + arrow(s) $ -Similarly, the probability of detection $P_d$ is the metric for a correct decision that a signal is present. It is defined as +The _detector operator_ $T$ acts on the acoustic data +to reduce the dimension of the data to a scalar +for comparison to a threshold value $h$ termed the _detector threshold_. +The transformed data is termed the _detector statistic_. +Under each hypothesis, the transformed data $X_0$ and $X_1$ are defined: $ -P_f = upright(Pr){X_1 >= h | H_1} +H_0 &: X_0 = T(arrow(x)_0) = T(arrow(n)) \ +H_1 &: X_1 = T(arrow(x)_1) = T(arrow(n) + arrow(s)) $ -The threshold value $h$ is termed the _detector threshold_, to be distinguished from the detection threshold. +The acoustic data $arrow(x)$ is modelled as having a distribution dependent on the noise and signal distributions. +Upon assuming a distribution for the signal and noise, +the distribution of the transformed data $X_0$ and $X_1$ can be inferred. + +The probability metrics of interest for applied detection theory +are the probabilities of false alarm $P_f$ and detection $P_d$, defined: + +$ +P_f = upright(Pr){X_0 >= h | H_0} \ +P_d = upright(Pr){X_1 >= h | H_1} +$ +That is, the $P_f$ is the probability of a false positive, +and $P_d$ is the probability of a true positive. The values of both probabilities need to be under consideration when making a decision as to the presence of a signal. -==== Gaussian Noise with Gaussian Signal +More generally, the probability metrics above can be expressed in terms of the cumulative distribution function $F_i$ for their respective hypotheses distributions $H_i$ for $i = 0, 1$: + +$ +H_0: P_f = 1 - F_0(h) \ +H_1: P_d = 1 - F_1(h) +$ + +In some special cases, the mathematical formulae for the detector statistics can be derived theoretically. +In general, this is not possible. +Literature approaches such a limitation either by approximations +or by Monte Carlo simulations. +Some theoretical results are complicated enough to also motivate approximations for particular conditions. +Theoretical derivations, approximations, and Monte Carlo simulations are explored in the following sub-sections. + +=== Gaussian Noise with Deterministic Signal The following assumptions are made: -- Noise is white, bandpass, wide-sense-stationary, ergodic, Gaussian. -- Noise spectral density is $N_0$ for $|f| in (f_c - W/2, f_c + W/2)$, $0$ elsewhere. -- $arrow(n) tilde cal(C)cal(N)(arrow(0), lambda_0 arrow(I))$. -- Signal is deterministic: $arrow(s) = ?$ +- The noise is white, meaning TODO. +- The noise is bandpass, meaning TODO. +- The noise is wide-sense-stationary ergodic, meaning TODO. +- The noise is Gaussian-distributed, i.e. $arrow(n) tilde cal("CN")(arrow(0), lambda_0 arrow(I))$. +- The noise spectral density is $N_0$ for $|f| in (f_c - W/2, f_c + W/2)$, $0$ elsewhere. +- The signal is deterministic: $arrow(s) = ?$ + +Under such conditions, a theoretical result is tractable. -When the signal is known exactly and the noise is Gaussian, a coherent matched filter (CMF) is the optimal processor. +When the signal is known exactly and the noise is Gaussian, +a coherent matched filter (CMF) is the optimal processor. +The CMF detector operator is defined as $ T(arrow(x)) = Re(arrow(s)^H arrow(X)) $ -This results in Gaussia-distributed detector decision statistics. +This results in Gaussian-distributed detector decision statistics. $ -arrow(s) + arrow(n) &tilde cal(N)(mu_N + mu_S, sigma_N) \ +arrow(s) + arrow(n) &tilde cal("CN")(?, ?) \ X_0 = T(arrow(n)) &tilde cal(N)(mu_0, sigma_0) \ -X_1 = T(arrow(s) + arrow(n)) &tilde cal(N)(mu_1, sigma_1) +X_1 = T(arrow(s) + arrow(n)) &tilde cal(N)(mu_1, sigma_0) $ with $ -mu_0 = T(mu_N) -mu_1 = T(mu_S + mu_N) +mu_0 = ? \ +mu_1 = ? \ +sigma_0 = ? $ +where we note that the standard deviation of the noise and signal-plus-noise are equivalent. + Probability metrics are $ P_f &= 1 - F_0(h) = 1 - Phi((h - mu_0) / sigma_0) \ P_d &= 1 - F_1(h) = 1 - Phi((h - mu_1) / sigma_0) -$ +$ + +where $Phi$ is the standard normal CDF, defined here as TODO. #jl(recompute: false, code: true, ``` @@ -100,3 +195,86 @@ lines!(axis, X, X1prob) fig ``` ) + +=== Gaussian Noise with Gaussian Signal + +The same assumptions are made here as in @prob_gaussian_deterministic, +bar the signal distribution: + +- The signal is fluctuating: $arrow(s) tilde cal("CN")(?, ?)$ + +Under such conditions, a theoretical result is tractable as in @prob_gaussian_deterministic, +with the addition of a fluctuating index $upright("FI")$ defined: + +$ +upright("FI") = sigma_1 / sigma_0 +$ + +where $sigma_i$ are the standard deviations of the detector statistics $X_i$ for hypotheses $H_i$, $i = 0, 1$. + +== Detection Index and Receiver Operating Characteristic Curves + +The _detection index_ $d$ is a measure of detectability. +A higher value implies an easier true positive over false positive detection. + +By definition: + +$ +d = ((upright(E)(X_1) - upright(E)(X_0)))^2 / (upright("Var")(X_0)) +$ + +The _receiver operating characteristic_ (ROC) curves are contour plots of $d$ or sometimes $10log_(10)(d)$ with respect to the probabilities of false alarm $P_f$ and detection $P_d$. + +=== Gaussian Noise with Deterministic Signal + +Manipulating @eqn_prob_gaussian_determinstic to make the detector threshold $h$ the subject: + +$ +h = mu_0 + sigma_0 Phi^(-1)(1 - P_f) \ +h = mu_1 + sigma_0 Phi^(-1)(1 - P_d) +$ + +then equating, yields + +$ +(mu_1 - mu_0) / sigma_0 = Phi^(-1)(1 - P_d) - Phi^(-1)(1 - P_f) +$ + +Noting that + +$ +upright(E)(X_0) = mu_0 \ +upright(E)(X_1) = mu_1 \ +upright("Var")(X_0) = sigma_0 +$ + +we then have + +$ +sqrt(d) = Phi^(-1)(1 - P_d) - Phi^(-1)(1 - P_f) +$ + +#jl(recompute: false, code: true, +``` +using Distributions, CairoMakie + +N = Normal(0, 1) + +Pf = range(10^-4, 10^-1, 301) +Pd = range(0.2, 0.9, 201) +d = @. (quantile(N, 1 - Pd') - quantile(N, 1 - Pf))^2 +d_decibels = 10log10.(d) + +contour(Pf, Pd, d_decibels, labels = true) +``` +) + +== Processing Gain + +== Detection Threshold + +== Signal Excess + +== Transition Detection Probability + +== Figure of Merit diff --git a/prototyping/convenience_structure.md b/prototyping/convenience_structure.md deleted file mode 100644 index 631b232..0000000 --- a/prototyping/convenience_structure.md +++ /dev/null @@ -1,31 +0,0 @@ -# Convenience Structure - -* Dispatchable structure. -* Function fields can be `OceanSonar.ModellingFunctor <: ... <: Function` or not. -* Design plotting recipes to enable conveniences on the flexibility of whether a user defined a field as an `OceanSonar.ModellingFunctor`. - -```julia -abstract type Position3D end - -@kwdef mutable struct RectangularPosition3D <: Position3D - -end - -@kwdef mutable struct - -@kwdef mutable struct Sensor - -end - -@kwdef mutable struct SensorArray{ - N, - Position3DType <: AbstractPosition, - BeamPatternFunctionType <: Function, - DirectivityIndexFunctionType <: Function, -} - sensors::NTuple{N, sensor::@NamedTuple{N, Sensor}, position::Position3DType} - - beam_pattern::BeamPatternFunctionType - directivity_index::DirectivityIndexFunctionType -end -``` diff --git a/prototyping/sonar_equations.jl b/prototyping/sonar_equations.jl new file mode 100644 index 0000000..80dbcdc --- /dev/null +++ b/prototyping/sonar_equations.jl @@ -0,0 +1,80 @@ +abstract type EmissionType end + +struct NoiseOnly <: EmissionType end +struct Signaling <: EmissionType end + +abstract type AbstractEntity{ET <: EmissionType} end + +struct AbsentEntity <: AbstractEntity end + +struct Entity{ET <: EmissionType} <: AbstractEntity + SL::Float64 + NL::Float64 + pos::RectangularCoordinate{2} +end + +source_level(own::Entity{Signaling}, tgt::Entity{NoiseOnly}, fac::AbsentEntity = AbsentEntity()) = own.SL +source_level(own::Entity{NoiseOnly}, tgt::Entity{NoiseOnly}, fac::AbsentEntity = AbsentEntity()) = tgt.NL +source_level(own::Entity{NoiseOnly}, tgt::Entity{Signaling}, fac::AbsentEntity = AbsentEntity()) = tgt.SL +source_level(own::Entity{NoiseOnly}, tgt::Entity{NoiseOnly}, fac::Entity{Signaling} ) = fac.SL + +transmission_loss(own::Entity{Signaling}, tgt::Entity{NoiseOnly}, fac::AbsentEntity = AbsentEntity()) = 2 * prop.from_ownship +transmission_loss(own::Entity{NoiseOnly}, tgt::Entity, fac::AbsentEntity = AbsentEntity()) = prop.from_ownship +transmission_loss(own::Entity{NoiseOnly}, tgt::Entity{NoiseOnly}, fac::Entity{Signaling} ) = prop.from_ownship + prop.from_facilitator + +target_strength(own::Entity{Signaling}, tgt::Entity{NoiseOnly}, fac::AbsentEntity = AbsentEntity()) = tgt.TS +target_strength(own::Entity{NoiseOnly}, tgt::Entity, fac::AbsentEntity = AbsentEntity()) = 0.0 +target_strength(own::Entity{NoiseOnly}, tgt::Entity{NoiseOnly}, fac::Entity{Signaling} ) = tgt.TS + +reverberation_level(env::Environment, prop::Propagation, own::Entity{Signaling}, tgt::Entity{NoiseOnly}, fac::AbsentEntity = AbsentEntity()) = reverberation_level(env, prop, own) +reverberation_level(env::Environment, prop::Propagation, own::Entity{NoiseOnly}, tgt::Entity, fac::AbsentEntity = AbsentEntity()) = 0.0 +reverberation_level(env::Environment, prop::Propagation, own::Entity{NoiseOnly}, tgt::Entity{NoiseOnly}, fac::Entity{Signaling} ) = reverberation(env, prop, own, fac) + +@kwdef struct AcousticConfig + beam_model::ModelName = ModelName("Gaussian") +end + +@kwdef struct PropagationConfig + acoustic_config::AcousticConfig = AcousticConfig("BeamTracing") +end + +struct Propagation + to_ownship + from_facilitator + + Propagation(config::PropagationConfig, env::Environment, grid::Grid, own::Entity, fac::AbsentEntity = AbsentEntity()) + Propagation(config::PropagationConfig, env::Environment, grid::Grid, own::Entity, fac::Entity{Signaling}) +end + +struct SonarEquation + SL + TL + TS + NL + RL + AG + SNR + DT + SE + POD + + function SonarEquation(env::Environment, proc::Processing, grid::Grid, own::Entity, tgt::Entity, fac::AbstractEntity = AbsentEntity(); + prop_config::PropagationConfig = PropagationConfig() + ) + SL = source_level(own, tgt, fac) + prop = Propagation(prop_config, env, grid, own, fac) + TL = transmission_loss(prop, own, tgt, fac) + TS = target_strength(own, tgt, fac) + NL = env.NL ⊕ own.NL + RL = reverberation_level(env, prop, own, tgt, fac) + AG = array_gain(env, own) + DT = detection_threshold(own, tgt, fac) + + SNR = SL - TL + TS - ((NL - AG) ⊕ RL) + SE = SNR - DT + + POD = probability_of_detection(proc, own, tgt, fac, SE) + + new(SL, TL, TS, NL, RL, AG, SNR, DT, SE, POD) + end +end diff --git a/prototyping/sonar_scenario.jl b/prototyping/sonar_scenario.jl deleted file mode 100644 index 1fcf175..0000000 --- a/prototyping/sonar_scenario.jl +++ /dev/null @@ -1,207 +0,0 @@ -abstract type VerticalDirection end -abstract type Downward <: VerticalDirection end -abstract type Upward <: VerticalDirection end - -struct RectangularHorizontalCoordinate2D - x::Float64 - y::Float64 -end - -struct RectangularVerticalCoordinate2D{VD <: VerticalDirection} - r::Float64 - z::Float64 -end - -struct PolarHorizontalCoordinate2D - r::Float64 - θ::Float64 -end - -struct SphericalCoordinate2D{VD <: VerticalDirection} - θ::Float64 - φ::Float64 -end - -struct PolarVerticalCoordinate2D{VD <: VerticalDirection} - r::Float64 - φ::Float64 -end - -struct RectangularCoordinate3D{VD <: VerticalDirection} - x::Float64 - y::Float64 - z::Float64 -end - -struct CylindricalCoordinate3D{VD <: VerticalDirection} - r::Float64 - θ::Float64 - z::Float64 -end - -struct SphericalCoordinate3D{VD <: VerticalDirection} - r::Float64 - θ::Float64 - φ::Float64 -end - -function RectangularCoordinate(pos::CylindricalCoordinate3D{VD}) where {VD <: VerticalDirection} - RectangularCoordinate3D{VD}(pos.r * cos(pos.θ), pos.r * sin(pos.θ), pos.z) -end - -struct AcousticSensor{ResponseFunctionType <: Function} - name::String - response::ResponseFunctionType -end - -struct AcousticArray{SourceLevelFunctionType <: Function} - name::String - acoustic_sensors::Vector{@NamedTuple{position::RectangularCoordinate3D{Upward}, sensor::AcousticSensor, weight::Float64}} - source_level::SourceLevelFunctionType -end - -struct Entity{ - NoiseLevelFunctionType <: Function, - TargetStrengthFunctionType <: Function -} - name::String - noise_level::NoiseLevelFunctionType - target_strength::TargetStrengthFunctionType - arrays::Dict{String, @NamedTuple{position::CylindricalCoordinate3D{Upward}, orientation::SphericalCoordinate2D{Upward}}} -end - -struct EntityPlacement{ - NoiseLevelFunctionType <: Function, - TargetStrengthFunctionType <: Function -} - entity::Entity{NoiseLevelFunctionType, TargetStrengthFunctionType} - position::CylindricalCoordinate3D{Upward} - orientation::SphericalCoordinate2D{Upward} -end - -@kwdef mutable struct AcousticMedium{ - DensityFunctionType <: Function, - WindFunctionType <: Function, - SoundSpeedFunctionType <: Function, - ShearSoundSpeedFunctionType <: Function, - AttenuationFunctionType <: Function, - ShearAttenuationFunctionType <: Function -} - density::DensityFunctionType = @initialise_function - wind_speed::WindFunctionType = @initialise_function - sound_speed::SoundSpeedFunctionType = @initialise_function - shear_sound_speed::ShearSoundSpeedFunctionType = @initialise_function - attenuation::AttenuationFunctionType = @initialise_function - shear_attenuation::ShearAttenuationFunctionType = @initialise_function -end - -""" -TODO: Figure out how to implement an arbitrary number of layers whilst keeping type stability. -""" -struct Environment{ - AtmosphereHeightFunctionType <: Function, - SedimentDepthFunctionType <: Function, - BasementDepthFunctionType <: Function, -} - atmosphere::@NamedTuple{height::AtmosphereHeightFunctionType, medium::AcousticMedium} - ocean::AcousticMedium - sediment::@NamedTuple{depth::SedimentDepthFunctionType, medium::AcousticMedium} - basement::@NamedTuple{depth::BasementDepthFunctionType, medium::AcousticMedium} -end - -abstract type WaveForm end - -struct LFM <: WaveForm end -struct HFM <: WaveForm end - -abstract type SignalProcessing end - -struct NarrowbandProcessing <: SignalProcessing end -struct BroadbandProcessing <: SignalProcessing end -struct InterceptProcessing <: SignalProcessing end -struct ContinuousWaveProcessing <: SignalProcessing end -struct FrequencyModulatedProcessing{WF <: WaveForm} <: SignalProcessing end - -struct Propagation3D - function Propagation3D( - environment::Environment, - entity::EntityPlacement - ) - - end -end - -""" -TODO: -* Incorporate beampatterns of transmitters, target reflection, and receivers. -""" -function propagation(environment, entities, transmitter_names, target_name, receiver_names, x_values, y_values, z_values) - p = zeros(ComplexF64, ((x_values, y_values, z_values) .|> length)...) - - return @. -20log10(p |> abs) -end - -struct SonarPerformance - environment::Environment - - processing::SignalProcessing - - entities::Dict{String, EntityPlacement} - target_name::String - transmitter_names::String - receiver_names::Vector{String} - - x::AbstractVector{Float64} - y::AbstractVector{Float64} - z::AbstractVector{Float64} - - PL::Array{3, Float64} - - function SonarPerformance( - environment::Environment, - processing::SignalProcessing, - entities::Dict{String, EntityPlacement}, - x_values::AbstractVector{<:Real}, - y_values::AbstractVector{<:Real}, - z_values::AbstractVector{<:Real}, - receiver_names::AbstractVector{<:AbstractString}, - target_name::AbstractString, - transmitter_names::AbstractString = String[] - ) - @assert allunique(entity -> entity.name, entities) - - target = entities_dict[target_name].entity - - "Weighted-pressure ray trace from each transmitter and receiver." - propagations = Dict( - name => Beams(environment, entities[name]) - for name in [transmitter_names; receiver_names] - ) - - SL_PL_TS, NL, RL = propagation(environment, entities, transmitter_names, target_name, receiver_names, x_values, y_values, z_values) - - DI = directivity_index(environment, processing, entities, transmitter_names, target_name, receiver_names, x_values, y_values, z_values) - - DT = detection(environment, processing, entities, transmitter_names, target_name, receiver_names, x_values, y_values, z_values) - - SNR = @. SL_PL_TS - NL_RL_DI - SE = @. SNR - DT - - new( - environment, - processing, - entities, - target_name, - transmitter_names, - receiver_names, - x_values, - y_values, - z_values, - SL_PL_TS, - NL_RL_DI, - DT, - SNR, - SE - ) - end -end diff --git a/test/99_postliminary/abstract_trees.jl b/test/99_postliminary/abstract_trees.jl index 3ab70ad..e89f940 100644 --- a/test/99_postliminary/abstract_trees.jl +++ b/test/99_postliminary/abstract_trees.jl @@ -1,8 +1,9 @@ using OceanSonar using Test +# Just test they don't error. buf = IOBuffer() @test print_tree(buf, SonarType) |> isnothing @test print_tree(buf, OceanSonar.ModellingType) |> isnothing -@test print_tree(buf, OceanSonar.EnvironmentComponent) |> isnothing -@test !isempty(buf |> take! |> String) \ No newline at end of file +@test print_tree(buf, OceanSonar.SpatialDimensionSize) |> isnothing +@test !isempty(buf |> take! |> String)