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.1"
+version = "0.3.0"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
@@ -13,6 +13,6 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 DSP = "0.7"
-LoopVectorization = "0.12.171"
+LoopVectorization = "0.12"
 StaticArrays = "1"
 julia = "1.6"

README.md

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

src/ISM.jl

@@ -36,36 +36,133 @@ function ISM(
     room::RectangularRoom,
     config::ISMConfig;
 )
-    h = zeros(config.N)
-    # TODO: h = zeros(T, config.N)
-    ISM!(h, tx, rx, room, config)
+    h = ISM_RectangularRoom_core(
+        tx.position,     # transmitter position
+        tx.directivity,  # transmitter directivity pattern
+        tx.B,            # transmitter orientation
+        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)
+        # TODO: Zmień to na funkcie operującą na zaalokowanym już h
+        # AllenBerkley_highpass100!(h, config.fs)
+        # return h
+    else
+        return h
+    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
+
+"""
+
+"""
+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!(
+    h::AbstractVector{<:T},
+    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
+    Wd::T,                              # Window width
+    ISD::T,                             # Random displacement of image source
+    lrng::AbstractRNG;                  # random number generator
+) where {T<:AbstractFloat, I<:Integer}
+
+    # Number of samples in impulose response
+    Nh = length(h)
 
     # Samples to distance coefficient [m]
-    Γ = room.c / config.fs
+    Γ = c / fs
 
     # Transform size of the room from meters to samples
-    Lₛ = room.L ./ Γ
+    Lₛ = L ./ Γ
 
     # Compute maximal wall reflection
     N = ceil.(Int, Nh ./ (2 .* Lₛ))
 
-    o_min, o_max = config.order
+    o_min, o_max = order
 
-    return (begin
+    # 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
+        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
@@ -73,91 +170,59 @@ function image_source_generator(
 
             if o_min <= o && (o <= o_max || o_max == -1)
                 # Compute Rp part
-            Rp = @. (1 - 2p) * tx.position - rx.position
+                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
+                isp .= Rp .+ Rr
 
-            if config.isd > 0.0 && o > 0
+                if ISD > 0.0 && o > 0
                     # generate random displacement for the image source
-                isp = isp .+ randn(config.lrng) .* config.isd
+                    randn!(lrng, rd)
+                    isp .+= rd .* ISD
                 end
 
                 # Distance between receiver and image source
                 d = norm(isp)
 
                 # Propagation time between receiver and image source
-            τ = d / room.c
+                τ = d / c
 
-            if τ <= Nh / config.fs # Check if it still in transfer function rangΩ
+                if τ <= Nh / 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 .= β .^ 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(rx_DoA, rx.B, rx.directivity)
+                    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
+                    tx_DoA .= -rx_DoA .* (-1, -1, -1).^perm
 
                     # Compute transmitter directivity gain
-                tx_DG = directivity_pattern(tx_DoA, tx.B, tx.directivity)
+                    tx_DG = directivity_pattern(SVector{3}(tx_DoA), tx_B, tx_dp)
 
                     # 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
+                    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) * Δ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
+                        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
 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!(
-    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 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
-    config.hp && AllenBerkley_highpass100!(h, config.fs)
-    return h
-end
-

src/RoomAcoustics.jl

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

src/TxRxModels.jl

@@ -1,6 +1,5 @@
 module TxRxModels
 
-using Statistics
 using LinearAlgebra
 using StaticArrays
 
@@ -199,26 +198,11 @@ 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                    # Sound wave velocity
-    L::Tuple{T,T,T}         # size of room (x, y, z)
-    β::Tuple{T,T,T,T,T,T}   # absorption coeffictions for walls
+    c::T
+    L::Tuple{T,T,T}
+    β::Tuple{T,T,T,T,T,T}
 end
 
 
@@ -19,16 +19,29 @@ abstract type AbstractRIRConfig end
 """
 
 """
-@kwdef struct ISMConfig{T<:Real,I<:Integer,R<:AbstractRNG} <: AbstractRIRConfig
-    order::Tuple{I,I} = (0, -1) # Order of reflection [low, high]
-    fs::T = 16e3                # Sampling frequency
-    N::I  = 4000                # Number of samples in impulse response
-    Wd::T = 8e-3                # Single impulse width
-    hp::Bool = true             # High pass filter
-    isd::T = 0.0                # Image source distortion (randomized image method)
-    lrng::R = GLOBAL_RNG        # Local randon number generator
+struct ISMConfig{T<:Real,I<:Integer,R<:AbstractRNG} <: AbstractRIRConfig
+    order::Tuple{I,I}  # Order of reflection [low, high]
+    fs::T               # Sampling frequency
+    N::I                # Number of samples in impulse response
+    Wd::T               # Single impulse width
+    hp::Bool            # High pass filter
+    isd::T              # Image source distortion (randomized image method)
+    lrng::R
 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,26 +35,29 @@ 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π / fs * 100
+    W = 2π*100/fs
 
     R1 = exp(-W)
     B1 = 2*R1*cos(W)
     B2 = -R1 * R1   
     A1 = -(1+R1)
 
-    for i = eachindex(x)
+    for i = 1:length(x)
         Y[3] = Y[2]
         Y[2] = Y[1]
         Y[1] = B2*Y[3] + B1*Y[2] + x[i]
-        x[i] = Y[1] + A1*Y[2] + R1*Y[3]
+        o[i] = Y[1] + A1*Y[2] + R1*Y[3]
     end
 
-    return x
+    return o
 end