Friends with LoopVectorization once again

Squashed commit of the following:

commit f00ca2cbc9
Author: zymon <s@zymon.org>
Date:   Mon Jun 10 16:29:11 2024 +0200

    insert_impuse

commit 9c922f79bf
Author: zymon <s@zymon.org>
Date:   Tue Jun 4 16:26:04 2024 +0200

    aktu

commit 88ba6ea0b0
Author: zymon <s@zymon.org>
Date:   Sun Jun 2 23:25:11 2024 +0200

    generator of image sources

Project.toml

@@ -1,7 +1,7 @@
 name = "RoomAcoustics"
 uuid = "9b22aa7e-b0d0-4fe8-9c3b-2b8bf774f735"
 authors = ["Szymon M. Woźniak <s@zymon.org>"]
-version = "0.3.0"
+version = "0.3.1"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
@@ -13,6 +13,6 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 DSP = "0.7"
-LoopVectorization = "0.12"
+LoopVectorization = "0.12.171"
 StaticArrays = "1"
 julia = "1.6"

README.md

@@ -4,7 +4,7 @@
 
 
 ```julia
-] add https://codeberg.org/zymon/RoomAcoustics.jl
+] add RoomAcoustics
 ```
 
 

src/ISM.jl

@@ -36,193 +36,128 @@ function ISM(
     room::RectangularRoom,
     config::ISMConfig;
 )
-    h = ISM_RectangularRoom_core(
+    h = zeros(config.N)
-        tx.position,     # transmitter position
+    # TODO: h = zeros(T, config.N)
-        tx.directivity,  # transmitter directivity pattern
+    ISM!(h, tx, rx, room, config)
-        tx.B,            # transmitter orientation
+end
-        rx.position,     # receiver position
-        rx.directivity,  # receiver directivity pattern
-        rx.B,            # receiver orientation
-        room.L,          # room size
-        room.β,          # reflection coefficients
-        room.c,          # velocity of the sound
-        config.fs,       # sampling frequency
-        config.order,    # order of reflections
-        config.N,        # length of response
-        config.Wd,       # single impulse width
-        config.isd,      # random image source displacement
-        config.lrng,     # Local random number generator
-    )
 
-    if config.hp
+
-        return AllenBerkley_highpass100(h, config.fs)
+function image_source_generator(
-        # TODO: Zmień to na funkcie operującą na zaalokowanym już h
+    tx::TxRx,
-        # AllenBerkley_highpass100!(h, config.fs)
+    rx::TxRx,
-        # return h
+    room::RectangularRoom,
-    else
+    config::ISMConfig;
-        return h
+)
+    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 = 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
+    @turbo for i ∈ i_s:i_e
+        t = (i - 1) * Δt - τ # time signature
+        w = 0.5 * (1.0 + cos(ω*t)) # Hann window
+        x = twid * t
+        sinc = ifelse(iszero(x), 1.0,  sin(x)/x)
+        h[i] +=  w * A * sinc
     end
 end
 
 
