Commit 069dc7ce authored by ollyl's avatar ollyl

Merge branch 'new_env' into 'master'

Various updates and improvements

See merge request !9
parents 13d4efb7 48f123cd
......@@ -17,7 +17,6 @@ export Environment,
steeringvectors,
steeringvectors!,
csm,
csm!,
beamforming,
psf,
damas,
......
......@@ -5,10 +5,10 @@ Calculate frequency-domain beamforming using cross-spectral matrix `csm` and st
`csm` must be a square (Hermitian) matrix optionally with a third (frequency) dimension.
First dimension of `v` and `csm` must be equal.
"""
function beamforming(csm::AbstractArray{Complex{T},N},v::AbstractArray{Complex{T},N}) where {T <: AbstractFloat,N}
function beamforming(csm::A,v::B) where {T1<:AbstractFloat,T2<:AbstractFloat, N, A<:AbstractArray{Complex{T1},N},B<:AbstractArray{Complex{T2},N}}
@assert size(csm,1) == size(csm,2) "csm must be square with dimensions M x M (x Nf)!"
@assert size(v,1) == size(csm,1) "First dimension of v and csm must be equal!"
b = Array{T, 2}(undef, size(v,2), size(csm,3))
b = Array{T2, 2}(undef, size(v,2), size(csm,3))
for j in axes(csm,3)
csmd = @view csm[:,:,j]
for i in axes(v,2)
......@@ -27,16 +27,15 @@ stored in Environment struct
"""
function beamforming(E::Environment)
@unpack CSM,steeringvec,M,N,Nf,fn,Cinds = E
T = eltype(CSM)
b = Array{Float64}(undef, N, Nf)
csmd = Array{T}(undef,M,M)
vd = Array{T}(undef,M)
for j in 1:Nf
csmd .= selectdim(CSM.arr[:,:,Cinds], 3, j) # a view
for i in 1:N
vd .= @view steeringvec.arr[:,i,j]
b[i,j] = real(vd'*csmd*vd)
end
@views @inbounds for j in 1:Nf
AeroAcoustics.bf_col!(b[:,j],steeringvec.arr[:,:,j],selectdim(CSM.arr[:,:,Cinds], 3, j))
end
return FreqArray(b,fn)
end
function bf_col!(b::A,st::B,csm::C) where {T1<:AbstractFloat,T2<:AbstractFloat, N, A<:AbstractArray{T1,1}, B<:AbstractArray{Complex{T1},N},C<:AbstractArray{Complex{T2},N}}
@views @inbounds for i in 1:length(b)
b[i] = real(st[:,i]'*csm*st[:,i])
end
end
......@@ -75,17 +75,17 @@ function csm(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(
end
"""
csm_slow(t;n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n))
csm_slow(t;n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum")
Calculate cross-spectral matrix from time series `t` which is `S x M` dimensional,
where `S` is the number of samples and `M`the number of microphones. This version is slower than
`csm` but allocates much less (in some cases), therefore `csm_slow` could be used for long time signals.
"""
function csm_slow(t::AbstractArray{T};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n)) where T <: AbstractFloat
csm_slow(flat_t(t);n=n,noverlap=noverlap,fs=fs,win=win)
function csm_slow(t::AbstractArray{T};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n),scaling="spectrum") where T <: AbstractFloat
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)) where T <: AbstractFloat
function csm_slow(t::Vector{Vector{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,19 +94,23 @@ 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)
ds = similar(S)
Pxy = Array{Complex{T}}(undef,Nf)
for m in 1:M
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)
broadcast!(*,S,Sc,S)
mean!(Pxy,S)
C[j,m,:] = Pxy
C[j,m,:] .= dropdims(mean(LazyArray(@~ Sc.*S),dims=2);dims=2)
end
end
C .*= 2/n/weight
if scaling == "density"
scale = 1/(fs*sum(abs2,win))
elseif scaling == "spectrum"
scale = 1/(sum(win).^2)
end
C .*= scale
C[:,:,2:end-1] .*= 2
for ω in 1:Nf
C[:,:,ω] = Hermitian(C[:,:,ω],:L)
......
......@@ -7,9 +7,8 @@ steeringvectors. Optionally, supply the index where the psf is centered, default
function psf(E::Environment,cent::Int64=floor(Int,E.N/2)+1)
@unpack steeringvec,M,N,Nf,fn = E
p = Array{Float64, 2}(undef, N, Nf)
for j in 1:Nf
vc = view(steeringvec.arr,:,:,j)
p[:,j] .= AeroAcoustics.psf_col!(p[:,j],vc,cent)
@views @inbounds for j in 1:Nf
AeroAcoustics.psf_col!(p[:,j],steeringvec.arr[:,:,j],cent)
end
return FreqArray(p,fn)
end
......@@ -19,7 +18,7 @@ end
Calculate single frequency point spread function as a column vector.
"""
function psf_col!(p,steeringvec,cent)
function psf_col!(p::A,steeringvec::B,cent::Int64) where {T<:AbstractFloat, N, A<:AbstractArray{T,1}, B<:AbstractArray{Complex{T},N}}
M = size(steeringvec,1)
@views p .= 0.5.*M^2 .*abs2.(steeringvec'*steeringvec[:,cent])
end
......@@ -19,13 +19,15 @@ import DSP
n = 1024
Pxx_d = csm(x;n=n,noverlap=div(n,2),win=DSP.hanning(n),fs=fs,scaling="density")
fc = Pxx_d.fc
Pxx_d = real.(Pxx_d[1,1,:])
Pxx_d1 = real.(Pxx_d[1,1,:])
Pxx_w = DSP.welch_pgram(x,n;onesided=true,fs=fs,window=DSP.hanning(n))
@test Pxx_d Pxx_w.power
@test Pxx_d1 Pxx_w.power
@test Pxx_d == AeroAcoustics.csm_slow(x;n=n,noverlap=div(n,2),win=DSP.hanning(n),fs=fs,scaling="density")
Pxx_s = csm(x;n=n,noverlap=div(n,2),win=DSP.hanning(n),fs=fs,scaling="spectrum")
fc = Pxx_s.fc
Pxx_s = real.(Pxx_s[1,1,:])
Pxx_s1 = real.(Pxx_s[1,1,:])
ENBW = enbw(fs,DSP.hanning(1024))
@test Pxx_s Pxx_w.power*ENBW
@test Pxx_s1 Pxx_w.power*ENBW
@test Pxx_s == AeroAcoustics.csm_slow(x;n=n,noverlap=div(n,2),win=DSP.hanning(n),fs=fs,scaling="spectrum")
end
......@@ -34,6 +34,8 @@ import DSP
steeringvectors!(env)
b = beamforming(env)
bd = beamforming(env.CSM.arr[:,:,env.Cinds],env.steeringvec.arr)
@test b.arr == bd
idx = 1 # Frequency index
s1,p1 = findmax(reshape(b[:,idx],env.Nx,env.Ny))
bmax = ceil(SPL(sqrt(2).*s1))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment