Add frequency estimators by Jacobsen and Quinn

src/estimation.jl

@@ -44,29 +44,19 @@ function esprit(x::AbstractArray, M::Integer, p::Integer, Fs::Real=1.0)
 end
 
 """
-    jacobsen(x::Vector, Fs::Real = 1.0)
+    jacobsen(x::AbstractVector, Fs::Real = 1.0)
 
-Estimate the largest frequency (in Hz) in the signal `x` using Jacobsen's
-algorithm [^Jacobsen2007].
+Estimate the largest frequency in the signal `x` using Jacobsen's algorithm
+[^Jacobsen2007]. Argument `Fs` is the sampling frequency.  All frequencies are
+expressed in Hz.
+
+If the sampling frequency `Fs` is not provided, then it is assumed that `Fs =
+1.0`.
 
 [^Jacobsen2007]: E Jacobsen and P Kootsookos, "Fast, Accurate Frequency Estimators", Chapter
 10 in "Streamlining Digital Signal Processing", edited by R. Lyons, 2007, IEEE Press.
 """
-function jacobsen(x::Vector{<:Real}, Fs::Real = 1.0)
-    N = length(x)
-    X = rfft(x)
-    k = argmax(abs.(X))  # index of DFT peak
-    if ((k != div(N,2)+1) && (k != 1))  # avoid out-of-bounds indexing
-        # Jacoben's formula
-        δ = -real((X[k+1] - X[k-1]) / (2X[k] - X[k-1] - X[k+1]) )
-    else
-        δ = 0.0
-    end
-
-    return (k - 1 + δ)*Fs/N
-end
-
-function jacobsen(x::Vector{<:Complex}, Fs::Real = 1.0)
+function jacobsen(x::AbstractVector, Fs::Real = 1.0)
     N = length(x)
     X = fftshift(fft(x))
     k = argmax(abs.(X))  # index of DFT peak
@@ -99,9 +89,8 @@ is assumed that `Fs = 1.0`.
 
 The following keyword arguments control the algorithm's behavior:
 
-- `tol`: () the algorithm stops when the absolut value of the
-  difference between two consecutive estimates is less than `tol`. Defaults to
-  `1e-6`.
+- `tol`: the algorithm stops when the absolut value of the difference between
+   two consecutive estimates is less than `tol`. Defaults to `1e-6`.
 - `maxiters`: the maximum number of iterations to run. Defaults to `20`.
 
 Returns a tuple `(estimate, reachedmaxiters)`, where `estimate` is the
@@ -125,7 +114,7 @@ quinn(x, Fs ; kwargs...) = quinn(x, jacobsen(x, Fs), Fs ; kwargs...)
 
 function quinn(x::Vector{<:Real}, f0::Real, Fs::Real ; tol = 1e-6, maxiters = 20)
     fₙ = Fs/2
-    T = length(x)
+    N = length(x)
 
     # Run a quick estimate of largest sinusoid in x
     ω̂ = π*f0/fₙ