-
 """
-
-"""
-function ISM_RectangularRoom_core(
-    tx_p::SVector{3,T},                 # transmitter position
-    tx_dp::AbstractDirectivityPattern,  # transmitter directivity pattern
-    tx_B::SMatrix{3,3,T},               # receiver orientation
-    rx_p::SVector{3,T},                 # reveiver position
-    rx_dp::AbstractDirectivityPattern,  # receiver directivity pattern
-    rx_B::SMatrix{3,3,T},               # receiver orientation
-    L::Tuple{T,T,T},                    # room size (Lx, Ly, Lz)
-    β::Tuple{T,T,T,T,T,T},              # reflection coefficients (βx1, βx2, βy1, βy2, βz1, βz2)
-    c::T,                               # velocity of the wave
-    fs::T,                              # sampling frequeyncy
-    order::Tuple{I,I},                  # order of reflections; (min, max)
-    Nh::Integer,                        # number of returned samples
-    Wd::T,                              # window width
-    ISD::T,                             # random displacement of image source
-    lrng::AbstractRNG                   # random number generator
-)::AbstractVector{T} where {T<:AbstractFloat, I<:Integer}
-    # Allocate memory for the impulose response
-    h = zeros(T, Nh)
-
-    # Call
-    ISM_RectangularRoom_core!(
-        h,      # container for impulse response
-        tx_p,   # transmitter position
-        tx_dp,  # transmitter directivity pattern
-        tx_B,   # transmitter orientation
-        rx_p,   # receiver position
-        rx_dp,  # receiver directivity pattern
-        rx_B,   # receiver orientation
-        L,      # room size
-        β,      # reflection coefficients
-        c,      # velocity of the sound
-        fs,     # sampling frequency
-        order,  # order of reflections
-        Wd,     # single impulse width
-        ISD,    # random image source displacement
-        lrng,   # Local random number generator
-    )
-    h
-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_RectangularRoom_core!(
+function ISM!(
-    h::AbstractVector{<:T},
+    h::AbstractVector{T},
-    tx_p::SVector{3,T},                 # transmitter position
+    tx::TxRx,
-    tx_dp::AbstractDirectivityPattern,  # transmitter directivity pattern
+    rx::TxRx,
-    tx_B::SMatrix{3,3,T},               # receiver orientation
+    room::RectangularRoom,
-    rx_p::SVector{3,T},                 # reveiver position
+    config::ISMConfig;
-    rx_dp::AbstractDirectivityPattern,  # receiver directivity pattern
+) where {T<:AbstractFloat}
-    rx_B::SMatrix{3,3,T},               # receiver orientation
+    # Impulse parameters
-    L::Tuple{T,T,T},                    # room size (Lx, Ly, Lz)
+    Δt = 1 / config.fs
-    β::Tuple{T,T,T,T,T,T},              # Reflection coefficients (βx1, βx2, βy1, βy2, βz1, βz2)
+    ω = 2π / config.Wd
-    c::T,                               # velocity of the wave
+    twid = π * config.fs
-    fs::T,                              # sampling frequeyncy
+    for image_source ∈ image_source_generator(tx, rx, room, config)
-    order::Tuple{I,I},                  # order of reflections; min max
+        image_source |> isnothing && continue
-    Wd::T,                              # Window width
+        A, τ, _, _ = image_source
-    ISD::T,                             # Random displacement of image source
+        insert_impuse!(h, A, τ, config.Wd, config.fs, Δt, twid, ω, config.N)
-    lrng::AbstractRNG;                  # random number generator
-) where {T<:AbstractFloat, I<:Integer}
-
-    # Number of samples in impulose response
-    Nh = length(h)
-
-    # Samples to distance coefficient [m]
-    Γ = c / fs
-
-    # Transform size of the room from meters to samples
-    Lₛ = L ./ Γ
-
-    # Compute maximal wall reflection
-    N = ceil.(Int, Nh ./ (2 .* Lₛ))
-
-    o_min, o_max = 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
-        Rr = 2 .* r .* 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_p[i] - rx_p[i]
-                end
-
-                # image source position for given permutation
-                isp .= Rp .+ Rr
-
-                if ISD > 0.0 && o > 0
-                    # generate random displacement for the image source
-                    randn!(lrng, rd)
-                    isp .+= rd .* ISD
-                end
-
-                # Distance between receiver and image source
-                d = norm(isp)
-
-                # Propagation time between receiver and image source
-                τ = d / c
-
-                if τ <= Nh / fs # Check if it still in transfer function rangΩ
-
-                    # Compute value of reflection coefficients
-                    b .= β .^ 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_dp)
-
-                    # 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_dp)
-
-                    # Compute attenuation coefficient
-                    A = tx_DG * rx_DG * prod(b) / (4π * d)
-
-                    # Compute range of samples in transfer function
-                    i_s = max(ceil(Int, (τ - Wd / 2) * fs) + 1, 1)  # start
-                    i_e = min(floor(Int, (τ + Wd / 2) * fs) + 1, Nh) # end
-
-                    # Insert yet another impulse into transfer function
-                    @turbo for i ∈ i_s:i_e
-                        t = (i - 1) / fs - τ # time signature
-                        w = 0.5 * (1.0 + cos(2π * t / Wd)) # Hann window
-                        x = π * fs * t + eps()
-                        h[i] += w * A * sin(x)/x
-                    end
-                end
-            end
-        end
     end
+    config.hp && AllenBerkley_highpass100!(h, config.fs)
+    return h
 end
+

src/RoomAcoustics.jl

@@ -12,7 +12,9 @@ using DSP: conv
 include("TxRxModels.jl")
 export TxRxModels
 using .TxRxModels
-using .TxRxModels: AbstractDirectivityPattern
+using .TxRxModels:
+    AbstractDirectivityPattern,
+    AbstractTxRx
 
 include("types.jl")
 include("utils.jl")

