From 727bed8060a5b37ad3802ac6f8e1c583fc909f6e Mon Sep 17 00:00:00 2001 From: Oliver Lylloff Date: Wed, 27 Nov 2024 15:46:16 +0100 Subject: [PATCH] added threads to csm and removed LazyArrays dep --- Project.toml | 5 ++--- src/AeroAcoustics.jl | 1 - src/csm.jl | 38 ++++++++++++++++++++------------------ 3 files changed, 22 insertions(+), 22 deletions(-) diff --git a/Project.toml b/Project.toml index 54248df..d5e2d64 100644 --- a/Project.toml +++ b/Project.toml @@ -4,13 +4,12 @@ keywords = ["AeroAcoustics"] license = "MIT" desc = "A package for post-processing of microphone array measurements" authors = ["Oliver Lylloff "] -version = "0.2.4" +version = "0.2.5" [deps] DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" -LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" @@ -21,9 +20,9 @@ ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d" DSP = "0.6.7, 0.7" Distances = "0.10" FFTW = "1.1" -LazyArrays = "0.16, 0.17, 0.18, 0.19, 0.20, 0.21, 0.22, 1" NLsolve = "4.4" Parameters = "0.12.3" +Statistics = "1" ThreadsX = "0.1.10" julia = "1.6" diff --git a/src/AeroAcoustics.jl b/src/AeroAcoustics.jl index 3a931b3..a9b1afe 100644 --- a/src/AeroAcoustics.jl +++ b/src/AeroAcoustics.jl @@ -4,7 +4,6 @@ using LinearAlgebra using Statistics using NLsolve using Parameters -using LazyArrays using FFTW import DSP using ThreadsX diff --git a/src/csm.jl b/src/csm.jl index 4f13f72..3ba8948 100644 --- a/src/csm.jl +++ b/src/csm.jl @@ -25,14 +25,10 @@ end # flat_t is a helper function to reshape a S x M matrix to a vector of vectors function flat_t(t::AbstractArray{T,N}) where {T <: AbstractFloat, N} - if size(t,1) < size(t,2) + if size(t, 1) < size(t, 2) t = permutedims(t) end - ta = Array{Vector{T}}(undef,size(t,2)) - for i in axes(t,2) - ta[i] = view(t,:,i) - end - return ta + return map(i -> view(t, :, i), axes(t, 2)) end """ @@ -45,7 +41,8 @@ function csm(t::AbstractArray{T};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n csm(flat_t(t);n=n,noverlap=noverlap,fs=fs,win=win,scaling=scaling) end -function csm(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum") where T <: AbstractFloat +function csm(t::Vector{<:AbstractVector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum",multi_thread=true) where T <: AbstractFloat + _foreach = AeroAcoustics.check_multithread(multi_thread) M = length(t) Nf = div(n,2)+1 nout = div((length(t[1]) - n), n - noverlap)+1 @@ -53,9 +50,12 @@ function csm(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning( C = Array{Complex{T}}(undef,M,M,Nf) S = DSP.stft.(t, n, noverlap; fs=fs, window=win, onesided=true) Sc = conj.(S) - for m in 1:M - for j in m:M - C[j,m,:] .= dropdims(mean(LazyArray(@~ Sc[m].*S[j]),dims=2);dims=2) + @views @inbounds _foreach(1:Nf) do ω + for m in 1:M + for j in m:M + C[j, m, ω] = mean(Sc[m][ω,:] .* S[j][ω,:]) + C[m, j, ω] = conj(C[j, m, ω]) + end end end @@ -68,9 +68,6 @@ function csm(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning( C .*= scale C[:,:,2:end-1] .*= 2 - for ω in 1:Nf - C[:,:,ω] = Hermitian(C[:,:,ω],:L) - end return FreqArray(C,fc) end @@ -85,7 +82,7 @@ function csm_slow(t::AbstractArray{T};n=1024,noverlap=div(n,2),fs=1,win=DSP.hann csm_slow(flat_t(t);n=n,noverlap=noverlap,fs=fs,win=win,scaling=scaling) end -function csm_slow(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum") where T <: AbstractFloat +function csm_slow(t::Vector{<:AbstractVector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum") where T <: AbstractFloat M = length(t) Nf = div(n,2)+1 nout = div((length(t[1]) - n), n - noverlap)+1 @@ -94,14 +91,19 @@ function csm_slow(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.han C = Array{Complex{T}}(undef,M,M,Nf) S = Array{Complex{T}}(undef,Nf,nout) Sc = similar(S) + for m in 1:M - stft!(Sc,t[m],n,noverlap; fs=fs, window=win, onesided=true) - conj!(Sc) + stft!(Sc, t[m], n, noverlap; fs=fs, window=win, onesided=true) + conj!(Sc) + for j in m:M - stft!(S,t[j],n,noverlap; fs=fs, window=win, onesided=true) - C[j,m,:] .= dropdims(mean(LazyArray(@~ Sc.*S),dims=2);dims=2) + stft!(S, t[j], n, noverlap; fs=fs, window=win, onesided=true) + product = Sc .* S + mean_product = mean(product, dims=2) + C[j, m, :] .= dropdims(mean_product; dims=2) # Drop the singleton dimension end end + if scaling == "density" scale = 1/(fs*sum(abs2,win))