source directivity for ISM rect core

README.md

@@ -4,7 +4,7 @@
 
 
 ```julia
-] add https://git.sr.ht/~zymon/RoomAcoustics.jl
+] add https://codeberg.org/zymon/RoomAcoustics.jl
 ```
 
 

src/ISM.jl

@@ -2,11 +2,11 @@ export ISM
 
 function ISM(
     tx::TxRx,
-    array::TxRxArray,
+    rx_array::TxRxArray,
     room::AbstractRoom,
     config::ISMConfig;
 )
-    rxs, origin, B = array.txrx, array.origin, array.B
+    rxs, origin, B = rx_array.txrx, rx_array.origin, rx_array.B
     l2g = rx -> TxRx(B * rx.position + origin, B * rx.B, rx.directivity)
     [ISM(tx, rx |> l2g, room, config) for rx in rxs]
 end
@@ -38,9 +38,11 @@ function ISM(
 )
     h = ISM_RectangularRoom_core(
         tx.position,     # transmitter position
+        tx.directivity,  # transmitter directivity pattern
+        tx.B,            # transmitter orientation
         rx.position,     # receiver position
-        rx.B,            # receiver orientation
         rx.directivity,  # receiver directivity pattern
+        rx.B,            # receiver orientation
         room.L,          # room size
         room.β,          # reflection coefficients
         room.c,          # velocity of the sound
@@ -68,31 +70,34 @@ end
 
 """
 function ISM_RectangularRoom_core(
-    tx::SVector{3,T},                  # transmitter position
+    tx_p::SVector{3,T},                 # transmitter position
-    rx::SVector{3,T},                  # reveiver position
+    tx_dp::AbstractDirectivityPattern,  # transmitter directivity pattern
-    B::SMatrix{3,3,T},                # receiver orientation
+    tx_B::SMatrix{3,3,T},               # receiver orientation
-    dp::AbstractDirectivityPattern,     # Receiver directivity pattern
+    rx_p::SVector{3,T},                 # reveiver position
-    L::Tuple{T,T,T},                  # room size (Lx, Ly, Lz)
+    rx_dp::AbstractDirectivityPattern,  # receiver directivity pattern
-    β::Tuple{T,T,T,T,T,T},         # Reflection coefficients (βx1, βx2, βy1, βy2, βz1, βz2)
+    rx_B::SMatrix{3,3,T},               # receiver orientation
+    L::Tuple{T,T,T},                    # room size (Lx, Ly, Lz)
+    β::Tuple{T,T,T,T,T,T},              # reflection coefficients (βx1, βx2, βy1, βy2, βz1, βz2)
     c::T,                               # velocity of the wave
     fs::T,                              # sampling frequeyncy
-    order::Tuple{<:Int,<:Int},         # order of reflections; min max
+    order::Tuple{I,I},                  # order of reflections; (min, max)
-    Nh::Integer,                        # h lenght in samples
+    Nh::Integer,                        # number of returned samples
-    Wd::T,                              # Window width
+    Wd::T,                              # window width
-    ISD::T,                             # Random displacement of image source
+    ISD::T,                             # random displacement of image source
     lrng::AbstractRNG                   # random number generator
-)::AbstractVector{T} where {T<:AbstractFloat}
+)::AbstractVector{T} where {T<:AbstractFloat, I<:Integer}
-
     # Allocate memory for the impulose response
     h = zeros(T, Nh)
 
     # Call
     ISM_RectangularRoom_core!(
-        h,      # Container for impulse response
+        h,      # container for impulse response
-        tx,     # transmitter position
+        tx_p,   # transmitter position
-        rx,     # receiver position
+        tx_dp,  # transmitter directivity pattern
-        B,      # receiver orientation
+        tx_B,   # transmitter orientation
-        dp,     # receiver directivity pattern
+        rx_p,   # receiver position
+        rx_dp,  # receiver directivity pattern
+        rx_B,   # receiver orientation
         L,      # room size
         β,      # reflection coefficients
         c,      # velocity of the sound
@@ -102,28 +107,34 @@ function ISM_RectangularRoom_core(
         ISD,    # random image source displacement
         lrng,   # Local random number generator
     )
-
+    h
-    return h
 end
 
 """
 
+References:
+[1] J. B. Allen and D. A. Berkley, “Image method for efficiently simulating small‐room acoustics,” The Journal of the Acoustical Society of America, vol. 65, no. 4, Art. no. 4, Apr. 1979, doi: 10.1121/1.382599.
+[2] P. M. Peterson, “Simulating the response of multiple microphones to a single acoustic source in a reverberant room,” The Journal of the Acoustical Society of America, vol. 80, no. 5, Art. no. 5, Nov. 1986, doi: 10.1121/1.394357.
+[3] E. De Sena, N. Antonello, M. Moonen, and T. van Waterschoot, “On the Modeling of Rectangular Geometries in Room Acoustic Simulations,” IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 23, no. 4, Art. no. 4, Apr. 2015, doi: 10.1109/TASLP.2015.2405476.
+[4] F. Brinkmann, V. Erbes, and S. Weinzierl, “Extending the closed form image source model for source directivity,” presented at the DAGA 2018, Munich, Germany, Mar. 2018.
 """
 function ISM_RectangularRoom_core!(
     h::AbstractVector{<:T},
-    tx::SVector{3,T},                  # transmitter position
+    tx_p::SVector{3,T},                 # transmitter position
-    rx::SVector{3,T},                  # reveiver position
+    tx_dp::AbstractDirectivityPattern,  # transmitter directivity pattern
-    B::SMatrix{3,3,T},                # receiver orientation
+    tx_B::SMatrix{3,3,T},               # receiver orientation
-    dp::AbstractDirectivityPattern,     # Receiver directivity pattern
+    rx_p::SVector{3,T},                 # reveiver position
-    L::Tuple{T,T,T},                  # room size (Lx, Ly, Lz)
+    rx_dp::AbstractDirectivityPattern,  # receiver directivity pattern
-    β::Tuple{T,T,T,T,T,T},         # Reflection coefficients (βx1, βx2, βy1, βy2, βz1, βz2)
+    rx_B::SMatrix{3,3,T},               # receiver orientation
+    L::Tuple{T,T,T},                    # room size (Lx, Ly, Lz)
+    β::Tuple{T,T,T,T,T,T},              # Reflection coefficients (βx1, βx2, βy1, βy2, βz1, βz2)
     c::T,                               # velocity of the wave
     fs::T,                              # sampling frequeyncy
-    order::Tuple{<:Int,<:Int},         # order of reflections; min max
+    order::Tuple{I,I},                  # order of reflections; min max
     Wd::T,                              # Window width
     ISD::T,                             # Random displacement of image source
-    lrng::AbstractRNG                   # random number generator
+    lrng::AbstractRNG;                  # random number generator
-)::AbstractVector{T} where {T<:AbstractFloat}
+) where {T<:AbstractFloat, I<:Integer}
 
     # Number of samples in impulose response
     Nh = length(h)
@@ -140,11 +151,12 @@ function ISM_RectangularRoom_core!(
     o_min, o_max = order
 
     # Allocate memory
-    rd = @MVector zeros(T, 3)       # Container for random displacement
+    rd = @MVector zeros(T, 3)     # Container for random displacement
-    tx_isp = @MVector zeros(T, 3)   # Container for relative image source position
+    isp = @MVector zeros(T, 3)    # Container for relative image source position
-    b = @MVector zeros(T, 6)        # Container for effective reflection coefficient
+    b = @MVector zeros(T, 6)      # Container for effective reflection coefficient
-    DoA = @MVector zeros(T, 3)      # Container for direction of incoming ray to receiver
+    rx_DoA = @MVector zeros(T, 3) # Container for direction of incoming ray to receiver
-    Rp = @MVector zeros(T, 3)       # 
+    tx_DoA = @MVector zeros(T, 3) # Container for direction of ray coming out from transmitter
+    Rp = @MVector zeros(T, 3)     #
 
     # Main loop
     for n = -N[1]:N[1], l = -N[2]:N[2], m = -N[3]:N[3]
@@ -159,37 +171,44 @@ function ISM_RectangularRoom_core!(
             if o_min <= o && (o <= o_max || o_max == -1)
                 # Compute Rp part
                 for i = 1:3
-                    Rp[i] = (1 .- 2 * p[i]) * tx[i] - rx[i]
+                    Rp[i] = (1 .- 2 * p[i]) * tx_p[i] - rx_p[i]
                 end
 
-                # Position of [randomized] image source for given permutation
+                # image source position for given permutation
-                tx_isp .= Rp .+ Rr
+                isp .= Rp .+ Rr
 
                 if ISD > 0.0 && o > 0
-                    # Generate random displacement for the image source
+                    # generate random displacement for the image source
                     randn!(lrng, rd)
-                    tx_isp .+= rd .* ISD
+                    isp .+= rd .* ISD
                 end
 
                 # Distance between receiver and image source
-                dist = norm(tx_isp)
+                d = norm(isp)
 
                 # Propagation time between receiver and image source
-                τ = dist / c
+                τ = d / c
 
-                if τ <= Nh / fs # Check if it still in transfer function range
+                if τ <= Nh / fs # Check if it still in transfer function rangΩ
 
                     # Compute value of reflection coefficients
                     b .= β .^ abs.((n - q, n, l - j, l, m - k, m))
 
-                    # Direction of Arrival of ray
+                    # Direction of Arrival of ray incoming from image source to receiver
-                    DoA .= tx_isp ./ dist
+                    rx_DoA .= isp ./ d
 
                     # Compute receiver directivity gain
-                    DG = directivity_pattern(SVector{3}(DoA), B, dp)
+                    rx_DG = directivity_pattern(SVector{3}(rx_DoA), rx_B, rx_dp)
+
+                    # Direction of Arrival of ray coming out from transmitter to wall
+                    perm = (abs(n - q) + abs(n), abs(l - j) + abs(l), abs(m - k) + abs(m))
+                    tx_DoA .= -rx_DoA .* (-1, -1, -1).^perm
+
+                    # Compute transmitter directivity gain
+                    tx_DG = directivity_pattern(SVector{3}(tx_DoA), tx_B, tx_dp)
 
                     # Compute attenuation coefficient
-                    A = DG * prod(b) / (4π * dist)
+                    A = tx_DG * rx_DG * prod(b) / (4π * d)
 
                     # Compute range of samples in transfer function
                     i_s = max(ceil(Int, (τ - Wd / 2) * fs) + 1, 1)  # start
@@ -206,5 +225,4 @@ function ISM_RectangularRoom_core!(
             end
         end
     end
-    h
 end

src/RoomAcoustics.jl

@@ -3,10 +3,12 @@ module RoomAcoustics
 using LinearAlgebra
 using StaticArrays
 using Statistics
-using DSP
 using Random
 using Random: GLOBAL_RNG
 
+using LoopVectorization
+using DSP: conv
+
 using TxRxModels
 using TxRxModels: AbstractDirectivityPattern
 
@@ -15,7 +17,6 @@ include("utils.jl")
 include("ISM.jl")
 include("moving_sources.jl")
 
-export WASN
+include("node_event.jl")
-include("WASN.jl")
 
 end # module

src/WASN.jl

@@ -1,77 +0,0 @@
-module WASN
-
-using LinearAlgebra
-using DSP: conv
-
-using TxRxModels
-
-using ..RoomAcoustics: AbstractRoom, AbstractRIRConfig, ISM
-
-export Node, Event
-export synth_events
-
-struct Node
-    rx::TxRxArray
-    fs::Real
-    δ::Real
-end
-
-struct Event # NOTE: TxNode może będzie lepszą nazwa?
-    tx::TxRx
-    emission::Real
-    fs::Real
-    signal::AbstractVector
-end
-
-function synth_events(
-    nodes::AbstractVector{<:Node},
-    events::AbstractVector{<:Event},
-    room::AbstractRoom,
-    rir_config::AbstractRIRConfig,
-)
-    hs = Matrix{Vector{Vector{Float64}}}(undef, nodes |> length, events |> length)
-    iter = Iterators.product(nodes |> eachindex, events |> eachindex) |> collect
-    Threads.@threads for (i,j) ∈ iter
-        hs[i,j] = ISM(events[j].tx, nodes[i].rx, room, rir_config)
-    end
-    s = [[conv(h, events[j].signal) for h in hs[i, j]] for i ∈ eachindex(nodes), j ∈ eachindex(events)]
-
-    (signals=s, hs=hs)
-end
-
-function synthesise(
-    nodes::AbstractVector{<:Node},
-    events::AbstractVector{<:Event},
-    room::AbstractRoom,
-    rir_config::AbstractRIRConfig,
-)
-    # WARN: THIS FUNCTION ASSUMES NOW THAT ALL sampling rates are the same.
-    fs = first(nodes).fs
-
-    # Synthesise individual events for given nodes
-    e_signal, h_s = synth_events(nodes, events, room, rir_config)
-
-    # Find length of the output signal
-    δ_max = [node.δ for node ∈ nodes] |> maximum;
-    N = [ceil(Int, (events[j].emission + δ_max)*fs) + length(e_signal[1, j][1]) for j ∈ eachindex(events)] |> maximum
-    N += floor(Int, fs)
-
-    # Allocate memory for output signals
-    output = [zeros(N, node.rx.txrx |> length) for node ∈ nodes]
-
-    for i ∈ eachindex(nodes), j ∈ eachindex(events)
-        shift_n = floor(Int, (events[j].emission + nodes[i].δ)*fs)
-        event = e_signal[i,j]
-        for (idx, channel) in enumerate(event)
-            output[i][shift_n:shift_n+length(channel)-1, idx] .= channel
-        end
-    end
-
-    return (
-        node   = output,
-        event  = e_signal,
-        h      = h_s,
-    )
-end
-
-end # module WASN

src/node_event.jl

@@ -0,0 +1,57 @@
+
+export synthesise_events4nodes, synthesise
+
+function synthesise_events4nodes(
+    nodes::AbstractVector{<:Node},
+    events::AbstractVector{<:Event},
+    room::AbstractRoom,
+    rir_config::AbstractRIRConfig,
+)
+    hs = Matrix{Vector{Vector{Float64}}}(undef, nodes |> length, events |> length)
+    s = Matrix{Vector{Vector{Float64}}}(undef, nodes |> length, events |> length)
+    iter = Iterators.product(nodes |> eachindex, events |> eachindex) |> collect
+    for (i,j) ∈ iter
+        hs[i,j] = ISM(events[j].tx, nodes[i].rx, room, rir_config)
+        s[i,j] = [conv(h, events[j].signal) for h in hs[i, j]]
+    end
+    return (
+        signals=s,
+        hs=hs
+    )
+end
+
+function synthesise(
+    nodes::AbstractVector{<:Node},
+    events::AbstractVector{<:Event},
+    room::AbstractRoom,
+    rir_config::AbstractRIRConfig;
+    pad::Real = 0.0
+)
+    # WARN: AT THIS MOMENT THIS FUNCTION ASSUMES THAT ALL sampling rates are the same.
+    fs = first(nodes).fs
+
+    # Synthesise individual events for ech node
+    nodes_events, h = synthesise_events4nodes(nodes, events, room, rir_config)
+
+    # Compute length of the output signal
+    δ_max = [node.δ for node ∈ nodes] |> maximum;
+    N = [ceil(Int, (events[j].emission + δ_max)*fs) + length(nodes_events[1, j][1]) for j ∈ eachindex(events)] |> maximum
+    N += floor(Int, pad*fs)
+
+    # Allocate memory for output signals
+    output = [zeros(N, node.rx.txrx |> length) for node ∈ nodes]
+
+    for i ∈ eachindex(nodes), j ∈ eachindex(events)
+        shift_n = floor(Int, (events[j].emission + nodes[i].δ)*fs)+1
+        event = nodes_events[i,j]
+        for (idx, channel) in enumerate(event)
+            output[i][shift_n:shift_n+length(channel)-1, idx] .+= channel
+        end
+    end
+
+    return (
+        node   = output,
+        event  = nodes_events,
+        h      = h,
+    )
+end

src/types.jl

@@ -1,7 +1,7 @@
 
 export AbstractRoom, RectangularRoom
 export AbstractRIRConfig, ISMConfig
-
+export Node, Event
 
 
 abstract type AbstractRoom end
@@ -42,3 +42,25 @@ function ISMConfig(
 end
 
 
+struct Node{T<:Real}
+    rx::TxRxArray{T}
+    fs::T
+    δ::T
+    function Node(rx, fs::T, δ::T = 0.0) where T
+        δ < 0.0 && error("δ < 0")
+        fs < 0.0 && error("fs < 0")
+        new{T}(rx, fs, δ)
+    end
+end
+
+struct Event{T<:Real, V<:AbstractVector{<:T}} # NOTE: TxNode może będzie lepszą nazwa?
+    tx::TxRx{T}
+    emission::T
+    fs::T
+    signal::V
+    function Event(tx, emission::T, fs::T, signal::V) where {T, V}
+        emission < 0.0 && error("emission < 0")
+        fs < 0.0 && error("fs < 0")
+        new{T, V}(tx, emission, fs, signal)
+    end
+end