adds simpler getperiod function for skew correction

Author: ssfrr

Project.toml

@@ -6,8 +6,8 @@ Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
-Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SampledSignals = "bd7594eb-a658-542f-9e75-4c4d8908c167"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

src/MeasureIR.jl

@@ -3,18 +3,19 @@ module MeasureIR
 using Compat: @warn
 using Unitful: s, Hz, Time, Frequency, uconvert, NoUnits
 using SampledSignals: SampleBuf, samplerate
-using Roots: find_zero
 using DSP: DSP, gaussian, hanning, hilbert
-using DSP: FIRFilter, FIRWindow, filt, db2amp, xcorr, nextfastfft
+using DSP: FIRFilter, FIRWindow, filt, db2amp, pow2db, xcorr, nextfastfft
 using DSP: digitalfilter, Bandpass, Lowpass
+using DSP: resample
 using FFTW: rfft, irfft
 using LinearAlgebra: dot
 using Statistics: mean, quantile
-using LsqFit: curve_fit, stderror
+using Optim: NewtonTrustRegion, optimize, converged, minimizer
+using Printf
 
 export stimulus, analyze, noisefloor, prepadding, snr
 export expsweep, golay, mls, rpms, impulse
-export pilot, findpilot, pilotsync
+export pilot, findpilot, pilotsync, getperiod
 
 """
 Abstract supertype for impulse response measurements. Subtypes should define a
@@ -84,6 +85,7 @@ include("expsweep.jl")
 include("mls.jl")
 include("rpms.jl")
 include("pilot.jl")
+include("period.jl")
 
 # workaround needed because DSP.jl doesn't handle SampleBufs
 # and Float32s well. We dispatch to an internal _analyze to avoid method

src/period.jl

@@ -0,0 +1,58 @@
+"""
+Estimate the period of `x`, from an initial guess `L`. The result will be the
+peak of the continuous cross-correlation between `0.5L` and `1.5L`, i.e. the
+period is not constrained to be an integer.
+
+The algorithm is to first compute a discrete cross-correlation and oversample by
+4x to account for inter-sample peaks. We use the peak of the oversampled
+cross-correlation as a starting point for a numerical optimization of the
+average-squared-distance function in the frequency domain.
+"""
+function getperiod(x::AbstractVector, L)
+    maxD = L÷2
+
+    # create two signals that are offset by L (our expected period length), and
+    # zero-pad by to do a linear (rather than circular) convolution. We set up
+    # the padding to pre-shift so the cross-correlation starts start at their
+    # left-most lag (`-maxD`), and the `maxD`th sample will actually be L lag)
+    # we also add extra padding to account for time aliasing from the fractional
+    # delay. Otherwise it can wrap around. 512 samples is enough for a sinc
+    # function to decay more than 60dB
+
+    # note this could be problematic if the signal is temporally compact, e.g.
+    # an impulse train. Innacuracy in `L` could chop off an impulse from one or
+    # the other of these.
+    x1 = [zeros(maxD); x[L+1:end]; zeros(512)]
+    x2 = [x[1:(end-L)]; zeros(maxD); zeros(512)]
+
+    m = length(x1)
+    X1 = rfft(x1)
+    X2 = rfft(x2)
+    M = length(X1)
+
+    # do a discrete cross-correlation and upsample. The maximum here should be
+    # a pretty good estimate of the true peak. Should be valid from 1:8L, which
+    # represents lags from -L:L (upsampled by 4)
+    xc = resample(irfft(X1 .* conj.(X2), m)[1:2maxD+1], 4)
+    D_est1 = (findmax(xc)[2]-1)/4-maxD
+
+    # note if `m` were guaranteed to be even the upper bound here would be π
+    ω = range(0, (M-1)*2π/m, length=M)
+
+    # error function for  frequency-domain delay estimation. This can be
+    # interpreted as a frequency-domain ASDF (average-squared distance function)
+    # as a function of the continuous delay
+    # err(D) = sum(abs2(X1[i] - X2[i] * exp(-im*ω[i]*(D+maxD)))
+    #              for i in 1:M)*(m-L)/(m-L-abs(maxD-D))
+    err(D) = sum(abs2(X1[i] - X2[i] * exp(-im*ω[i]*(D+maxD)))
+                 for i in 1:M)
+
+    # opt = optimize(D->err(first(D)), [D_est1], NewtonTrustRegion(), autodiff=:forward)
+    opt = optimize(err, D_est1-0.5, D_est1+0.5)
+    if !converged(opt)
+        throw(ErrorException("Delay estimation failed to converge"))
+    end
+
+    # return the delay
+    L+minimizer(opt)[1]
+end

src/pilot.jl

@@ -14,48 +14,28 @@ function pilotfilt(sig, ω)
 end
 
 """