@@ -137,18 +126,18 @@ function quinn(x::Vector{<:Real}, f0::Real, Fs::Real ; tol = 1e-6, maxiters = 20
     α = 2cos(ω̂)
 
     # iteration
-    ξ = zeros(eltype(x), T)
+    ξ = zeros(eltype(x), N)
     β = zero(eltype(x))
     iter = 0
     @inbounds while iter < maxiters
         iter += 1
         ξ[1] = x[1]
         ξ[2] = x[2] + α*ξ[1]
-        for t in 3:T
+        for t in 3:N
             ξ[t] = x[t] + α*ξ[t-1] - ξ[t-2]
         end
         β = ξ[2]/ξ[1]
-        for t = 3:T
+        for t = 3:N
             β += (ξ[t]+ξ[t-2])*ξ[t-1]
         end
         β = β/(ξ[1:end-1]'*ξ[1:end-1])
@@ -161,26 +150,26 @@ end
 
 function quinn(x::Vector{<:Complex}, f0::Real, Fs::Real ; tol = 1e-6, maxiters = 20)
     fₙ = Fs/2
-    T = length(x)
+    N = length(x)
 
     ω̂ = π*f0/fₙ
 
     # Remove any DC term in x
     x .= x .- mean(x)
 
     # iteration
-    ξ = zeros(eltype(x), T)
+    ξ = zeros(eltype(x), N)
     iter = 0
     @inbounds while iter < maxiters
         iter += 1
         # step 2
         ξ[1] = x[1]
-        for t in 2:T
+        for t in 2:N
             ξ[t] = x[t] + exp(complex(0,ω̂))*ξ[t-1]
         end
         # step 3
         S = let s = 0.0
-                for t=2:T
+                for t=2:N
                     s += x[t]*conj(ξ[t-1])
                 end
                 s

test/estimation.jl

@@ -21,25 +21,9 @@ using DSP, Test
 end
 
 @testset "jacobsen" begin
-    # real input: test at two arbitrary frequencies
+    # test at two arbitrary frequencies
     fs = 100
     t = range(0, 5, step = 1/fs)
-    fr = 28.3
-    sr = cos.(2π*fr*t .+ π/4.2)
-    f_est_real = jacobsen(sr, fs)
-    @test isapprox(f_est_real, fr, atol = 1e-5)
-    fr = 12.45
-    sr = sin.(2π*fr*t .+ 3π/2.2)
-    f_est_real = jacobsen(sr, fs)
-    @test isapprox(f_est_real, fr, atol = 1e-5)
-    # test at higher extreme of DFT
-    fr = 49.9002
-    sr = cos.(2π*fr*t)
-    f_est_real = jacobsen(sr, fs)
-    @test isapprox(f_est_real, fr, atol = 1e-5)
-    # test at lower extreme of DFT
-    @test jacobsen(zeros(10)) == 0.0
-    # complex input: test at two arbitrary frequencies
     fc = -40.3
     sc = cis.(2π*fc*t .+ π/1.4)
     f_est_complex = jacobsen(sc, fs)
@@ -61,21 +45,20 @@ end
 end
 
 @testset "quinn" begin
-    # real input
+    ### real input
     fs = 100
     t = range(0, 5, step = 1/fs)
     fr = 28.3
     sr = cos.(2π*fr*t .+ π/4.2)
-    ### real input
     (f_est_real, maxiter) = quinn(sr, 50, fs)
     @test maxiter == false
     @test isapprox(f_est_real, fr, atol = 1e-3)
     # use default initial guess
-    (f_est_real, maxiter) = quinn(sr, fs) # assumes f0 = 0.0
+    (f_est_real, maxiter) = quinn(sr, fs) # initial guess given by Jacobsen
     @test maxiter == false
     @test isapprox(f_est_real, fr, atol = 1e-3)
     # use default fs
-    (f_est_real, maxiter) = quinn(sr) # assumes fs = 1.0, f0 = 0.0
+    (f_est_real, maxiter) = quinn(sr) # fs = 1.0, initial guess given by Jacobsen
     @test maxiter == false
     @test isapprox(f_est_real, fr/fs, atol = 1e-3)
     ### complex input
@@ -85,11 +68,11 @@ end
     @test maxiter == false
     @test isapprox(f_est_real, fc, atol = 1e-3)
     # use default initial guess
-    (f_est_real, maxiter) = quinn(sc, fs) # assumes f0 = 0.0
+    (f_est_real, maxiter) = quinn(sc, fs) # initial guess given by Jacobsen
     @test maxiter == false
     @test isapprox(f_est_real, fc, atol = 1e-3)
     # use default fs
-    (f_est_real, maxiter) = quinn(sc) # assumes fs = 1.0, f0 = 0.0
+    (f_est_real, maxiter) = quinn(sc) # fs = 1.0, initial guess by Jacobsen
     @test maxiter == false
     @test isapprox(f_est_real, fc/fs, atol = 1e-3)
 end