src/TxRxModels.jl

@@ -1,5 +1,6 @@
 module TxRxModels
 
+using Statistics
 using LinearAlgebra
 using StaticArrays
 
@@ -198,11 +199,26 @@ function circular_array(N::Integer, r::Real)
     [SVector{3}([r * cos(α), r * sin(α), 0]) for α ∈ 0:Δα:2π-Δα]
 end
 
+function square_array(N::I, M::I, Lx::T, Ly::T) where {I<:Integer, T<:Real}
+    Δx, Δy = Lx / N, Ly / N
+    p = [SVector{3}([Δx*i, Δy*j, 0.0]) for i = 1:N for j = 1:M]
+    p .- [mean(p)]
+end
+
 function fibonacci_array(N::Integer, r::Real)
     P = fibonacci_sphere(N)
     [SVector{3}(r .* p) for p in P]
 end
 
+function tetrahedron_array(r=1.0)
+    [
+        SVector{3}([+1,+1,+1] * r / √3 ),
+        SVector{3}([+1,-1,-1] * r / √3 ),
+        SVector{3}([-1,+1,-1] * r / √3 ),
+        SVector{3}([-1,-1,+1] * r / √3 ),
+    ]
+end
+
 physical_array = (
     matrix_voice = (
         cartesian = [

src/types.jl

@@ -7,9 +7,9 @@ export Node, Event
 abstract type AbstractRoom end
 
 struct RectangularRoom{T<:Real} <: AbstractRoom
-    c::T
+    c::T                    # Sound wave velocity
-    L::Tuple{T,T,T}
+    L::Tuple{T,T,T}         # size of room (x, y, z)
-    β::Tuple{T,T,T,T,T,T}
+    β::Tuple{T,T,T,T,T,T}   # absorption coeffictions for walls
 end
 
 
@@ -19,29 +19,16 @@ abstract type AbstractRIRConfig end
 """
 
 """
-struct ISMConfig{T<:Real,I<:Integer,R<:AbstractRNG} <: AbstractRIRConfig
+@kwdef struct ISMConfig{T<:Real,I<:Integer,R<:AbstractRNG} <: AbstractRIRConfig
-    order::Tuple{I,I}  # Order of reflection [low, high]
+    order::Tuple{I,I} = (0, -1) # Order of reflection [low, high]
-    fs::T               # Sampling frequency
+    fs::T = 16e3                # Sampling frequency
-    N::I                # Number of samples in impulse response
+    N::I  = 4000                # Number of samples in impulse response
-    Wd::T               # Single impulse width
+    Wd::T = 8e-3                # Single impulse width
-    hp::Bool            # High pass filter
+    hp::Bool = true             # High pass filter
-    isd::T              # Image source distortion (randomized image method)
+    isd::T = 0.0                # Image source distortion (randomized image method)
-    lrng::R
+    lrng::R = GLOBAL_RNG        # Local randon number generator
 end
 
-function ISMConfig(
-    order=(0, -1),
-    fs=16000.0,
-    N=4000,
-    Wd=8e-3,
-    hp=true,
-    isd=0.0,
-    lrng=GLOBAL_RNG
-)
-    ISMConfig(order, fs, N, Wd, hp, isd, lrng)
-end
-
-
 struct Node{T<:Real}
     rx::TxRxArray{T}
     fs::T

src/utils.jl

@@ -35,29 +35,26 @@ function Sabine_RT60(T60, L::Tuple, c)
     sqrt(1-α)
 end
 
-
-
 """
 b = [1, -B1, -B2]
 a = [1, A1, R1] 
 """
-function AllenBerkley_highpass100(x, fs)
+function AllenBerkley_highpass100!(x, fs)
-    o = x .* 0
     Y = zeros(3)
 
-    W = 2π*100/fs
+    W = 2π / fs * 100
 
     R1 = exp(-W)
     B1 = 2*R1*cos(W)
-    B2 = -R1 * R1   
+    B2 = -R1 * R1
     A1 = -(1+R1)
 
-    for i = 1:length(x)
+    for i = eachindex(x)
         Y[3] = Y[2]
         Y[2] = Y[1]
         Y[1] = B2*Y[3] + B1*Y[2] + x[i]
-        o[i] = Y[1] + A1*Y[2] + R1*Y[3]
+        x[i] = Y[1] + A1*Y[2] + R1*Y[3]
     end
 
-    return o
+    return x
 end