From a5116cbcddb9a12c8eb3ed760d52c56115bb3e11 Mon Sep 17 00:00:00 2001 From: Tom Krauss Date: Mon, 23 Apr 2018 06:16:21 -0500 Subject: [PATCH] Add esprit function to util (#193) * Add esprit function to util esprit finds the frequencies of sinusoids. It can also be used for direction finding in array processing and time-delay estimation in the frequency domain. * add test for esprit * fix bugs in esprit test * increase tolerance in esprit test * incorporate all review comments - added a reference to ESPRIT paper - changed types of function - moved "esprit" out of util into a new "estimation.jl" module - Keeping the "X" matrix implementation as is where a large index matrix is formed first. In my simple tests this was faster although it used more memory. A double for-loop implementation may still be better but is left "for futher study." * get rid of duplicate test for esprit. Add the esprit doc string to the docs. --- docs/make.jl | 1 + docs/src/contents.md | 1 + docs/src/estimation.md | 5 +++++ src/DSP.jl | 3 ++- src/estimation.jl | 39 +++++++++++++++++++++++++++++++++++++++ src/util.jl | 1 - test/estimation.jl | 22 ++++++++++++++++++++++ test/runtests.jl | 2 +- test/util.jl | 1 - 9 files changed, 71 insertions(+), 4 deletions(-) create mode 100644 docs/src/estimation.md create mode 100644 src/estimation.jl create mode 100644 test/estimation.jl diff --git a/docs/make.jl b/docs/make.jl index 772126a8c..2113b87c5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,6 +7,7 @@ makedocs( pages = Any[ "Contents" => "contents.md", "periodograms.md", + "estimation.md", "windows.md", "filters.md", "util.md", diff --git a/docs/src/contents.md b/docs/src/contents.md index 0ee0d263a..9c598d746 100644 --- a/docs/src/contents.md +++ b/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", diff --git a/docs/src/estimation.md b/docs/src/estimation.md new file mode 100644 index 000000000..767bc8b73 --- /dev/null +++ b/docs/src/estimation.md @@ -0,0 +1,5 @@ +# `Estimation` - parametric estimation functions + +```@docs +esprit +``` diff --git a/src/DSP.jl b/src/DSP.jl index 2e020c7e9..068d0d74d 100644 --- a/src/DSP.jl +++ b/src/DSP.jl @@ -50,9 +50,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 +@reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Estimation include("deprecated.jl") end diff --git a/src/estimation.jl b/src/estimation.jl new file mode 100644 index 000000000..2d585a8cd --- /dev/null +++ b/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 diff --git a/src/util.jl b/src/util.jl index 2082e49d3..002865071 100644 --- a/src/util.jl +++ b/src/util.jl @@ -371,5 +371,4 @@ function shiftin!(a::AbstractVector{T}, b::AbstractVector{T}) where T return a end - end # end module definition diff --git a/test/estimation.jl b/test/estimation.jl new file mode 100644 index 000000000..4a3811b53 --- /dev/null +++ b/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 + diff --git a/test/runtests.jl b/test/runtests.jl index b536c968f..92624dcf4 100644 --- a/test/runtests.jl +++ b/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"] + "periodograms.jl", "resample.jl", "lpc.jl", "estimation.jl"] for testfile in testfiles eval(@__MODULE__, :(@testset $testfile begin include($testfile) end)) diff --git a/test/util.jl b/test/util.jl index f68971cc3..0b4754895 100644 --- a/test/util.jl +++ b/test/util.jl @@ -41,7 +41,6 @@ end 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