diff --git a/src/ISM.jl b/src/ISM.jl index d851d87..f38f547 100644 --- a/src/ISM.jl +++ b/src/ISM.jl @@ -51,108 +51,6 @@ function ISM( 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. -[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!( - 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] - Γ = 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 - - # 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 - 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(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 - end -end - function image_source_generator( tx::TxRx, @@ -175,7 +73,6 @@ function image_source_generator( o_min, o_max = config.order - # Main loop return (begin r = (n, l, m) # Wall reflection indicator Rr = 2 .* r .* room.L # Wall lattice @@ -229,7 +126,15 @@ function image_source_generator( 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 ISM_slower!( + +""" +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!( h::AbstractVector{T}, tx::TxRx, rx::TxRx, @@ -243,20 +148,20 @@ function ISM_slower!( twid = π * config.fs for a in image_source_generator(tx, rx, room, config) - if nothing !== a - A, τ, _, _ = a - # 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, config.N) # end + a |> isnothing && continue - # 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 + A, τ, _, _ = a + # 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, config.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 end