机器学习>> LPVSpectral.jl>> 返回
项目作者： baggepinnen

项目描述 ：
Least-squares (sparse) spectral estimation and (sparse) LPV spectral decomposition.
高级语言： Julia
项目地址: git://github.com/baggepinnen/LPVSpectral.jl.git
创建时间： 2017-01-23T12:23:37Z
项目社区：https://github.com/baggepinnen/LPVSpectral.jl
开源协议：Other
下载


LPVSpectral

CI
codecov

A toolbox for least-squares spectral estimation, sparse spectral estimation and Linear Parameter-Varying (LPV) spectral estimation. Contains an implementation of the spectral estimation method presented in
Bagge Carlson et al. “Linear Parameter-Varying Spectral Decomposition.” 2017 American Control Conference.

  1. @inproceedings{bagge2017spectral,
  2.   title        = {Linear Parameter-Varying Spectral Decomposition},
  3.   author       = {Bagge Carlson, Fredrik and Robertsson, Anders and Johansson, Rolf},
  4.   booktitle    = {2017 American Control Conference (ACC)},
  5.   year         = {2017},
  6. }

Extensions (sparse estimation methods) to the above article were developed in
Bagge Carlson, F., “Machine Learning and System Identification for Estimation in Physical Systems” (PhD Thesis 2018).

  1. @thesis{bagge2018,
  2.   title        = {Machine Learning and System Identification for Estimation in Physical Systems},
  3.   author       = {Bagge Carlson, Fredrik},
  4.   keyword      = {Machine Learning,System Identification,Robotics,Spectral estimation,Calibration,State estimation},
  5.   month        = {12},
  6.   type         = {PhD Thesis},
  7.   number       = {TFRT-1122},
  8.   institution  = {Dept. Automatic Control, Lund University, Sweden},
  9.   year         = {2018},
  10.   url          = {https://lup.lub.lu.se/search/publication/ffb8dc85-ce12-4f75-8f2b-0881e492f6c0},
  11. }

Installation

import Pkg; Pkg.add("LPVSpectral")

List of functions

This package provides tools for general least-squares spectral analysis, check out the functions

  1. ls_spectral             # Least-squares spectral analysis
  2. ls_sparse_spectral      # Least-squares sparse (L0) spectral analysis (uses ADMM)
  3. tls_spectral            # Total Least-squares spectral analysis
  4. ls_windowpsd            # Windowed Least-squares spectral analysis (sparse estimates available, see kwarg `estimator`)
  5. ls_windowcsd            # Windowed Least-squares cross-spectral density estimation (sparse estimates available, see kwarg `estimator`)
  6. ls_cohere               # Least-squares cross coherence estimation (sparse estimates available, see kwarg `estimator`)
  7. ls_spectral_lpv         # LPV spectral decomposition
  8. ls_sparse_spectral_lpv  # LPV spectral decomposition with group-lasso penalty on frequencies (uses ADMM)
  9. ls_windowpsd_lpv        # Windowed power spectral density estimation with LPV method
  11. mel                     # Compute Mel projection matrix
  12. melspectrogram          # Standard Mel spectrogram
  13. mfcc                    # Mel cepstrum spectrogram

The functions that estimate sparse spectra require the user to manually import using ProximalOperators.

All functions have docstrings available in the REPL. The general pattern is

  1. x,f = ls_XXX(y,t,f=default_freqs(t) [, W]; kwargs...)

where x are the complex Fourier coefficients and f are the frequency points. If no frequency vector is supplied, the default is to assume a sample time of 1 and use an equidistant grid from 0 to 0.5 of length(t)÷2.
W is an optional weight vector of length(y) for weighted least-squares estimation. Some methods accept keyword arguments, these methods are ls_windowpsd, ls_windowcsd, ls_cohere and the keywords and their defaults are
nw = 10, noverlap = -1, window_func=rect, estimator=ls_spectral.

Sparse spectral estimation

We provide a number of ways to estimate spare spectra. These functions require the user to manually load using ProximalOperators.

L₁ regularized spectral estimation

Minimize ||y-Ax||₂² + λ||x||₁ where x are the Fourier coefficients. Promotes a sparse spectrum

  1. x = ls_sparse_spectral(y,t,ω; proxg=NormL1(λ), tol=1e-9, printerval=1000, iters=30000, μ=0.000001)

L₀ regularized spectral estimation

Minimize ||y-Ax||₂² + λ||x||₀ where x are the Fourier coefficients. Promotes a sparse spectrum

  1. x = ls_sparse_spectral(y,t,ω; proxg=NormL0(λ), tol=1e-9, printerval=1000, iters=30000, μ=0.000001)

L₀ constrained spectral estimation

Minimize ||y-Ax||₂² s.t. ||x||₀ ≦ r where x are the Fourier coefficients. Enforces an r-sparse spectrum

  1. x = ls_sparse_spectral(y,t,ω; proxg=IndBallL0(r), tol=1e-9, printerval=1000, iters=30000, μ=0.000001)

Sparse LPV spectral estimation

See detailed example below and Bagge 2018.

  1. se = ls_sparse_spectral_lpv(Y,X,V,ω_test,Nv; λ = 0.1, normalize = normal, tol=1e-8, printerval=100, iters=6000)

LPV spectral estimation

We demonstrate the usage of the package with a simple example using simulated data, details can be found in the paper.

Signal generation

  1. using LPVSpectral, Plots, LaTeXStrings, DSP
  3. """
  4. `y,v,x = generate_signal(f,w,N)`
  6. `f` is a vector of functions `f(v)` that determine the functional dependence of the spectrum upon the velocity, one function for each frequency in `w`  both the amplitude and the phase are determined from these functions
  8. `w` is a vector of frequencies for which to estimate the spectrum
  10. `y,v,x` are output signal, sample points and scheduling variable respectively
  11. """
  12. function generate_signal(f,w,N, modphase=false)
  13.     x = sort(10rand(N)) # Sample points
  14.     v = range(0, stop=1, length=N) # Scheduling variable
  16.     # generate output signal
  17.     dependence_matrix = Float64[f[(i-1)%length(f)+1](v) for v in v, i in eachindex(w)] # N x nw
  18.     frequency_matrix  = [cos(w*x -0.5modphase*(dependence_matrix[i,j])) for (i,x) in enumerate(x), (j,w) in enumerate(w)] # N x nw
  19.     y = sum(dependence_matrix.*frequency_matrix,dims=2)[:] # Sum over all frequencies
  20.     y += 0.1randn(size(y))
  21.     y,v,x,frequency_matrix, dependence_matrix
  22. end
  24. N      = 500 # Number of training data points
  25. f      = [v->2v^2, v->2/(5v+1), v->3exp(-10*(v-0.5)^2),] # Functional dependences on the scheduling variable
  26. w      = 2π.*[2,10,20] # Frequency vector
  27. w_test = 2π.*(2:2:25) # Test Frequency vector, set w_test = w for a nice function visualization
  29. Y,V,X,frequency_matrix, dependence_matrix = generate_signal(f,w,N, true)

Signal analysis

We now make use of the spectral estimation method presented in the paper:

  1. # Options for spectral estimation
  2. λ      = 0.02 # Regularization parameter
  3. λs     = 1    # Regularization parameter group-lasso
  4. normal = true # Use normalized basis functions
  5. Nv     = 50   # Number of basis functions
  7. se  = ls_spectral_lpv(Y,X,V,w_test,Nv; λ = λ, normalize = normal) # Perform LPV spectral estimation
  8. ses = ls_sparse_spectral_lpv(Y,X,V,w_test,Nv; λ = λs, normalize = normal, tol=1e-8, printerval=100, iters=6000) # Same as above but with a group-lasso penalty on frequencies, promoting a solution with a sparse set of frequencies. Can be used to identify a sparse spectrum, i.e. to find w among w_test.

All that remains now is to visualize the result, along with the result of standard spectral estimation methods.

  1. plot(X,[Y V], linewidth=[1 2], lab=["\$y_t\$" "\$v_t\$"], xlabel=L"$x$ (sampling points)", title=L"Test signal $y_t$ and scheduling signal $v_t$", legend=true, xlims=(0,10), grid=false, c=[:cyan :blue])
  2. plot(se; normalization=:none, dims=2, l=:solid, c = [:red :green :blue], fillalpha=0.5, nMC = 5000, fillcolor=[RGBA(1,.5,.5,.5) RGBA(.5,1,.5,.5) RGBA(.5,.5,1,.5)], linewidth=2, bounds=true, lab=reshape(["Est. \$\\omega = $(round(w/π))\\pi \$" for w in w_test],1,:), phase = false)
  3. plot!(V,dependence_matrix, title=L"Functional dependencies $A(\omega,v)$", xlabel=L"$v$", ylabel=L"$A(\omega,v)$", c = [:red :green :blue], l=:dot, linewidth=2,lab=reshape(["True \$\\omega = $(round(w/π))\\pi\$" for w in w],1,:), grid=false)
  5. # Plot regular spectrum
  6. spectrum_lpv   = psd(se) # Calculate power spectral density
  7. spectrum_lpvs  = psd(ses) # Calculate sparse power spectral density
  8. fs             = N/(X[end]-X[1]) # This is the (approximate) sampling freqency of the generated signal
  9. spectrum_per   = DSP.periodogram(Y, fs=fs)
  10. spectrum_welch = DSP.welch_pgram(Y, fs=fs)
  11. plot(2π*collect(spectrum_per.freq), spectrum_per.power, lab="Periodogram", l=:path, m=:none, yscale=:log10, c=:cyan)
  12. plot!(2π*collect(spectrum_welch.freq), spectrum_welch.power, lab="Welch", l=:path, m=:none, yscale=:log10, linewidth=2, c=:blue)
  13. plot!(w_test,spectrum_lpv/fs, xlabel=L"$\omega$ [rad/s]", ylabel="Spectral density", ylims=(-Inf,Inf), grid=false, lab="LPV", l=:scatter, m=:o, yscale=:log10, c=:orange)
  14. plot!(w_test,spectrum_lpvs/fs, lab="Sparse LPV", l=:scatter, m=:o, c=:green)

window
window
window

When the three frequencies in w have been identified, w_test can be replaced by w for a nicer plot. As indicated by the last figure, the sparse estimate using group-lasso is better at identifying the three frequency components present (with a small bias in the estimation of the true frequencies).

Plotting

This package defines a recipe for plotting of periodogram types from DSP.jl. You can thus type

  1. using LPVSpectral, DSP, Plots
  2. plot(periodogram(y))
  3. plot(welch_pgram(y))
  4. plot(melspectrogram(y)) # melspectrogram, mel, mfcc are defined in this package