Refactor of the ISM method for rectangular room.
Squashed commit of the following: commitf00ca2cbc9Author: zymon <s@zymon.org> Date: Mon Jun 10 16:29:11 2024 +0200 insert_impuse commit9c922f79bfAuthor: zymon <s@zymon.org> Date: Tue Jun 4 16:26:04 2024 +0200 aktu commit88ba6ea0b0Author: zymon <s@zymon.org> Date: Sun Jun 2 23:25:11 2024 +0200 generator of image sources
This commit is contained in:
parent
7505a5f5fc
commit
a9200e1106
1 changed files with 98 additions and 90 deletions
188
src/ISM.jl
188
src/ISM.jl
|
|
@ -51,8 +51,100 @@ function ISM(
|
||||||
end
|
end
|
||||||
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:
|
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.
|
[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.
|
[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;
|
config::ISMConfig;
|
||||||
) where {T<:AbstractFloat}
|
) 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
|
# Impulse parameters
|
||||||
Δt = 1 / config.fs
|
Δt = 1 / config.fs
|
||||||
ω = 2π / config.Wd
|
ω = 2π / config.Wd
|
||||||
twid = π * config.fs
|
twid = π * config.fs
|
||||||
|
|
||||||
# Compute maximal wall reflection
|
for image_source ∈ image_source_generator(tx, rx, room, config)
|
||||||
N = ceil.(Int, Nh ./ (2 .* Lₛ))
|
image_source |> isnothing && continue
|
||||||
|
A, τ, _, _ = image_source
|
||||||
o_min, o_max = config.order
|
insert_impuse!(h, A, τ, config.Wd, config.fs, Δt, twid, ω, config.N)
|
||||||
|
|
||||||
# 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
|
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Reference in a new issue