From f00ca2cbc973c0ccbd32c1217099968299ab1f2e Mon Sep 17 00:00:00 2001
From: zymon <s@zymon.org>
Date: Mon, 10 Jun 2024 16:29:11 +0200
Subject: [PATCH] insert_impuse

---
 src/ISM.jl | 41 +++++++++++++++++++++++------------------
 1 file changed, 23 insertions(+), 18 deletions(-)

diff --git a/src/ISM.jl b/src/ISM.jl
index f38f547..6c7cf34 100644
--- a/src/ISM.jl
+++ b/src/ISM.jl
@@ -108,14 +108,14 @@ function image_source_generator(
                 rx_DoA = isp ./ d
 
                 # Compute receiver directivity gain
-                rx_DG = directivity_pattern(SVector{3}(rx_DoA), rx.B, rx.directivity)
+                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(SVector{3}(tx_DoA), tx.B, tx.directivity)
+                tx_DG = directivity_pattern(tx_DoA, tx.B, tx.directivity)
 
                 # Compute attenuation coefficient
                 A = tx_DG * rx_DG * prod(b) / (4π * d)
@@ -127,6 +127,23 @@ function image_source_generator(
 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.
@@ -147,22 +164,10 @@ function ISM!(
     ω = 2π / config.Wd
     twid = π * config.fs
 
-    for a in image_source_generator(tx, rx, room, config)
-        a |> isnothing && continue
-
-        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
+    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