Edit on GitHub

Models and root manipulations

Overview

Most distances available in this package operate on linear models estimated from a time-domain signal. This package supports two kind of LTI models, AR and ARMA. AR represents a model with only poles whereas ARMA has zeros as well. These types are subtypes of ControlSystemsBase.LTISystem, so many of the functions from the ControlSystemsBase.jl toolbox work on these models as well. When acting like a ControlSystemsBase.LTISystem, the default is to use the continuous-time representation of the model. The discrete-time representation can be obtained by tf(m, 1) where 1 is the sample time. More on the time-domain representation below.

Note

This package makes the assumption that the sample time is 1 everywhere. When an AbstractModel is constructed, one must thus take care to rescale the frequency axis accordingly if this does not hold. If the discrete-time representation is never used, this is of no concern.

To fit a model to data, one first has to specify a FitMethod, the options are

IRLS
LS
PLR
TLS
TimeWindow

For example, to estimate an AR model of order 6 using least-squares, we can do the following

julia> data = randn(1000);
julia> fitmethod = LS(na=6)LS(6, 0.01)
julia> model = fitmethod(data)AR{coeff type: Float64, root type: ComplexF64}(
b: 103.84062840482422
Cont. poles: 6-element ContinuousRoots{ComplexF64, Vector{ComplexF64}}:
 -0.7443030205850414 - 2.4899449018437965im
 -0.7752684097766342 - 1.6738261171647648im
 -0.5688112262426465 - 0.6720601480368995im
 -0.5688112262426465 + 0.6720601480368995im
 -0.7752684097766342 + 1.6738261171647648im
 -0.7443030205850414 + 2.4899449018437965im
Abs disc: [0.4750652953096256, 0.4605801401939561, 0.56619812025405, 0.56619812025405, 0.4605801401939561, 0.4750652953096256]
Cont. coeffs: 7-element Vector{Float64}:
  1.0
  4.176765313208644
 16.69729426336235
 32.073450812713986
 50.319854279001234
 38.18893433645754
 17.8154619523519
Spectral energy: 0.9850350554830362
julia> change_precision(Float32, model) # Can be useful to reduce the computational cost of some distancesAR{coeff type: Float32, root type: ComplexF32}(
b: 103.84062840482422
Cont. poles: 6-element ContinuousRoots{ComplexF32, Vector{ComplexF32}}:
 -0.74430305f0 - 2.489945f0im
 -0.77526844f0 - 1.6738261f0im
 -0.56881124f0 - 0.67206013f0im
 -0.56881124f0 + 0.67206013f0im
 -0.77526844f0 + 1.6738261f0im
 -0.74430305f0 + 2.489945f0im
Abs disc: Float32[0.4750653, 0.46058014, 0.5661981, 0.5661981, 0.46058014, 0.4750653]
Cont. coeffs: 7-element Vector{Float32}:
  1.0
  4.1767654
 16.697294
 32.073452
 50.319855
 38.188934
 17.815462
Spectral energy: 0.98503506
julia> pzmap(model)Plot{Plots.PlotlyBackend() n=1}

Time domain

This package allows you to view a model through two different lenses: as a continuous-time model that models the differential properties of the signal, or as a discrete-time model that models the difference properties of the signal. Signals are inevetably sampled before the computer interacts with them, and are thus natively in the discrete domain. Theory, on the other hand, is slightly more intuitive in the continuous time domain. The two domains are related by the conformal mapping $p_c = \log(p_d)$ where $p$ denotes a pole of a transfer function and subscripts $c,d$ denote the continuous and discrete domains respectively. When creating a distance, the default domain is Continuous. Some functions require you to be explicit regarding which domain you have in mind, such as when creating models from vectors or when asking for the roots/poles of a model.

Sometimes you may get a message saying "Roots on the negative real axis, no continuous time representation exist." when estimating a model from a signal. This means that one of the poles in the discrete time model, which is what is being estimated from data, landed on the legative real axis. No continuous-time system can ever create such a discrete-time model through sampling, and the some features of this package will work slightly worse if such a model is used, notably the EuclideanRootDistance and embedding. Optimal-transport based distances will not have a problem with this scenario.

To reduce the likelihood of this occurring, you may try to bandpass filter the signal before estimating the model, reduce the regularization factor if regularization was used, change the model order, or consider using the TLS fit method.

The difference between the pole locations in continuous and discrete time is vizualized in the pole diagrams below

pzmap(model, layout=2, sp=1, xlabel="Re", ylabel="Im", title="Continuous")
vline!([0], primary=false, l=(:black, :dash), sp=1)
pzmap!(tf(model,1),    sp=2, xlabel="Re", ylabel="Im", title="Discrete")

Computational performance improvements

Estimating a lot (1000s) of models might take a while. The bulk of the time is spent performing a matrix factorization, something that can be significantly sped up by

  • Using Flaot32 instead of Float64
  • Use MKL.jl instead of the default OpenBLAS (can yield about 2x performance improvement).

Further reading and examples

This notebook illustrates how the estimation of models in the presence of strong impulsive noise can be improved.

Type reference

Function reference

Docstrings

SpectralDistances.AR — Type
struct AR <: AbstractModel

Represents an all-pole transfer function, i.e., and AR model

Arguments:

  • a: denvec
  • ac: denvec cont. time
  • p: discrete time poles
  • pc: continuous time poles
  • b: Numerator scalar
source
SpectralDistances.ARMA — Type
struct ARMA{T} <: AbstractModel

Represents an ARMA model, i.e., transfer function

Arguments:

  • b: numvec
  • bc: numvec cont. time
  • a: denvec
  • ac: denvec cont. time
  • z: zeros
  • p: poles
source
SpectralDistances.IRLS — Type
IRLS <: FitMethod

Iteratively reqeighted least squares. This fitmethod is currently not recommended, it does not appear to produce nice spectra. (feel free to try it out though).

Arguments:

  • na::Int: number of roots (order of the system)
source
SpectralDistances.LS — Type
LS <: FitMethod

This fitmethod is a good default option.

Arguments:

  • na::Int: number of roots (order of the system). The number of peaks in the spectrum will be na÷2.
  • λ::Float64 = 0.01: reg factor

Computational performance improvements

Estimating a lot (1000s) of models might take a while. The bulk of the time is spent performing a matrix factorization, something that can be significantly sped up by

  • Using Flaot32 instead of Float64
  • Use MKL.jl instead of the default OpenBLAS (can yield about 2x performance improvement).
source
SpectralDistances.PLR — Type
PLR <: FitMethod

Pseudo linear regression. Estimates the noise components by performing an initial fit of higher order. Tihs fitmethod produces an ARMA model. Support for ARMA models is not as strong as for AR models.

Arguments:

  • nc::Int: order of numerator
  • na::Int: order of denomenator
  • initial::T = TLS(na=80): fitmethod for the initial fit. Can be, e.g., LS, TLS or any function that returns a coefficient vector
  • λ::Float64 = 0.0001: reg factor
source
SpectralDistances.TLS — Type
TLS <: FitMethod

Total least squares. This fit method is good if the spectrum has sharp peaks, in particular if the number of peaks is known in advance.

Arguments:

  • na::Int: number of roots (order of the system). The number of peaks in the spectrum will be na÷2.
  • λ::Float64 = 0: reg factor
source
SpectralDistances.fitmodel — Method
fitmodel(fm::LS, X::AbstractArray)

Computational performance improvements

Estimating a lot (1000s) of models might take a while. The bulk of the time is spent performing a matrix factorization, something that can be significantly sped up by

  • Using Flaot32 instead of Float64
  • Use MKL.jl instead of the default OpenBLAS (can yield about 2x performance improvement).
source
SpectralDistances.plr — Method
plr(y, na, nc, initial; λ=0.01)

Performs pseudo-linear regression to estimate an ARMA model.

Arguments:

  • y: signal
  • na: denomenator order
  • nc: numerator order
  • initial: fitmethod for the initial fit. Can be, e.g., LS, TLS or any function that returns a coefficient vector
  • λ: reg
source
SpectralDistances.residues — Method
residues(a::AbstractVector, b, r=roots(reverse(a)))

Returns a vector of residues for the system represented by denominator polynomial a

Ref: slide 21 https://stanford.edu/~boyd/ee102/rational.pdf Tihs methid is numerically sensitive. Note that if two poles are on top of each other, the residue is infinite.

source
SpectralDistances.spectralenergy — Method
spectralenergy(G::LTISystem)

Calculates the energy in the spectrum associated with G Ref: Robust and optimal control, Kemin Zhou, John C. Doyle, Keith Glover Lemma 4.6

source
SpectralDistances.change_precision — Method
change_precision(F, m::AbstractModel)

Changes the precision of all fields in m to F, e.g., F=Float64. This can be useful since the default precision for many operations in this package is Double64. This ensures that roots are calculated with high accuracy, but the high precision might not be required to evaluate distances etc.

source
SpectralDistances.hungariansort — Function
hungariansort(p1, p2)

takes two vectors of numbers and sorts and returns p2 such that it is in the order of the best Hungarian assignement between p1 and p2. Uses abs for comparisons, works on complex numbers.

source
SpectralDistances.reflect — Method
reflect(r::AbstractRoots)

Reflects unstable roots to a corresponding stable position (in unit circle for disc. in LHP for cont.)

source
ControlSystemsBase.tf — Method
ControlSystemsBase.tf(m::AR, ts)

Convert model to a transfer function compatible with ControlSystemsBase.jl

source
ControlSystemsBase.tf — Method
ControlSystemsBase.tf(m::AR)

Convert model to a transfer function compatible with ControlSystemsBase.jl

source