diff --git a/src/ISM.jl b/src/ISM.jl
index d8601eb..d851d87 100644
--- a/src/ISM.jl
+++ b/src/ISM.jl
@@ -86,14 +86,6 @@ function ISM!(
 
     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
@@ -106,17 +98,14 @@ function ISM!(
 
             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
+                Rp = @. (1 - 2p) * tx.position - rx.position
 
                 # image source position for given permutation
-                isp .= Rp .+ Rr
+                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
+                    isp .+= randn(config.lrng) * config.isd
                 end
 
                 # Distance between receiver and image source
@@ -128,17 +117,17 @@ function ISM!(
                 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))
+                    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
+                    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
+                    tx_DoA = @. -rx_DoA * (-1, -1, -1)^perm
 
                     # Compute transmitter directivity gain
                     tx_DG = directivity_pattern(SVector{3}(tx_DoA), tx.B, tx.directivity)
@@ -163,3 +152,112 @@ function ISM!(
         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
+
+    # Main loop
+    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(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)
+
+                (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 ISM_slower!(
+    h::AbstractVector{T},
+    tx::TxRx,
+    rx::TxRx,
+    room::RectangularRoom,
+    config::ISMConfig;
+) where {T<:AbstractFloat}
+
+    # Impulse parameters
+    Δt = 1 / config.fs
+    ω = 2π / config.Wd
+    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
+
+            # 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
+