feat(dsp): add mfo option to findsignal()

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,15 +726,15 @@ 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)
   pp = v[1:length(p)] .+ i
   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
 
 """

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