init

.gitignore

@@ -0,0 +1 @@
+Manifest.toml

LICENSE

@@ -0,0 +1,19 @@
+Copyright (c) 2022 Szymon M. Woźniak
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

Project.toml

@@ -0,0 +1,10 @@
+name = "RoomAcoustics"
+uuid = "9b22aa7e-b0d0-4fe8-9c3b-2b8bf774f735"
+authors = ["zymon"]
+version = "0.1.0"
+
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

README.md

@@ -0,0 +1,46 @@
+# RoomAcoustics.jl
+
+`RoomAcoustics.jl` is a Julia package for acoustics simulations of the rooms.
+
+
+```julia
+] add https://git.sr.ht/~zymon/RoomAcoustics.jl
+```
+
+
+Currently, supported methods:
+* Image Source for rectangular (shoebox) rooms
+
+
+# Example
+
+```julia
+using StaticArrays
+using LinearAlgebra
+
+using RoomAcoustics
+
+
+c = 343.0;
+fs = 16000.0;
+rir_Nsamples = 4000;
+β = 0.75;
+room_β = (β, β, β, β, β, β);
+room_L = (10., 10., 3.);
+
+
+mic = SVector{3}([5., 5., 1.]);
+source = SVector{3}([1., 9., 2.]);
+
+
+# Setup configuration
+room = RectangularRoom(c, room_L, room_β);
+rir_config = ISMConfig((0, -1), fs, rir_Nsamples, 8e-3, true, 0.0);
+rx = TxRx(mic);
+tx = TxRx(source);
+
+# Compute transfer function using Image Source Method
+h = ISM(rx, tx, room, rir_config);
+
+
+```

src/ISM.jl

@@ -0,0 +1,198 @@
+export ISM
+
+
+"""
+
+"""
+function ISM(
+    rxs::AbstractVector{<:TxRx},
+    tx::TxRx,
+    room::Room,
+    config::ISMConfig;
+)
+	[ISM(rx, tx, room, config) for rx in rxs]
+end
+
+
+
+
+"""
+
+"""
+function ISM(
+    rx::TxRx,
+    tx::TxRx,
+    room::RectangularRoom,
+    config::ISMConfig;
+)
+    h = ISM_RectangularRoom_core(
+        tx.position,     # transmitter position
+        rx.position,     # receiver position
+        rx.B,            # receiver orientation
+        rx.directivity,  # receiver directivity pattern
+        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)
+        # Zmień to na funkcie operującą na zaalokowanym już h
+        # AllenBerkley_highpass100!(h, config.fs)
+        #return h
+    else
+        return h
+    end
+end
+
+
+
+"""
+
+"""
+function ISM_RectangularRoom_core(
+    tx::SVector{3, T},                  # transmitter position
+    rx::SVector{3, T},                  # reveiver position
+    B::SMatrix{3, 3, T},                # receiver orientation
+    dp::DirectivityPattern,             # Receiver directivity pattern
+    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{<:Int, <:Int}, # order of reflections; min max
+    Nh::Integer,                        # h lenght in samples
+    Wd::T,                              # Window width
+    ISD::T,                             # Random displacement of image source
+    lrng::AbstractRNG                   # random number generator
+)::AbstractVector{T} where {T<:AbstractFloat}
+
+    # Allocate memory for impulose response
+    h = zeros(T, Nh)
+
+    # Call
+    ISM_RectangularRoom_core!(
+        h,      # Container for impulse response
+        tx,     # transmitter position
+        rx,     # receiver position
+        B,      # receiver orientation
+        dp,     # receiver directivity pattern
+        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
+    )
+
+    return h
+end
+
+"""
+
+"""
+function ISM_RectangularRoom_core!(
+    h::AbstractVector{<:T},
+    tx::SVector{3, T},                  # transmitter position
+    rx::SVector{3, T},                  # reveiver position
+    B::SMatrix{3, 3, T},                # receiver orientation
+    dp::DirectivityPattern,             # Receiver directivity pattern
+    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{<:Int, <:Int},         # order of reflections; min max
+    Wd::T,                              # Window width
+    ISD::T,                             # Random displacement of image source
+    lrng::AbstractRNG                   # random number generator
+)::AbstractVector{T} where {T<:AbstractFloat}
+
+    # 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
+    tx_isp = @MVector zeros(T, 3)   # Container for relative image source position
+    b = @MVector zeros(T, 6)        # Container for effective reflection coefficient
+    DoA = @MVector zeros(T, 3)      # Container for direction of incoming ray to receiver
+    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[i] - rx[i]
+                end
+
+                # Position of [randomized] image source for given permutation
+                tx_isp .= Rp .+ Rr
+
+                if ISD > 0.0 && o > 0
+                    # Generate random displacement for the image source
+                    randn!(lrng, rd)
+                    tx_isp .+= rd .* ISD
+                end
+
+                # Distance between receiver and image source
+                dist = norm(tx_isp)
+
+                # Propagation time between receiver and image source
+                τ = dist/c
+
+                if τ <= Nh/fs # Check if it still in transfer function range
+
+                    # Compute value of reflection coefficients
+                    b .=  β .^ abs.((n-q, n, l-j, l, m-k, m))
+
+                    # Direction of Arrival of ray
+                    DoA .= tx_isp./dist
+
+                    # Compute receiver directivity gain
+                    DG = directivity_pattern(SVector{3}(DoA), B, dp)
+
+                    # Compute attenuation coefficient
+                    A = DG * prod(b)/(4π * dist)
+
+                    # 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
+                    for i ∈ i_s:i_e
+                        t = (i-1)/fs - τ # time signature
+                        w = 0.5 * (1.0 + cos(2π*t/Wd)) # Hann window
+                        h[i] += w * A * sinc(fs * t) # sinc
+                    end
+                end
+            end
+        end
+    end
+    h
+end

