From 5bb29f8bb1179f7fa06f25a0e2bbd3b31b9aac5b Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Sat, 26 Aug 2023 23:35:21 +0800 Subject: [PATCH 1/5] doc(dsp): minor doc fix --- src/dsp.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dsp.jl b/src/dsp.jl index 91b2557..ce3453c 100644 --- a/src/dsp.jl +++ b/src/dsp.jl @@ -671,7 +671,7 @@ julia> y .+= 0.1 * randn(length(y)) julia> findsignal(x, y, 3; coarse=true) (time = Float32[0.000781, 0.001538, 0.003125], amplitude = ComplexF64[...]) julia> findsignal(x, y, 3) -(time = Float32[33, 64, 129], [0.000775, 0.001545, 0.003124], amplitude = ComplexF64[...]) +(time = Float32[0.000775, 0.001545, 0.003124], amplitude = ComplexF64[...]) ``` """ function findsignal(r, s, n=1; prominence=0.2, coarse=false) From da1974be438e2d8f22a9733cb94f4e3008aaaccd Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 28 Aug 2023 03:06:34 +0800 Subject: [PATCH 2/5] chore: add TODO to gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index eb05528..28fe48b 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ Manifest.toml docs/build .vscode/settings.json +TODO From 8ac3c104f68641a092ce4fe7247354e94c275bd2 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 28 Aug 2023 03:07:23 +0800 Subject: [PATCH 3/5] fix(dsp): improve robustness of mfilter for views and multichannel signals --- src/dsp.jl | 17 +++++++++++++++-- test/runtests.jl | 6 ++++++ 2 files changed, 21 insertions(+), 2 deletions(-) diff --git a/src/dsp.jl b/src/dsp.jl index ce3453c..671c69d 100644 --- a/src/dsp.jl +++ b/src/dsp.jl @@ -426,14 +426,27 @@ DSP.Filters.resample(x::SampledSignal, rate::Real, coef::Vector) = sresample(x, $(SIGNATURES) Matched filter looking for reference signal `r` in signal `s`. """ -function mfilter(r, s) +function mfilter(r, s::AbstractVector) issamerate(r, s) || throw(ArgumentError("signals `r` and `s` must have the same sampling rate")) - r̄, s̄ = promote(samples(r), samples(s)) + if eltype(r) == eltype(s) + r̄ = samples(r) + s̄ = samples(s) + else + T = promote_type(eltype(r), eltype(s)) + r̄ = convert(AbstractArray{T}, samples(r)) + s̄ = convert(AbstractArray{T}, samples(s)) + end f = conj.(reverse(r̄)) n = length(r) - 1 sfilt(f, padded(samerateas(s, s̄), (0, n)))[n+1:end] end +function mfilter(r, s::AbstractMatrix) + mapreduce(hcat, eachcol(s)) do x + mfilter(r, x) + end +end + """ $(SIGNATURES) Compute the inverse short time Fourier transform (ISTFT) of one-sided STFT coefficients `X` which is based diff --git a/test/runtests.jl b/test/runtests.jl index a4e1b14..c1686cf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -580,6 +580,12 @@ end x2 = mfilter(x, x1) @test !isanalytic(x2) + x = randn(rng, ComplexF64, 100) + x1 = signal(randn(rng, 1000, 2), 10kHz) + x2 = mfilter(x, x1) + @test isanalytic(x2) + @test nchannels(x2) == 2 + onesideds = [true, false] nfft = 256 windows = [nothing, rect, tukey(nfft, 0.5), hanning, hamming] From aede6201b951a114e2f4304fa25ce72707127b55 Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 28 Aug 2023 03:07:57 +0800 Subject: [PATCH 4/5] feat(dsp): add mfo option to findsignal() --- src/dsp.jl | 36 +++++++++++++++++++----------------- test/runtests.jl | 8 +++++--- 2 files changed, 24 insertions(+), 20 deletions(-) diff --git a/src/dsp.jl b/src/dsp.jl index 671c69d..a1bb22f 100644 --- a/src/dsp.jl +++ b/src/dsp.jl @@ -661,15 +661,16 @@ delay!(x, t::Units.Unitful.Time; kwargs...) = delay!(x, inseconds(t) * framerate """ $(SIGNATURES) -Finds up to `n` copies of reference signal `r` in signal `s`. The reference +Finds up to `n` strongest copies of reference signal `r` in signal `s`. The reference signal `r` should have a delta-like autocorrelation for this function to work well. If the keyword parameter `coarse` is set to `true`, approximate arrival times are computed based on a matched filter. If it is set to `false`, an iterative optimization is performed to find more accruate arrival times. -Returns named tuple `(time=t, amplitude=a)` where `t` is a vector of arrival -times and `a` is a vector of complex amplitudes of the arrivals. The arrivals -are sorted in ascending order of arrival times. +Returns named tuple `(time=t, amplitude=a, mfo=m)` where `t` is a vector of arrival +times, `a` is a vector of complex amplitudes of the arrivals, and `m` is the complex +matched filter output (if `mfo` is set to `true`). The arrivals are sorted in +ascending order of arrival times. # Examples: ```julia-repl @@ -682,29 +683,30 @@ julia> y4[513:512+length(x4)] += 0.6 * real(x4) # time 0.003125𝓈, index 129 julia> y = resample(y4, 1//4) julia> y .+= 0.1 * randn(length(y)) julia> findsignal(x, y, 3; coarse=true) -(time = Float32[0.000781, 0.001538, 0.003125], amplitude = ComplexF64[...]) -julia> findsignal(x, y, 3) -(time = Float32[0.000775, 0.001545, 0.003124], amplitude = ComplexF64[...]) +(time = [0.000781, 0.001538, 0.003125], amplitude = [...], mfo=[]) +julia> findsignal(x, y, 3; mfo=true) +(time = [0.000775, 0.001545, 0.003124], amplitude = [...], mfo=[...]) ``` """ -function findsignal(r, s, n=1; prominence=0.2, coarse=false) +function findsignal(r, s, n=1; prominence=0.2, coarse=false, mfo=false) # coarse arrival time estimation r = analytic(r) r = r / std(r) s = analytic(s) - mfo = mfilter(r, s) / length(r) - absmfo = abs.(samples(mfo)) - p, _ = findmaxima(absmfo) - peakproms!(p, absmfo; minprom=prominence*maximum(absmfo)) - length(p) > length(s)/10 && return (time=Float64[], amplitude=ComplexF64[]) - h = absmfo[p] + m = mfilter(r, s) / length(r) + T = eltype(m) + absm = abs.(samples(m)) + p, _ = findmaxima(absm) + peakproms!(p, absm; minprom=prominence*maximum(absm)) + length(p) > length(s)/10 && return (time=Float64[], amplitude=T[], mfo=(mfo ? m : T[])) + h = absm[p] ndx = sortperm(h; rev=true) length(ndx) > n && (ndx = ndx[1:n]) p = p[ndx] if coarse t = time(Float64.(p), s) ndx = sortperm(t) - return (time=t[ndx], amplitude=samples(mfo[p[ndx]])) + return (time=t[ndx], amplitude=samples(m[p[ndx]]), mfo=(mfo ? m : T[])) end # iterative fine arrival time estimation margin = 5 # arrival time may vary up to margin from coarse estimates @@ -724,7 +726,7 @@ function findsignal(r, s, n=1; prominence=0.2, coarse=false) end @view real(ifft(Z))[1:N] end - v0 = [p .- i; real.(mfo[p]); imag.(mfo[p])] + v0 = [p .- i; real.(m[p]); imag.(m[p])] optimize(v -> sum(abs2, reconstruct(v) .- s[i:i+N-1]), v0) end v = minimizer(soln) @@ -732,7 +734,7 @@ function findsignal(r, s, n=1; prominence=0.2, coarse=false) t = time(pp, s) a = complex.(v[length(p)+1:2*length(p)], v[2*length(p)+1:3*length(p)]) ndx = sortperm(t) - (time=t[ndx], amplitude=a[ndx]) + (time=t[ndx], amplitude=a[ndx], mfo=(mfo ? m : T[])) end """ diff --git a/test/runtests.jl b/test/runtests.jl index c1686cf..efeeaa9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -636,15 +636,17 @@ end @test energy(z) ≈ energy(x) y .+= 0.1 * randn(rng, length(y)) - t, a = findsignal(x, y, 3; coarse=false) + t, a, m = findsignal(x, y, 3; coarse=false) @test t ≈ [0.000775, 0.001545, 0.003124] atol=2e-6 @test real(a) / real(a[1]) ≈ [1.0, -0.8, 0.6] atol=1e-2 - t, a = findsignal(x, y, 3; coarse=true) + @test length(m) == 0 + t, a, m = findsignal(x, y, 3; coarse=true, mfo=true) @test t ≈ [0.000775, 0.001545, 0.003124] atol=1e-5 @test real(a) / real(a[1]) ≈ [1.0, -0.8, 0.6] atol=1e-2 + @test length(m) > 0 y = compose(real(x), time([32.75, 64.25, 129.0], x), [0.8, -0.7, 1.0]; duration=0.2) - t, a = findsignal(x, y, 3; coarse=false) + t, a, m = findsignal(x, y, 3; coarse=false) @test t ≈ [0.000775, 0.001545, 0.003124] atol=2e-6 @test real(a) / real(a[3]) ≈ [0.8, -0.7, 1.0] atol=1e-2 From 359412c9377cdb0243ec229d9634da36c3b7db9f Mon Sep 17 00:00:00 2001 From: Mandar Chitre Date: Mon, 28 Aug 2023 03:09:12 +0800 Subject: [PATCH 5/5] chore: version bump to 0.7.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6be3036..4e2735f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SignalAnalysis" uuid = "df1fea92-c066-49dd-8b36-eace3378ea47" authors = ["Mandar Chitre "] -version = "0.6.0" +version = "0.7.0" [deps] DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"