-Compute reasonable onset and offset thresholds to find a pulse of length L
-within `sig`.
-"""
-function pulsethreshs(sig, L)
-    if length(sig) < 4L
-        @warn "`pulsethreshs` assumes `length(sig) > 4L` to estimate noise statistics"
-    end
-    lfrac = L/length(sig)
-    qs = quantile(sig, [0.75, 1-lfrac/2])
-    qs[1] .+ (0.55, 0.45) .* (qs[2]-qs[1])
-end
-
-"""
-Find a pilot tone in the given signal. Note that the total signal length should
-be at least 4 times longer than the expected length of the tone, so we can
-estimate noise statistics.
+Find a pilot tone in the given signal.
 """
 function findpilot(sig, freq, len; samplerate=1)
+    # TODO: this isn't robust to a wideband energy increase. It looks for an
+    # energy pulse in the target band, but that could also get triggered if
+    # everything just gets louder
     ω = freqnorm(freq, samplerate)
     L = durnorm(len, samplerate)
     filtered = pilotfilt(sig, ω)
+
     env = abs2.(hilbert(filtered))
-    onthresh, offthresh = pulsethreshs(env, L)
-    # lpf = gaussian(round(Int, 0.25*L), 0.15)
-    lpf = digitalfilter(Lowpass(1/4L),
-                        FIRWindow(;transitionwidth=30/L))
-    smoothed = filt(lpf, env)
-    onset = findfirst((x->x>onthresh), smoothed)
-    onset === nothing && return nothing, nothing
-    offset = onset+findfirst(x->x<offthresh, smoothed[onset+1:end])
-    filtdelay = length(lpf)÷2
-    onset -= filtdelay
-    offset === nothing && return onset, nothing
-    offset -= filtdelay
-    if !isapprox(offset-onset, L; rtol=0.5)
-        @warn "Found pilot with length $(offset-onset) instead of $L"
-    end
+    _, onset = findmax(xcorr(env, ones(L); padmode=:none))
+    onset -= L
 
-    # only pull the middle 80% of the tone to account for any other onset/offset
-    # irregularities
-    slack = 0.05(offset-onset)
-    round.(Int, (onset+slack, offset-slack))
+    prerollpower = sum(x->x^2, filtered[1:onset-1]) / (onset-1)
+    pilotpower = sum(x->x^2, filtered[(0:L-1).+onset]) / L
+
+    snr = pow2db(pilotpower/prerollpower)
+    if snr < 6
+        @warn @sprintf "Narrowband Pilot SNR only %.2f dB" snr
+    end
+    onset
 end
 
 """

test/Project.toml

@@ -1,8 +1,9 @@
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
-MeasureIR = "1d92ad0d-f520-5482-bb8f-e5ddddeec5c1"
+MeasureIR = "4e68e94c-cd2d-5c3b-bae2-1718e3447baa"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 SampledSignals = "bd7594eb-a658-542f-9e75-4c4d8908c167"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

test/runtests.jl

@@ -1,10 +1,55 @@
 using MeasureIR
 using Unitful: s, kHz
 using SampledSignals: SampleBuf, samplerate
-# using Suppressor
+using Suppressor
 using DSP
 using Test
 
+"""
+Simulate running the given stimulus through a real system, including a linear
+IR, nonlinear distortion, noise, and clock skew.
+
+If `dist` is given, it should be a multiplier that scales the stimulus before
+passing into `tanh`. The distorted signal is scaled to have the same energy as
+the original. A multiplier of 1 gives pretty mild distortion, and it goes up
+from there.
+
+Returns the simulated response, as well as the IR used.
+"""
+function simIR(stim; ir=nothing, dist=nothing, noisepow=-Inf, skew=1.0)
+    if dist !== nothing
+        stimenergy = sum(x->x^2, stim)
+        stim .= tanh.(stim)
+        distenergy = sum(x->x^2, stim)
+
+        stim .*= sqrt(sigenergy/distenergy)
+    end
+    if ir === nothing
+        ir = [zeros(100); 1; randn.() .* exp.(range(-1.5, -8, length=7500))]
+    end
+    resp = conv(stim, ir)
+    if noisepow > -Inf
+        resp .+= db2amp(noisepow) .* randn.()
+    end
+    if skew != 1.0
+        resp = resample(resp, skew)
+    end
+
+    resp, ir
+end
+
+"""
+Generates a broadband random signal of length `N` that's periodic with period
+`P`. `P` does not need to be an integer.
+"""
+function randperiodic(N, P)
+    f = 1/P # fundamental freq
+    # we add eps() here so we don't get hit by floating-point error when f
+    # divides evenly into 1
+    H = Int(fld(0.5+eps(), f)) # number of harmonics we can fit below nyquist
+    sum(cos.(2π.*(h.*f.*(0:N-1).+rand()))./H for h in 1:H)
+end
+
 function testmeasure(measfunc)
     meas = measfunc(8192)
     stim = stimulus(meas)