src/RoomAcoustics.jl

@@ -0,0 +1,14 @@
+module RoomAcoustics
+
+using LinearAlgebra
+using StaticArrays
+using Statistics
+using Random
+using Random: GLOBAL_RNG
+
+include("types.jl")
+include("utils.jl")
+include("directivity.jl")
+include("ISM.jl")
+
+end # module

src/directivity.jl

@@ -0,0 +1,74 @@
+
+"""
+
+"""
+function cardioid_pattern(
+    d::SVector{3, <:Real},
+    B::SMatrix{3, 3, <:Real},
+    ρ::Real,
+)::Real
+    r = [1., 0., 0.]
+    ρ + (1-ρ) * r' * B' * d
+end
+
+
+
+"""
+
+"""
+function directivity_pattern(
+    d::SVector{3, <:Real},
+    B::SMatrix{3, 3, <:Real},
+    pattern::OmnidirectionalPattern,
+)::Real
+    1
+end
+
+
+"""
+
+"""
+function directivity_pattern(
+    d::SVector{3, <:Real},
+    B::SMatrix{3, 3, <:Real},
+    pattern::SubcardioidPattern,
+)::Real
+    cardioid_pattern(d, B, 0.75)
+end
+
+
+"""
+
+"""
+function directivity_pattern(
+    d::SVector{3, <:Real},
+    B::SMatrix{3, 3, <:Real},
+    pattern::CardioidPattern,
+)::Real
+    cardioid_pattern(d, B, 0.50)
+end
+
+
+"""
+
+"""
+function directivity_pattern(
+    d::SVector{3, <:Real},
+    B::SMatrix{3, 3, <:Real},
+    pattern::HypercardioidPattern,
+)::Real
+    cardioid_pattern(d, B, 0.25)
+end
+
+
+"""
+
+"""
+function directivity_pattern(
+    d::SVector{3, <:Real},
+    B::SMatrix{3, 3, <:Real},
+    pattern::BidirectionalPattern,
+)::Real
+    cardioid_pattern(d, B, 0.00)
+end
+    

src/types.jl

@@ -0,0 +1,79 @@
+export  Omnidirectional,
+        Subcardioid,
+        Cardioid,
+        Hypercardioid,
+        Bidirectional
+
+export TxRx
+
+export Room, RectangularRoom
+export RIRConfig, ISMConfig
+
+
+
+
+abstract type DirectivityPattern end
+struct OmnidirectionalPattern <: DirectivityPattern end
+struct SubcardioidPattern     <: DirectivityPattern end
+struct CardioidPattern        <: DirectivityPattern end
+struct HypercardioidPattern   <: DirectivityPattern end
+struct BidirectionalPattern   <: DirectivityPattern end
+
+
+const Omnidirectional = OmnidirectionalPattern()
+const Subcardioid     = SubcardioidPattern()
+const Cardioid        = CardioidPattern()
+const Hypercardioid   = HypercardioidPattern()
+const Bidirectional   = BidirectionalPattern()
+
+
+
+struct TxRx{T<:AbstractFloat}
+    position::SVector{3, T}            # Position
+    B::SMatrix{3, 3, T}                # Orientation
+    directivity::DirectivityPattern    # Directivity pattern
+end
+
+function TxRx(position, B=SMatrix{3,3}(1.0I), d=Omnidirectional)
+    TxRx(position, B, d)
+end
+
+
+abstract type Room end
+
+struct RectangularRoom{T<:AbstractFloat} <: Room
+    c::T
+    L::Tuple{T, T, T}
+    β::Tuple{T, T, T, T, T, T}
+end
+
+
+
+abstract type RIRConfig end
+
+"""
+
+"""
+struct ISMConfig{T<:AbstractFloat, I <: Integer} <: RIRConfig
+    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::AbstractRNG
+end
+
+function ISMConfig(
+    order=(0, -1),
+    fs=16000,
+    N=4000,
+    Wd=8e-3,
+    hp=true,
+    isd=0.0,
+    lrng=GLOBAL_RNG
+)
+    ISMConfig(order, fs, N, Wd, hp, isd, lrng)
+end
+
+

src/utils.jl

@@ -0,0 +1,62 @@
+export h2RT60, Sabine_RT60
+
+
+"""
+
+"""
+function h2RT60(h::AbstractVector{<:Number}, Fs::Real)
+    cs = cumsum(reverse(h.^2))
+    edc = 10*log10.(reverse(cs./cs[end])) # energy decay curve
+
+    ind = findfirst(edc .<= -60. )
+    if ind == nothing 
+        rt = length(h)/Fs 
+    else
+        rt = ind/Fs
+    end
+    rt, edc
+end
+
+
+"""
+"""
+function Sabine_RT60(T60, L::Tuple, c)
+    # Compute volume of the room
+    V = prod(L)
+
+    # Compute surface of the room
+    S = 2*(L[1]*L[2] + L[1]*L[3] + L[2]*L[3])
+
+    #
+    α = 24 * V * log(10)/(c * S * T60)
+
+    #
+    sqrt(1-α)
+end
+
+
+
+"""
+b = [1, -B1, -B2]
+a = [1, A1, R1] 
+"""
+function AllenBerkley_highpass100(x, fs)
+    o = x .* 0
+    Y = zeros(3)
+
+    W = 2π*100/fs
+
+    R1 = exp(-W)
+    B1 = 2*R1*cos(W)
+    B2 = -R1 * R1   
+    A1 = -(1+R1)
+
+    for i = 1:length(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]
+    end
+
+    return o
+end