From 8d22ab308d901ad63eceb3947d6b8ac7fb392edc Mon Sep 17 00:00:00 2001 From: zymon Date: Sun, 2 Jun 2024 19:00:02 +0200 Subject: [PATCH] ISM refactor --- src/ISM.jl | 126 ++++++++++++--------------------------------------- src/types.jl | 13 ------ 2 files changed, 28 insertions(+), 111 deletions(-) diff --git a/src/ISM.jl b/src/ISM.jl index 5a2a982..f27eac8 100644 --- a/src/ISM.jl +++ b/src/ISM.jl @@ -36,26 +36,12 @@ function ISM( room::RectangularRoom, config::ISMConfig; ) - h = ISM_RectangularRoom_core( - tx.position, # transmitter position - tx.directivity, # transmitter directivity pattern - tx.B, # transmitter orientation - rx.position, # receiver position - rx.directivity, # receiver directivity pattern - rx.B, # receiver orientation - room.L, # room size - room.β, # reflection coefficients - room.c, # velocity of the sound - config.fs, # sampling frequency - config.order, # order of reflections - config.N, # length of response - config.Wd, # single impulse width - config.isd, # random image source displacement - config.lrng, # Local random number generator - ) + h = config.N + ISM!(h, tx, rx, room, config) if config.hp return AllenBerkley_highpass100(h, config.fs) + # TODO: Tego nie ma w ISM! # TODO: Zmień to na funkcie operującą na zaalokowanym już h # AllenBerkley_highpass100!(h, config.fs) # return h @@ -64,52 +50,6 @@ function ISM( end end - - -""" - -""" -function ISM_RectangularRoom_core( - tx_p::SVector{3,T}, # transmitter position - tx_dp::AbstractDirectivityPattern, # transmitter directivity pattern - tx_B::SMatrix{3,3,T}, # receiver orientation - rx_p::SVector{3,T}, # reveiver position - rx_dp::AbstractDirectivityPattern, # receiver directivity pattern - 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{I,I}, # order of reflections; (min, max) - Nh::Integer, # number of returned samples - Wd::T, # window width - ISD::T, # random displacement of image source - lrng::AbstractRNG # random number generator -)::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 - tx_p, # transmitter position - tx_dp, # transmitter directivity pattern - tx_B, # transmitter orientation - rx_p, # receiver position - rx_dp, # receiver directivity pattern - rx_B, # receiver orientation - L, # room size - β, # reflection coefficients - c, # velocity of the sound - fs, # sampling frequency - order, # order of reflections - Wd, # single impulse width - ISD, # random image source displacement - lrng, # Local random number generator - ) - h -end - """ References: @@ -118,42 +58,32 @@ References: [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}, # container for impulse reponse - tx_p::SVector{3,T}, # transmitter position - tx_dp::AbstractDirectivityPattern, # transmitter directivity pattern - tx_B::SMatrix{3,3,T}, # transmitter orientation - rx_p::SVector{3,T}, # reveiver position - rx_dp::AbstractDirectivityPattern, # receiver directivity pattern - 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{I,I}, # order of reflections; min max - Wd::T, # Window width - ISD::T, # Random displacement of image source position - lrng::AbstractRNG; # random number generator -) where {T<:AbstractFloat, I<:Integer} +function ISM!( + h::AbstractVector{T}, + tx::TxRx, + rx::TxRx, + room::RectangularRoom, + config::ISMConfig; +) where {T<:AbstractFloat} # Number of samples in impulse response Nh = length(h) # Samples to distance coefficient [m] - Γ = c / fs + Γ = room.c / config.fs # Transform size of the room from meters to samples - Lₛ = L ./ Γ + Lₛ = room.L ./ Γ # Impulse parameters - Δt = 1 / fs - ω = 2π / Wd - twid = π * fs + Δt = 1 / config.fs + ω = 2π / config.Wd + twid = π * config.fs # Compute maximal wall reflection N = ceil.(Int, Nh ./ (2 .* Lₛ)) - o_min, o_max = order + o_min, o_max = config.order # Allocate memory rd = @MVector zeros(T, 3) # Container for random displacement @@ -166,7 +96,7 @@ function ISM_RectangularRoom_core!( # Main loop for n = -N[1]:N[1], l = -N[2]:N[2], m = -N[3]:N[3] r = (n, l, m) # Wall reflection indicator - Rr = 2 .* r .* L # Wall lattice + Rr = 2 .* r .* room.L # Wall lattice for q ∈ 0:1, j ∈ 0:1, k ∈ 0:1 p = @SVector [q, j, k] # Permutation tuple @@ -176,48 +106,48 @@ 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_p[i] - rx_p[i] + Rp[i] = (1 .- 2 * p[i]) * tx.position[i] - rx.position[i] end # image source position for given permutation isp .= Rp .+ Rr - if ISD > 0.0 && o > 0 + if config.isd > 0.0 && o > 0 # generate random displacement for the image source - randn!(lrng, rd) - isp .+= rd .* ISD + randn!(config.lrng, rd) + isp .+= rd .* config.isd end # Distance between receiver and image source d = norm(isp) # Propagation time between receiver and image source - τ = d / c + τ = d / room.c - if τ <= Nh / fs # Check if it still in transfer function rangΩ + if τ <= Nh / config.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)) + b .= room.β .^ abs.((n - q, n, l - j, l, m - k, m)) # Direction of Arrival of ray incoming from image source to receiver rx_DoA .= isp ./ d # Compute receiver directivity gain - rx_DG = directivity_pattern(SVector{3}(rx_DoA), rx_B, rx_dp) + rx_DG = directivity_pattern(SVector{3}(rx_DoA), rx.B, rx.directivity) # 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) + tx_DG = directivity_pattern(SVector{3}(tx_DoA), tx.B, tx.directivity) # Compute attenuation coefficient 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 - i_e = min(floor(Int, (τ + Wd / 2) * fs) + 1, Nh) # end + i_s = max(ceil(Int, (τ - config.Wd / 2) * config.fs) + 1, 1) # start + i_e = min(floor(Int, (τ + config.Wd / 2) * config.fs) + 1, Nh) # end # Insert yet another impulse into transfer function for i ∈ i_s:i_e diff --git a/src/types.jl b/src/types.jl index e1353b6..1967a45 100644 --- a/src/types.jl +++ b/src/types.jl @@ -29,19 +29,6 @@ abstract type AbstractRIRConfig end lrng::R = GLOBAL_RNG # Local randon number generator end -function ISMConfig( - order=(0, -1), - fs=16000.0, - N=4000, - Wd=8e-3, - hp=true, - isd=0.0, - lrng=GLOBAL_RNG -) - ISMConfig(order, fs, N, Wd, hp, isd, lrng) -end - - struct Node{T<:Real} rx::TxRxArray{T} fs::T