@@ -103,10 +148,10 @@ end
     meas = golay(L)
     ir = randn(L÷2) .* exp.(linspace(0,-8, L÷2))
     stim = stimulus(meas)
-    out = @color_output false @capture_err analyze(meas, conv(stim, ir))
+    out = @capture_err analyze(meas, conv(stim, ir))
     @test out == ""
     ir = randn(3L÷2) .* exp.(linspace(0,-8, 3L÷2))
-    out = @color_output false @capture_err analyze(meas, conv(stim, ir))
+    out = @capture_err analyze(meas, conv(stim, ir))
     @test out == string("Warning: nonsilent samples past end of analysis ",
                         "window. Check your test signal is long enough for ",
                         "the response\n")
@@ -117,17 +162,17 @@ end
 #     L = 4096
 #     meas = expsweep(L)
 #     stim = stimulus(meas)
-#     out = @color_output false @capture_err analyze(meas, stim)
+#     out = @capture_err analyze(meas, stim)
 #     @test out == ""
-#     out = @color_output false @capture_err analyze(meas, tanh.(stim*2))
+#     out = @capture_err analyze(meas, tanh.(stim*2))
 #     @test out == "Warning: Energy in noncausal IR is above noise floor. Check for nonlinearity\n"
 #
 #     # now try with some noise and a more realistic IR
 #     testir = randn(512) .* exp.(linspace(0,-8, 512)) / 2
 #     resp = conv(stim, testir)[1:length(stim)]
-#     out = @color_output false @capture_err analyze(meas, stim.+randn.().*0.1)
+#     out = @capture_err analyze(meas, stim.+randn.().*0.1)
 #     @test out == ""
-#     out = @color_output false @capture_err analyze(meas, tanh.(stim*2).+randn.().*0.1)
+#     out = @capture_err analyze(meas, tanh.(stim*2).+randn.().*0.1)
 #     @test out == "Warning: Energy in noncausal IR is above noise floor. Check for nonlinearity\n"
 # end
 
@@ -189,34 +234,46 @@ end
 # plt = plot(welch_pgram(sig[(1:10000) .+ 50000], 256).power)
 @testset "findpilot" begin
     sr = 48kHz
+    dur = 1s
     @testset "f=$f" for f in 0.99kHz:0.005kHz:1.01kHz
         sig = [zeros(50000); pilot(dur, f, samplerate=sr); zeros(48000*10)]
         sig .+= 8 .* randn.() # SNR of -21.1dB - not bad!
-        onset, offset = findpilot(sig, 1kHz, 1s; samplerate=sr)
-        # make sure we don't over-shoot
-        @test offset-onset <= 48000
-        # undershooting by a bit is OK
-        @test offset-onset > 48000*0.7
-        @test onset >= 50000
-        @test onset <= 50000+48000*0.3
+        onset = findpilot(sig, 1kHz, 1s; samplerate=sr)
+        @test isapprox(onset, 50001, rtol=0.1)
     end
 end
 
-# TODO: for some reason this is failing when skew == 1
 @testset "pilotsync" begin
-    L = 5s
+    dur = 5s
     sr = 48kHz
     f = 1kHz
+    L = 5*48000
     @testset "skew=$skew" for skew in 0.99:0.005:1.01
-        p = pilot(L, f; samplerate=sr)
+        p = pilot(dur, f; samplerate=sr)
         sig = [zeros(50000);
                isapprox(skew, 1) ? p : resample(p, skew);
                zeros(48000*15)]
         sig .+= randn.() ./ sqrt(2) # SNR of 0dB
-        onset, offset = findpilot(sig, f, L; samplerate=sr)
-        @test isapprox(pilotsync(sig[onset:offset], f; samplerate=sr),
+        onset = findpilot(sig, f, dur; samplerate=sr)
+        @test isapprox(pilotsync(sig[(0:L-1) .+ onset], f; samplerate=sr),
                        1/skew, rtol=5e-7)
     end
 end
 
+@testset "getperiod" begin
+    N = 65535
+    meas = mls(N, 20)
+
+    # test at a variety of skews in ppm
+    skews = 1 .+ (-400:80:400) ./ 1e6
+
+    @testset "skew=$skew" for skew in skews
+        resp, ir = simIR(stimulus(meas), skew=skew, noisepow=-8.2+20) # -20dB SNR
+        period = getperiod(resp, N)
+        err = (period-skew*N)/period
+        @show skew, err
+        @test isapprox(period, skew*N, rtol=1e-7)
+    end
+end
+
 end # @testset MeasureIR