diff --git a/src/ISM.jl b/src/ISM.jl index e260344..2ec2bb5 100644 --- a/src/ISM.jl +++ b/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 - [ISM(TxRx(B*p.position+o, B*p.B, p.directivity), tx, room, config) for p in Pn] + rxs, origin, B = array.txrx, array.origin, 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 + """ """ 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 - rx::SVector{3, T}, # reveiver position - B::SMatrix{3, 3, T}, # receiver orientation + tx::SVector{3,T}, # transmitter position + rx::SVector{3,T}, # reveiver position + B::SMatrix{3,3,T}, # receiver orientation dp::AbstractDirectivityPattern, # Receiver directivity pattern - 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) + 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{<: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 - rx::SVector{3, T}, # reveiver position - B::SMatrix{3, 3, T}, # receiver orientation + tx::SVector{3,T}, # transmitter position + rx::SVector{3,T}, # reveiver position + B::SMatrix{3,3,T}, # receiver orientation dp::AbstractDirectivityPattern, # Receiver directivity pattern - 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) + 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{<: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_e = min(floor(Int, (τ+Wd/2)*fs)+1, Nh) # end + i_s = max(ceil(Int, (τ - Wd / 2) * fs) + 1, 1) # start + 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 - w = 0.5 * (1.0 + cos(2π*t/Wd)) # Hann window + t = (i - 1) / fs - τ # time signature + w = 0.5 * (1.0 + cos(2π * t / Wd)) # Hann window h[i] += w * A * sinc(fs * t) # sinc end end diff --git a/src/RoomAcoustics.jl b/src/RoomAcoustics.jl index 1cb752c..5431f21 100644 --- a/src/RoomAcoustics.jl +++ b/src/RoomAcoustics.jl @@ -15,6 +15,7 @@ include("utils.jl") include("ISM.jl") include("moving_sources.jl") +export WASN include("WASN.jl") end # module diff --git a/src/WASN.jl b/src/WASN.jl index fdf0c03..c4ed7cd 100644 --- a/src/WASN.jl +++ b/src/WASN.jl @@ -16,11 +16,7 @@ struct Node δ::Real end -function Node(rx, δ=0.0) - Node(rx, δ) -end - -struct Event +struct Event # NOTE: TxNode może będzie lepszą nazwa? 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)] - s = [[conv(h, events[j].signal) for h in hs[i, j]] for i in eachindex(nodes), j in eachindex(events)] + 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 @@ -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) diff --git a/src/types.jl b/src/types.jl index 5b53ec8..75675ea 100644 --- a/src/types.jl +++ b/src/types.jl @@ -8,8 +8,8 @@ abstract type AbstractRoom end struct RectangularRoom{T<:Real} <: AbstractRoom c::T - L::Tuple{T, T, T} - β::Tuple{T, T, T, T, T, T} + L::Tuple{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 - order::Tuple{I, I} # Order of reflection [low, high] +struct ISMConfig{T<:Real,I<:Integer,R<:AbstractRNG} <: AbstractRIRConfig + 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 diff --git a/src/utils.jl b/src/utils.jl index 34fd73f..f08823b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -5,11 +5,12 @@ export h2RT60, Sabine_RT60 """ function h2RT60(h::AbstractVector{<:Number}, Fs::Real) - cs = cumsum(reverse(h.^2)) - edc = 10*log10.(reverse(cs./cs[end])) # energy decay curve + dB(x) = 10log10(x) + 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