diff --git a/src/ISM.jl b/src/ISM.jl index d8601eb..6c7cf34 100644 --- a/src/ISM.jl +++ b/src/ISM.jl @@ -51,8 +51,100 @@ function ISM( end end -""" +function image_source_generator( + tx::TxRx, + rx::TxRx, + room::RectangularRoom, + config::ISMConfig; +) + T = Float64 + # Number of samples in impulse response + Nh = config.N + + # Samples to distance coefficient [m] + Γ = room.c / config.fs + + # Transform size of the room from meters to samples + Lₛ = room.L ./ Γ + + # Compute maximal wall reflection + N = ceil.(Int, Nh ./ (2 .* Lₛ)) + + o_min, o_max = config.order + + return (begin + r = (n, l, m) # Wall reflection indicator + Rr = 2 .* r .* room.L # Wall lattice + p = @SVector [q, j, k] # Permutation tuple + + # Order of reflection generated by image source + o = sum(abs.(2 .* r .- p)) + + if o_min <= o && (o <= o_max || o_max == -1) + # Compute Rp part + Rp = @. (1 - 2p) * tx.position - rx.position + + # image source position for given permutation + isp = Rp .+ Rr + + if config.isd > 0.0 && o > 0 + # generate random displacement for the image source + isp .+= randn(config.lrng) * config.isd + end + + # Distance between receiver and image source + d = norm(isp) + + # Propagation time between receiver and image source + τ = d / room.c + + if τ <= Nh / config.fs # Check if it still in transfer function rangΩ + + # Compute value of reflection coefficients + 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(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(tx_DoA, tx.B, tx.directivity) + + # Compute attenuation coefficient + A = tx_DG * rx_DG * prod(b) / (4π * d) + + (A, τ, isp, o) + end + end + end for n = -N[1]:N[1], l = -N[2]:N[2], m = -N[3]:N[3], q ∈ 0:1, j ∈ 0:1, k ∈ 0:1) +end + + +function insert_impuse!(h, A, τ, impulse_width, fs, Δt, twid, ω, N) + # Compute range of samples in transfer function + a = impulse_width / 2 + i_s = max(ceil(Int, (τ - a) * fs) + 1, 1) # start + i_e = min(floor(Int, (τ + a) * fs) + 1, N) # end + + # Insert yet another impulse into transfer function + for i ∈ i_s:i_e + t = (i - 1) * Δt - τ # time signature + w = fma(cos_fast(ω*t), 0.5, 0.5) # Hann window + x = twid * t + sinc = ifelse(iszero(x), 1.0, sin_fast(x)/x) + @inbounds h[i] += w * A * sinc + end +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. @@ -67,99 +159,15 @@ function ISM!( config::ISMConfig; ) where {T<:AbstractFloat} - # Number of samples in impulse response - Nh = length(h) - - # Samples to distance coefficient [m] - Γ = room.c / config.fs - - # Transform size of the room from meters to samples - Lₛ = room.L ./ Γ - # Impulse parameters Δt = 1 / config.fs ω = 2π / config.Wd twid = π * config.fs - # Compute maximal wall reflection - N = ceil.(Int, Nh ./ (2 .* Lₛ)) - - o_min, o_max = config.order - - # Allocate memory - rd = @MVector zeros(T, 3) # Container for random displacement - isp = @MVector zeros(T, 3) # Container for relative image source position - b = @MVector zeros(T, 6) # Container for effective reflection coefficient - rx_DoA = @MVector zeros(T, 3) # Container for direction of incoming ray to receiver - 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] - r = (n, l, m) # Wall reflection indicator - Rr = 2 .* r .* room.L # Wall lattice - for q ∈ 0:1, j ∈ 0:1, k ∈ 0:1 - p = @SVector [q, j, k] # Permutation tuple - - # Order of reflection generated by image source - o = sum(abs.(2 .* r .- p)) - - if o_min <= o && (o <= o_max || o_max == -1) - # Compute Rp part - for i = 1:3 - Rp[i] = (1 .- 2 * p[i]) * tx.position[i] - rx.position[i] - end - - # image source position for given permutation - isp .= Rp .+ Rr - - if config.isd > 0.0 && o > 0 - # generate random displacement for the image source - 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 / room.c - - if τ <= Nh / config.fs # Check if it still in transfer function rangΩ - - # Compute value of reflection coefficients - 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.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.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, (τ - 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 - t = (i - 1) * Δt - τ # time signature - w = fma(cos_fast(ω*t), 0.5, 0.5) # Hann window - x = twid * t - sinc = ifelse(iszero(x), 1.0, sin_fast(x)/x) - @inbounds h[i] += w * A * sinc - end - end - end - end + for image_source ∈ image_source_generator(tx, rx, room, config) + image_source |> isnothing && continue + A, τ, _, _ = image_source + insert_impuse!(h, A, τ, config.Wd, config.fs, Δt, twid, ω, config.N) end end +