Skip to content

Commit

Permalink
fix(dsp): improve type stability of findsignal()
Browse files Browse the repository at this point in the history
  • Loading branch information
mchitre committed Aug 19, 2023
1 parent b2f5958 commit 3f10bed
Showing 1 changed file with 14 additions and 13 deletions.
27 changes: 14 additions & 13 deletions src/dsp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -682,35 +682,36 @@ function findsignal(r, s, n=1; prominance=0.2, coarse=false)
absmfo = abs.(samples(mfo))
p, _ = findmaxima(absmfo)
peakproms!(p, absmfo; minprom=prominance*maximum(absmfo))
length(p) > length(s)/10 && return nothing
length(p) > length(s)/10 && return (time=Float64[], amplitude=ComplexF64[])
h = absmfo[p]
ndx = sortperm(h; rev=true)
length(ndx) > n && (ndx = ndx[1:n])
p = p[ndx]
if coarse
t = time(p, s)
return (time=t, amplitude=mfo[p])
t = time(Float64.(p), s)
return (time=t, amplitude=samples(mfo[p]))
end
# iterative fine arrival time estimation
margin = 5 # arrival time may vary up to margin from coarse estimates
i = minimum(p)
i::Int = minimum(p)
n = maximum(p) - i + length(r) + 2 * margin
n = nextfastfft(n)
i = max(1, i - margin)
N = n
i + N - 1 > length(s) && (N = length(s) - i + 1)
X = fft(vcat(samples(r), zeros(n-length(r))))
f = fftfreq(n)
function reconstruct(v)
ii = @view v[1:length(p)]
aa = @views complex.(v[length(p)+1:2*length(p)], v[2*length(p)+1:3*length(p)])
Z = mapreduce(+, zip(ii, aa)) do (i, a)
a .* X .* cis.(-2π .* i .* f)
soln = let p=p, f=fftfreq(n)
function reconstruct(v)
ii = @view v[1:length(p)]
aa = @views complex.(v[length(p)+1:2*length(p)], v[2*length(p)+1:3*length(p)])
Z = mapreduce(+, zip(ii, aa)) do (i, a)
a .* X .* cis.(-2π .* i .* f)
end
@view real(ifft(Z))[1:N]
end
@view real(ifft(Z))[1:N]
v0 = [p .- i; real.(mfo[p]); imag.(mfo[p])]
optimize(v -> sum(abs2, reconstruct(v) .- s[i:i+N-1]), v0)
end
v0 = [p .- i; real.(mfo[p]); imag.(mfo[p])]
soln = optimize(v -> sum(abs2, reconstruct(v) .- s[i:i+N-1]), v0)
v = minimizer(soln)
pp = v[1:length(p)] .+ i
t = time(pp, s)
Expand Down

0 comments on commit 3f10bed

Please sign in to comment.