commit 5c81c7802592a41cecdf3f75bdbcbaa8a4317c92 Author: zymon Date: Thu Oct 13 15:50:42 2022 +0200 init diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..ba39cc5 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +Manifest.toml diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..b464027 --- /dev/null +++ b/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. diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..3b41c7b --- /dev/null +++ b/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" diff --git a/README.md b/README.md new file mode 100644 index 0000000..ddc25b6 --- /dev/null +++ b/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); + + +``` diff --git a/src/ISM.jl b/src/ISM.jl new file mode 100644 index 0000000..0bf25f3 --- /dev/null +++ b/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 diff --git a/src/RoomAcoustics.jl b/src/RoomAcoustics.jl new file mode 100644 index 0000000..37ea793 --- /dev/null +++ b/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 diff --git a/src/directivity.jl b/src/directivity.jl new file mode 100644 index 0000000..fc02e65 --- /dev/null +++ b/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 + diff --git a/src/types.jl b/src/types.jl new file mode 100644 index 0000000..ff7a829 --- /dev/null +++ b/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 + + diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..34fd73f --- /dev/null +++ b/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