Commit 27f91568 authored by ollyl's avatar ollyl

Julia07fixes

parent 79edc77f
......@@ -17,3 +17,7 @@
*.tmp
.DS_Store
test/data/
readpti.jl
image: julia:0.6
image: julia:1.0 # image comes from Docker hub
before_script:
- apt-get update -qq && apt-get install git hdf5-tools -y
- julia build_script.jl
- apt-get -qq update; apt-get -y install git hdf5-tools
- julia --project=@. -e "import Pkg; Pkg.build()"
default:
script:
- julia -e "Pkg.test(\"AeroAcoustics.jl\"; coverage = true)"
- julia coverage_script.jl
- julia --project=@. -e "import Pkg; Pkg.test(; coverage = true)"
- julia --project=test/coverage -e 'import Pkg; Pkg.instantiate()'
- julia --project=test/coverage test/coverage/coverage-summary.jl
# uncomment to build docs
#pages:
# stage: deploy
# script:
# - julia --project=docs -e 'using Pkg; Pkg.instantiate(); Pkg.develop(PackageSpec(path=pwd()))'
# - julia --project=docs --color=yes docs/make.jl
# - mv docs/build public # move to the directory picked up by Gitlab pages
# artifacts:
# paths:
# - public
# only:
# - master
The AeroAcoustics.jl package is licensed under the MIT "Expat" License:
> Copyright (c) 2017: ollyl.
>
> Copyright (c) 2019: ollyl.
>
> 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
......@@ -19,4 +19,4 @@ The AeroAcoustics.jl package is licensed under the MIT "Expat" License:
> 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.
>
>
[[AbstractFFTs]]
deps = ["Compat", "LinearAlgebra"]
git-tree-sha1 = "8d59c3b1463b5e0ad05a3698167f85fac90e184d"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "0.3.2"
[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
[[BinDeps]]
deps = ["Compat", "Libdl", "SHA", "URIParser"]
git-tree-sha1 = "12093ca6cdd0ee547c39b1870e0c9c3f154d9ca9"
uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee"
version = "0.8.10"
[[BinaryProvider]]
deps = ["Libdl", "Pkg", "SHA", "Test"]
git-tree-sha1 = "055eb2690182ebc31087859c3dd8598371d3ef9e"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.3"
[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "ec61a16eed883ad0cfa002d7489b3ce6d039bb9a"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "1.4.0"
[[Conda]]
deps = ["Compat", "JSON", "VersionParsing"]
git-tree-sha1 = "fb86fe40cb5b35990e368709bfdc1b46dbb99dac"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.1.1"
[[DSP]]
deps = ["AbstractFFTs", "Compat", "FFTW", "LinearAlgebra", "Polynomials", "Random", "Reexport", "SpecialFunctions"]
git-tree-sha1 = "5ec38ebc4ddf6320ad50b826eb8fd7fb521993a5"
uuid = "717857b8-e6f2-59f4-9121-6e50c889abd2"
version = "0.5.2"
[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
[[Distances]]
deps = ["LinearAlgebra", "Printf", "Random", "Statistics", "Test"]
git-tree-sha1 = "0e37d8a95bafbeb9c800ef27ab6f443d29e17610"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.7.4"
[[Distributed]]
deps = ["LinearAlgebra", "Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
[[FFTW]]
deps = ["AbstractFFTs", "BinaryProvider", "Compat", "Conda", "Libdl", "LinearAlgebra", "Reexport", "Test"]
git-tree-sha1 = "29cda58afbf62f35b1a094882ad6c745a47b2eaa"
uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
version = "0.2.4"
[[InteractiveUtils]]
deps = ["LinearAlgebra", "Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
[[JSON]]
deps = ["Dates", "Distributed", "Mmap", "Sockets", "Test", "Unicode"]
git-tree-sha1 = "1f7a25b53ec67f5e9422f1f551ee216503f4a0fa"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.20.0"
[[LibGit2]]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"
[[OrderedCollections]]
deps = ["Random", "Serialization", "Test"]
git-tree-sha1 = "85619a3f3e17bb4761fe1b1fd47f0e979f964d5b"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.0.2"
[[Parameters]]
deps = ["Markdown", "OrderedCollections", "REPL", "Test"]
git-tree-sha1 = "70bdbfb2bceabb15345c0b54be4544813b3444e4"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.10.3"
[[Pkg]]
deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
[[Polynomials]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "62142bd65d3f8aeb2226ec64dd8493349147df94"
uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
version = "0.5.2"
[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
[[Reexport]]
deps = ["Pkg"]
git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "0.2.0"
[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
[[SpecialFunctions]]
deps = ["BinDeps", "BinaryProvider", "Libdl", "Test"]
git-tree-sha1 = "0b45dc2e45ed77f445617b99ff2adf0f5b0f23ea"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "0.7.2"
[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
[[Test]]
deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[[URIParser]]
deps = ["Test", "Unicode"]
git-tree-sha1 = "6ddf8244220dfda2f17539fa8c9de20d6c575b69"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.0"
[[UUIDs]]
deps = ["Random"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
[[VersionParsing]]
deps = ["Compat"]
git-tree-sha1 = "c9d5aa108588b978bd859554660c8a5c4f2f7669"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.1.3"
name = "AeroAcoustics"
uuid = "03ede2b8-90d7-56cb-9af1-359a44c6802d"
keywords = ["AeroAcoustics"]
license = "MIT"
desc = "A package for post-processing af microphone array measurements"
authors = ["Oliver Lylloff <ollyl@dtu.dk>"]
version = "0.0.1"
[deps]
DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
[extras]
BinaryProvider = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
[targets]
test = ["Test", "Dates", "BinaryProvider", "HDF5"]
......@@ -8,5 +8,3 @@
A [julia](http://julialang.org) package for Aeroacoustics.
The roadmap and todos are tracked in the [roadmap meta-issue](https://gitlab.windenergy.dtu.dk/ollyl/AeroAcoustics.jl/issues/1).
Here is a small [demonstration](https://nbviewer.jupyter.org/urls/gitlab.windenergy.dtu.dk/ollyl/AeroAcoustics.jl/raw/master/examples/benchmark_b7_showcase.ipynb) of the current API (subject to change).
julia 0.6
julia 0.7
NLsolve
Distances
DSP
JuMP
HDF5
ProximalOperators
Statistics
LinearAlgebra
\ No newline at end of file
Pkg.update()
Pkg.clone(pwd(), "AeroAcoustics")
#Pkg.clone("https://gitlab.windenergy.dtu.dk/ollyl/AeroAcoustics.jl","AeroAcoustics")
Pkg.build("AeroAcoustics")
## This script prints coverage information.
Pkg.add("Coverage")
cd(Pkg.dir("AeroAcoustics"))
using Coverage
covered_lines, total_lines = get_summary(process_folder())
percentage = covered_lines / total_lines * 100
println("($(percentage)%) covered")
#__precompile__()
module AeroAcoustics
using NLsolve
using Distances
using DSP
using JuMP
#using SCS
using HDF5
using ProximalOperators
using LinearAlgebra # LinAlg in julia 1.0
using Statistics
using Parameters
import DSP
import Base.length,
Base.push!,
Base.Threads
export Constants,
export Environment,
CrossSpectralMatrix,
Environment,
SteeringMatrix,
cmf,
SPL,
octavebandlimits,
beamformer,
beamformer_corr,
beamformersetup,
parseHDF5data,
steeringvectors,
pointspreadfunction,
sourceintegration,
fista,
fistaprox!,
nonneg!,
csm,
# diagrm!,
fftnnls,
shear
beamforming
include("types.jl")
include("utils.jl")
include("steeringvectors.jl")
include("csm.jl")
include("beamformer.jl")
include("pointspreadfunction.jl")
#include("beamformer_corr.jl")
include("cmf.jl")
include("fista.jl")
include("fistaprox.jl")
include("fftnnls.jl")
include("beamforming.jl")
function SPL(p::Array{T}) where T <: Real
s = similar(p)
s[p.>0] = 10*log10.(p[p.>0]/4e-10)
s[p.<=0] = -350
return s
end
SPL(p::Number) = p > 0.0 ? 10*log10(p/4e-10) : -350.0
end
function beamformer(E::Environment{T},C::Constants,V::SteeringMatrix{T,S}) where {T <: AbstractFloat, S <: SteeringVectorType}
b = Array{T}(E.N,E.Nf)
vd = Array{Complex{T}}(E.M)
csm = Array{Complex{T}}(E.M,E.M)
# Compute beamforming
for j in 1:E.Nf
csm = E.CSM.csmReal[j,:,:] + im*E.CSM.csmImag[j,:,:]
for i in eachindex(E.D0)
vd = V.v[:,i,j]
b[i,j] = real(vd'*csm*vd)
end
end
return reshape(b,E.Nx,E.Ny,E.Nf)
end
function beamformer_old{T,C}(
Nx::Int64,
Ny::Int64,
X::Array{T,2},
Y::Array{T,2},
z0::T,
f::T,
rn::Array{T,2},
CSM::Array{C,2};
psf::Bool=false)
const M = size(rn,1) # Number of microphones
const c::Float64 = 343.0 # Speed of sound
const kw::T = 2pi*f/c # wavenumber
# CSM[eye(Bool,M)] = 0; # Naive diagonal removal
# Allocation of arrays
gj = Array{C}(M)
gjs = Array{C}(Nx,Ny,M)
b = Array{T}(Nx,Ny)
# Compute transfer functions
Threads.@threads for i in 1:Nx
for j in 1:Ny
r0::Float64 = sqrt(X[i,j]^2 + Y[i,j]^2 + z0^2)
#rsum = 0.0 # Type III Steering vector
for m in 1:M
rm::Float64 = sqrt((X[i,j]-rn[m,1])^2+(Y[i,j]-rn[m,2])^2 + z0^2)
#gj[m] = exp(-im*kw*(rm-r0)) # TYPE I Steering vector
gj[m] = (1/M)*(rm/r0)*exp(-im*kw*(rm-r0)) # TYPE II Steering vector
#gj[m] = (1/(r0*rm))*exp(-im*kw*(rm-r0)) # TYPE III Steering vector
#rsum += 1/rm^2
end
#gj *= 1/rsum # Type III Steering vector
gjs[i,j,:] = gj
b[i,j] = real(gj'*CSM*gj)
end
end
if psf
PSF = Array{T}(Nx,Ny)
grs = Array{C}(M)
midx::Int64 = 0
midy::Int64 = 0
if iseven(Nx)
midx = Nx/2
elseif isodd(Nx)
midx = round(Int64,Nx/2)+1
end
if iseven(Ny)
midy = Ny/2
elseif isodd(Ny)
midy = round(Int64,Ny/2)+1
end
grs = vec(gjs[midx,midy,:])
Threads.@threads for i in 1:Nx
for j in 1:Ny
#PSF[i,j] = abs(vec(gjs[i,j,:])'*grs*grs'*vec(gjs[i,j,:]))/M^2
PSF[i,j] = M^2*abs(vec(gjs[i,j,:])'*grs)^2 # Needs to multiply with (r0/rm)^2 to correct level
end
end
return b,gjs,PSF
else
return b,gjs
end
end
function beamformer_old{T,C}(
Nx::Int64,
Ny::Int64,
X::Array{T,2},
Y::Array{T,2},
z0::T,
f::Array{T,1},
rn::Array{T,2},
CSM::Array{C,3};
psf::Bool=false)
const M = size(rn,1) # Number of microphones
Nf = length(f)
b = Array{T}(Nx,Ny,Nf)
gjs = Array{C}(Nx,Ny,M,Nf)
if psf
PSF = Array{T}(Nx,Ny,Nf)
Threads.@threads for i in 1:Nf
b[:,:,i],gjs[:,:,:,i],PSF[:,:,i] = beamformer(Nx,Ny,X,Y,z0,f[i],rn,CSM[i,:,:];psf=true)
end
return b,gjs,PSF
else
Threads.@threads for i in 1:Nf
b[:,:,i],gjs[:,:,:,i] = beamformer(Nx,Ny,X,Y,z0,f[i],rn,CSM[i,:,:];psf=false)
end
return b,gjs
end
end
"""
beamforming(csm,v)
Calculate frequency-domain beamforming using cross-spectral matrix `csm` and steering vector `v`.
`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}
@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))
for j in axes(csm,3)
csmd = @view csm[:,:,j]
for i in axes(v,2)
vd = @view v[:,i,j]
b[i,j] = real(vd'*csmd*vd)
end
end
return b
end
function csm(t::AbstractArray{T};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n)) where T <: AbstractFloat
const nout = div((size(t,1) - n), n - noverlap)+1
const M = size(t,2)
const Nf = div(n,2)+1
const weight = sum(abs2,win)
const df = div(fs,n)
const nspec = div(n,2)
# 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}
ta = Array{Vector{T}}(undef,size(t,2))
for i in axes(t,2)
ta[i] = view(t,:,i)
end
return ta
end
Pxy = Array{Complex{T}}(Nf)
C = Array{Complex{T}}(Nf,M,M)
"""
csm(t;n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n))
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.
"""
function csm(t::AbstractArray{T};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n)) where T <: AbstractFloat
csm(flat_t(t);n=n,noverlap=noverlap,fs=fs,win=win)
end
function csm(t::Vector{Vector{T}};n=1024,noverlap=div(n,2),fs=1,win=DSP.hanning(n)) where T <: AbstractFloat
M = length(t)
Nf = div(n,2)+1
nout = div((length(t[1]) - n), n - noverlap)+1
weight = sum(abs2,win)
fc = [k*div(fs,n) for k in range(0,stop=noverlap)]
Pxy = Array{Complex{T}}(undef,Nf)
C = Array{Complex{T}}(undef,M,M,Nf)
ds = Array{Complex{T}}(undef,Nf,nout)
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
mean!(Pxy,conj!(DSP.stft(t[:,m], n, noverlap; fs=fs, window=win, onesided=true)).*DSP.stft(t[:,j], n, noverlap; fs=fs, window=win, onesided=true))
for ω in 1:Nf
C[ω,j,m] = Pxy[ω]
end
C[j,m,:] = mean!(Pxy,broadcast!(*,ds,Sc[m],S[j]))
#C[j,m,:] = mean!(Pxy,Sc[m].*S[j])
end
end
C .*= 2/n/weight
for ω in 1:Nf
C[ω,:,:] = Hermitian(C[ω,:,:],:L)
end
return CrossSpectralMatrix(real.(C),imag.(C),Array{T,1}([k*df for k in 0:nspec]),false)
#return C
end
function csm(::Type{T},filename::AbstractString) where T<:AbstractFloat
time_data = h5open(filename, "r") do file
read(file, "MicrophoneData/microphoneDataPa")
end
if T != eltype(time_data)
csm(Array{T,2}(time_data))
else
csm(time_data)
end
end
function csm(filename::AbstractString)
csm(Float64,filename)
end
# TODO: diagrm! Not finished!!
function diagrm!(csm::CrossSpectralMatrix{T}) where T <: AbstractFloat
Nf,M,M = size(csm.csmReal)
m = JuMP.Model(solver = SCS.SCSSolver())
JuMP.@variable(m,d[1:2M])
JuMP.@objective(m,Min,sum(d))
D = diagm(d)
for ω in 1:Nf
# TODO: Make imag_expand!() and call that here:
C = [csm.csmReal[ω,:,:] -csm.csmImag[ω,:,:]; csm.csmImag[ω,:,:] csm.csmReal[ω,:,:]]
JuMP.@SDconstraint(m,(C+D>=0))
JuMP.solve(m)
dopt = JuMP.getvalue(d)
for i in 1:M
csm.csmReal[ω,i,i] += dopt[i]
end
#csm.csmImag[ω,:,:] += diagm(dopt[M+1:2M])
C[:,:,ω] = Hermitian(C[:,:,ω],:L)
end
return CrossSpectralMatrix(csm.csmReal,csm.csmImag,csm.binCenterFrequenciesHz,true)
end
function diagrm_naive!(csm::CrossSpectralMatrix{T}) where T <: AbstractFloat
Nf,M,M = size(csm.csmReal)
for ω in 1:Nf
for i in 1:M
csm.csmReal[ω,i,i] = zero(T)
csm.csmImag[ω,i,i] = zero(T)
end