another dev commit

src/ISM.jl

@@ -1,37 +1,38 @@
 export ISM
 
 function ISM(
+    tx::TxRx,
     array::TxRxArray,
-    tx::TxRx,
     room::AbstractRoom,
     config::ISMConfig;
 )
-    Pn, o, B = array.p, array.origin, array.B
+    rxs, origin, B = array.txrx, array.origin, array.B
-    [ISM(TxRx(B*p.position+o, B*p.B, p.directivity), tx, room, config) for p in Pn]
+    l2g = rx -> TxRx(B * rx.position + origin, B * rx.B, rx.directivity)
+    [ISM(tx, rx |> l2g, room, config) for rx in rxs]
 end
 
 
+
 """
 
 """
 function ISM(
+    tx::TxRx,
     rxs::AbstractVector{<:TxRx},
-    tx::TxRx,
     room::AbstractRoom,
     config::ISMConfig;
 )
-	[ISM(rx, tx, room, config) for rx in rxs]
+    [ISM(rx, tx, room, config) for rx in rxs]
 end
 
 
 
-
 """
 
 """
 function ISM(
-    rx::TxRx,
     tx::TxRx,
+    rx::TxRx,
     room::RectangularRoom,
     config::ISMConfig;
 )
@@ -53,9 +54,9 @@ function ISM(
 
     if config.hp
         return AllenBerkley_highpass100(h, config.fs)
-        # Zmień to na funkcie operującą na zaalokowanym już h
+        # TODO: Zmień to na funkcie operującą na zaalokowanym już h
         # AllenBerkley_highpass100!(h, config.fs)
-        #return h
+        # return h
     else
         return h
     end
@@ -67,15 +68,15 @@ end
 
 """
 function ISM_RectangularRoom_core(
-    tx::SVector{3, T},                  # transmitter position
+    tx::SVector{3,T},                  # transmitter position
-    rx::SVector{3, T},                  # reveiver position
+    rx::SVector{3,T},                  # reveiver position
-    B::SMatrix{3, 3, T},                # receiver orientation
+    B::SMatrix{3,3,T},                # receiver orientation
     dp::AbstractDirectivityPattern,     # Receiver directivity pattern
-    L::Tuple{T, T, T},                  # room size (Lx, Ly, Lz)
+    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)
+    β::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{<:Int,<:Int},         # order of reflections; min max
     Nh::Integer,                        # h lenght in samples
     Wd::T,                              # Window width
     ISD::T,                             # Random displacement of image source
@@ -110,15 +111,15 @@ end
 """
 function ISM_RectangularRoom_core!(
     h::AbstractVector{<:T},
-    tx::SVector{3, T},                  # transmitter position
+    tx::SVector{3,T},                  # transmitter position
-    rx::SVector{3, T},                  # reveiver position
+    rx::SVector{3,T},                  # reveiver position
-    B::SMatrix{3, 3, T},                # receiver orientation
+    B::SMatrix{3,3,T},                # receiver orientation
     dp::AbstractDirectivityPattern,     # Receiver directivity pattern
-    L::Tuple{T, T, T},                  # room size (Lx, Ly, Lz)
+    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)
+    β::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{<:Int,<:Int},         # order of reflections; min max
     Wd::T,                              # Window width
     ISD::T,                             # Random displacement of image source
     lrng::AbstractRNG                   # random number generator
@@ -128,7 +129,7 @@ function ISM_RectangularRoom_core!(
     Nh = length(h)
 
     # Samples to distance coefficient [m]
-    Γ = c/fs
+    Γ = c / fs
 
     # Transform size of the room from meters to samples
     Lₛ = L ./ Γ
@@ -158,7 +159,7 @@ 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[i] - rx[i]
                 end
 
                 # Position of [randomized] image source for given permutation
@@ -174,30 +175,30 @@ function ISM_RectangularRoom_core!(
                 dist = norm(tx_isp)
 
                 # Propagation time between receiver and image source
-                τ = dist/c
+                τ = dist / c
 
-                if τ <= Nh/fs # Check if it still in transfer function range
+                if τ <= Nh / fs # Check if it still in transfer function range
 
                     # Compute value of reflection coefficients
-                    b .=  β .^ abs.((n-q, n, l-j, l, m-k, m))
+                    b .= β .^ abs.((n - q, n, l - j, l, m - k, m))
 
                     # Direction of Arrival of ray
-                    DoA .= tx_isp./dist
+                    DoA .= tx_isp ./ dist
 
                     # Compute receiver directivity gain
                     DG = directivity_pattern(SVector{3}(DoA), B, dp)
 
                     # Compute attenuation coefficient
-                    A = DG * prod(b)/(4π * dist)
+                    A = DG * prod(b) / (4π * dist)
 
                     # Compute range of samples in transfer function
-                    i_s = max(ceil(Int, (τ-Wd/2)*fs)+1, 1)  # start
+                    i_s = max(ceil(Int, (τ - Wd / 2) * fs) + 1, 1)  # start
-                    i_e = min(floor(Int, (τ+Wd/2)*fs)+1, Nh) # end
+                    i_e = min(floor(Int, (τ + Wd / 2) * fs) + 1, Nh) # end
 
                     # Insert yet another impulse into transfer function
                     for i ∈ i_s:i_e
-                        t = (i-1)/fs - τ # time signature
+                        t = (i - 1) / fs - τ # time signature
-                        w = 0.5 * (1.0 + cos(2π*t/Wd)) # Hann window
+                        w = 0.5 * (1.0 + cos(2π * t / Wd)) # Hann window
                         h[i] += w * A * sinc(fs * t) # sinc
                     end
                 end

src/RoomAcoustics.jl

@@ -15,6 +15,7 @@ include("utils.jl")
 include("ISM.jl")
 include("moving_sources.jl")
 
+export WASN
 include("WASN.jl")
 
 end # module

src/WASN.jl

@@ -16,11 +16,7 @@ struct Node
     δ::Real
 end
 
-function Node(rx, δ=0.0)
+struct Event # NOTE: TxNode może będzie lepszą nazwa?
-    Node(rx, δ)
-end
-
-struct Event
     tx::TxRx
     emission::Real
     fs::Real
@@ -33,8 +29,12 @@ function synth_events(
     room::AbstractRoom,
     rir_config::AbstractRIRConfig,
 )
-    hs = [ISM(nodes[i].rx, events[j].tx, room, rir_config) for i in eachindex(nodes), j in eachindex(events)]
+    hs = Matrix{Vector{Vector{Float64}}}(undef, nodes |> length, events |> length)
-    s = [[conv(h, events[j].signal) for h in hs[i, j]] for i in eachindex(nodes), j in eachindex(events)]
+    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
@@ -54,9 +54,10 @@ function synthesise(
     # 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.p |> length) for node ∈ nodes]
+    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)

src/types.jl

@@ -8,8 +8,8 @@ abstract type AbstractRoom end
 
 struct RectangularRoom{T<:Real} <: AbstractRoom
     c::T
-    L::Tuple{T, T, T}
+    L::Tuple{T,T,T}
-    β::Tuple{T, T, T, T, T, T}
+    β::Tuple{T,T,T,T,T,T}
 end
 
 
@@ -19,8 +19,8 @@ abstract type AbstractRIRConfig end
 """
 
 """
-struct ISMConfig{T<:Real, I<:Integer, R<:AbstractRNG} <: AbstractRIRConfig
+struct ISMConfig{T<:Real,I<:Integer,R<:AbstractRNG} <: AbstractRIRConfig
-    order::Tuple{I, I}  # Order of reflection [low, high]
+    order::Tuple{I,I}  # Order of reflection [low, high]
     fs::T               # Sampling frequency
     N::I                # Number of samples in impulse response
     Wd::T               # Single impulse width

src/utils.jl

@@ -5,11 +5,12 @@ export h2RT60, Sabine_RT60
 
 """
 function h2RT60(h::AbstractVector{<:Number}, Fs::Real)
-    cs = cumsum(reverse(h.^2))
+    dB(x) = 10log10(x)
-    edc = 10*log10.(reverse(cs./cs[end])) # energy decay curve
+    normalize(x) = x ./ x[begin]
+    edc = h .|> abs2 |> reverse |> cumsum |> reverse |> normalize .|> dB
 
     ind = findfirst(edc .<= -60. )
-    if ind == nothing 
+    if ind === nothing 
         rt = length(h)/Fs 
     else
         rt = ind/Fs