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 ControlSystems.LTISystem, so many of the functions from the ControlSystems.jl toolbox work on these models as well. When acting like a ControlSystems.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)
┌ Warning: Roots on negative real axis, no corresponding continuous time representation exists. Consider prefiltering the signal or decreasing the regularization factor.
└ @ SpectralDistances ~/work/SpectralDistances.jl/SpectralDistances.jl/src/ltimodels.jl:83
AR{coeff type: Float64, root type: Complex{Float64}}(
b: 24.1879920231178
Cont. poles: 6-element ContinuousRoots{Complex{Float64},Array{Complex{Float64},1}}:
-0.543568696277846 - 2.235125945694924im
-0.6109325622076213 - 0.9074245948536653im
-1.0012533692830874 + 0.0im
-0.6109325622076213 + 0.9074245948536653im
-0.543568696277846 + 2.235125945694924im
-0.9922236868551086 + 3.141592653589793im
Abs disc: [0.5806723072668087, 0.5428443967836869, 0.3674186412165771, 0.5428443967836869, 0.5806723072668087, 0.37075133882107425]
Cont. coeffs: 7-element Array{Float64,1}:
1.0
4.302479573109131
13.412659026707326
25.641561424234055
29.578614490165485
20.3377406098773
6.290458646503813
Spectral energy: 1.0442449114635568

julia> change_precision(Float32, model) # Can be useful to reduce the computational cost of some distances
AR{coeff type: Float32, root type: Complex{Float32}}(
b: 24.1879920231178
Cont. poles: 6-element ContinuousRoots{Complex{Float32},Array{Complex{Float32},1}}:
-0.5435687f0 - 2.235126f0im
-0.6109326f0 - 0.90742457f0im
-1.0012534f0 + 0.0f0im
-0.6109326f0 + 0.90742457f0im
-0.5435687f0 + 2.235126f0im
-0.9922237f0 + 3.1415927f0im
Abs disc: Float32[0.5806723, 0.5428444, 0.36741865, 0.5428444, 0.5806723, 0.37075135]
Cont. coeffs: 7-element Array{Float32,1}:
1.0
4.3024797
13.412659
25.641562
29.578615
20.33774
6.2904587
Spectral energy: 1.0442449

julia> pzmap(model)
┌ Warning: markershape x is unsupported with Plots.PlotlyBackend().  Choose from: [:none, :auto, :circle, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross, :pentagon, :hexagon, :octagon, :vline, :hline]
└ @ Plots ~/.julia/packages/Plots/uCh2y/src/args.jl:1195
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")
┌ Warning: markershape x is unsupported with Plots.PlotlyBackend().  Choose from: [:none, :auto, :circle, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross, :pentagon, :hexagon, :octagon, :vline, :hline]
└ @ Plots ~/.julia/packages/Plots/uCh2y/src/args.jl:1195
┌ Warning: markershape x is unsupported with Plots.PlotlyBackend().  Choose from: [:none, :auto, :circle, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross, :pentagon, :hexagon, :octagon, :vline, :hline]
└ @ Plots ~/.julia/packages/Plots/uCh2y/src/args.jl:1195

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).

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

Docstrings

SpectralDistances.ARType
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.ARMAType
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.IRLSType
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.LSType
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.PLRType
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.TLSType
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.fitmodelMethod
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.plrMethod
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.residuesMethod
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.change_precisionMethod
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.hungariansortFunction
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