Merge branch 'unwrap' of https://github.com/platawiec/DSP.jl into unwrap

docs/make.jl

@@ -7,6 +7,7 @@ makedocs(
     pages = Any[
         "Contents" => "contents.md",
         "periodograms.md",
+        "estimation.md",
         "windows.md",
         "filters.md",
         "util.md",

docs/src/contents.md

@@ -6,6 +6,7 @@ DSP.jl provides a number of common DSP routines in Julia. So far, the following
 
 ```@contents
 Pages = ["periodograms.md",
+    "estimation.md",
     "windows.md",
     "filters.md",
     "util.md",

docs/src/estimation.md

@@ -0,0 +1,5 @@
+# `Estimation` - parametric estimation functions
+
+```@docs
+esprit
+```

src/DSP.jl

@@ -51,9 +51,10 @@ include("windows.jl")
 include("periodograms.jl")
 include("Filters/Filters.jl")
 include("lpc.jl")
+include("estimation.jl")
 
 using Reexport
-@reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Unwrap
+@reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Unwrap, .Estimation
 
 include("deprecated.jl")
 end

src/estimation.jl

@@ -0,0 +1,39 @@
+module Estimation
+import Base: *
+
+export esprit
+
+"""
+    esprit(x::AbstractArray, M::Integer, p::Integer, Fs::Real=1.0)
+
+ESPRIT [^Roy1986] algorithm for frequency estimation.
+Estimation of Signal Parameters via Rotational Invariance Techniques
+    
+Given length N signal "x" that is the sum of p sinusoids of unknown frequencies,
+estimate and return an array of the p frequencies.
+    
+# Arguments
+- `x::AbstractArray`: complex length N signal array
+- `M::Integer`: size of correlation matrix, must be <= N.
+      The signal subspace is computed from the SVD of an M x (N-M+1) signal matrix
+      formed from N-M+1 length-M shifts of the signal x in its columns.
+      For best performance for 1 sinusoid, use M = (N+1)/3 (according to van der Veen and Leus).
+      For faster execution (due to smaller SVD), use small M or small N-M
+- `p::Integer`: number of sinusoids to estimate.
+- `Fs::Float64`: sampling frequency, in Hz.
+ 
+# Returns
+length p real array of frequencies in units of Hz.
+
+[^Roy1986]: R Roy, A Paulraj and T Kailath, ESPRIT - A subspace approach to estimation of parameters of cisoids in noise, IEEE Trans. Acoustics, Speech, Signal Process., 34, 1340-1342 (1986). [url](http://ieeexplore.ieee.org/abstract/document/1164935/).
+"""
+function esprit(x::AbstractArray, M::Integer, p::Integer, Fs::Real=1.0)
+    count(v->v != 1, size(x)) <= 1 || error("`x` must be a 1D signal")
+    N = length(x)
+    X = x[ (1:M) .+ (0:N-M)' ]
+    U,s,V = svd(X)
+    D,_ = eig( U[1:end-1,1:p] \ U[2:end,1:p] )
+    angle.(D)*Fs/2π
+end
+
+end # end module definition

src/util.jl

@@ -329,5 +329,4 @@ function shiftin!(a::AbstractVector{T}, b::AbstractVector{T}) where T
     return a
 end
 
-
 end # end module definition

test/estimation.jl

@@ -0,0 +1,22 @@
+using DSP, Compat, Compat.Test
+
+@testset "esprit" begin
+    # create a sum of sinusoids in noise, and estimate their frequencies
+    Fs = 10000.0           # sampling frequency in Hz
+    duration = 1           # length of signal, in seconds
+    n = Int(Fs * duration) # number of samples
+    t = collect((1:n)/Fs)  # time vector
+    frequencies = [1000.0 1250.0]
+    amplitudes = [2.0 1.5]
+    phases = [0.7 -1.0]
+    x = exp.( 1im*2*π*t*frequencies .+ phases) * amplitudes'
+    noise = randn(n, 2)*[1;1im]
+    sigma = 0.1
+    noise *= sigma
+    x += noise
+    M = 300
+    p = 2                  # number of sinusoids to estimate
+    frequencies_estimated = esprit(x, M, p, Fs)
+    @test isapprox(frequencies', frequencies_estimated; atol = 1e-2)
+end
+

test/runtests.jl

@@ -2,7 +2,7 @@ using Compat, DSP, AbstractFFTs, FFTW, Compat.Test
 
 testfiles = [ "dsp.jl", "util.jl", "windows.jl", "filter_conversion.jl",
     "filter_design.jl", "filter_response.jl", "filt.jl", "filt_stream.jl",
-    "periodograms.jl", "resample.jl", "lpc.jl", "unwrap.jl"]
+    "periodograms.jl", "resample.jl", "lpc.jl", "estimation.jl", "unwrap.jl"]
 
 for testfile in testfiles
     eval(@__MODULE__, :(@testset $testfile begin include($testfile) end))

test/util.jl

@@ -13,7 +13,6 @@ using DSP, Compat, Compat.Test
     h_abs = abs.(h)
     h_angle = angle.(h)
     h_real = real.(h)
-
     #The real part should be equal to the original signals:
     @test h_real ≈ a