This commit is contained in:
zymon 2024-06-04 16:26:04 +02:00
parent 88ba6ea0b0
commit 9c922f79bf

View file

@ -51,108 +51,6 @@ function ISM(
end end
end end
"""
References:
[1] J. B. Allen and D. A. Berkley, “Image method for efficiently simulating smallroom 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( function image_source_generator(
tx::TxRx, tx::TxRx,
@ -175,7 +73,6 @@ function image_source_generator(
o_min, o_max = config.order o_min, o_max = config.order
# Main loop
return (begin return (begin
r = (n, l, m) # Wall reflection indicator r = (n, l, m) # Wall reflection indicator
Rr = 2 .* r .* room.L # Wall lattice 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 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 end
function ISM_slower!(
"""
References:
[1] J. B. Allen and D. A. Berkley, “Image method for efficiently simulating smallroom 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}, h::AbstractVector{T},
tx::TxRx, tx::TxRx,
rx::TxRx, rx::TxRx,
@ -243,7 +148,8 @@ function ISM_slower!(
twid = π * config.fs twid = π * config.fs
for a in image_source_generator(tx, rx, room, config) for a in image_source_generator(tx, rx, room, config)
if nothing !== a a |> isnothing && continue
A, τ, _, _ = a A, τ, _, _ = a
# Compute range of samples in transfer function # Compute range of samples in transfer function
i_s = max(ceil(Int, (τ - config.Wd / 2) * config.fs) + 1, 1) # start i_s = max(ceil(Int, (τ - config.Wd / 2) * config.fs) + 1, 1) # start
@ -258,6 +164,5 @@ function ISM_slower!(
@inbounds h[i] += w * A * sinc @inbounds h[i] += w * A * sinc
end end
end end
end
end end