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

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 distances`

`AR{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

`SpectralDistances.AR`

`SpectralDistances.ARMA`

`SpectralDistances.AbstractModel`

`SpectralDistances.ContinuousRoots`

`SpectralDistances.ContinuousRoots`

`SpectralDistances.DiscreteRoots`

`SpectralDistances.DiscreteRoots`

`SpectralDistances.FitMethod`

`SpectralDistances.HungarianAssignement`

`SpectralDistances.IRLS`

`SpectralDistances.LS`

`SpectralDistances.PLR`

`SpectralDistances.SortAssignement`

`SpectralDistances.TLS`

## Function reference

`ControlSystemsBase.denvec`

`ControlSystemsBase.tf`

`ControlSystemsBase.tf`

`PolynomialRoots.roots`

`SpectralDistances.change_precision`

`SpectralDistances.checkroots`

`SpectralDistances.coefficients`

`SpectralDistances.domain_transform`

`SpectralDistances.examplemodels`

`SpectralDistances.fitmodel`

`SpectralDistances.fitmodel`

`SpectralDistances.fitmodel`

`SpectralDistances.fitmodel`

`SpectralDistances.hungariansort`

`SpectralDistances.ls`

`SpectralDistances.normalization_factor`

`SpectralDistances.normalize_energy`

`SpectralDistances.plr`

`SpectralDistances.polar`

`SpectralDistances.reflect`

`SpectralDistances.residues`

`SpectralDistances.residues`

`SpectralDistances.residueweight`

`SpectralDistances.roots2poly`

`SpectralDistances.simplex_residueweight`

`SpectralDistances.spectralenergy`

`SpectralDistances.spectralenergy`

`SpectralDistances.unitweight`

## 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

`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

`SpectralDistances.AbstractModel`

— Typeabstract type AbstractModel <: ControlSystemsBase.LTISystem end

`SpectralDistances.FitMethod`

— TypeAbstract type that represents a way to fit a model to data

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

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

`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:**

`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

`PolynomialRoots.roots`

— Function`roots(m::AbstractModel)`

Returns the roots of a model

`SpectralDistances.checkroots`

— Method`checkroots(r::DiscreteRoots)`

prints a warning if there are roots on the negative real axis.

`SpectralDistances.coefficients`

— Method`coefficients(::Domain, m::AbstractModel)`

Return all fitted coefficients

`SpectralDistances.examplemodels`

— Function`examplemodels(n = 10)`

Return `n`

random models with 6 complex poles each.

`SpectralDistances.fitmodel`

— Method`fitmodel(fm::IRLS, X::AbstractArray)`

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

`SpectralDistances.fitmodel`

— Method`fitmodel(fm::PLR, X::AbstractArray)`

`SpectralDistances.fitmodel`

— Method`fitmodel(fm::TLS, X::AbstractArray)`

`SpectralDistances.ls`

— Function`ls(A, y, λ=0)`

Regularized Least-squares

`SpectralDistances.normalization_factor`

— Method`normalize_energy(r)`

Returns the factor that, when used to multiply the poles, results in a system with unit spectral energy.

`SpectralDistances.normalize_energy`

— Method`normalize_energy(r)`

Returns poles scaled to achieve unit spectral energy.

`SpectralDistances.plr`

— Method`plr(y, na, nc, initial; λ=0.01)`

Performs pseudo-linear regression to estimate an ARMA model.

**Arguments:**

`SpectralDistances.residues`

— Method`residues(r::AbstractRoots)`

Returns a vector of residues for the system represented by roots `r`

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

`SpectralDistances.roots2poly`

— Method`roots2poly(roots)`

Accepts a vector of complex roots and returns the polynomial with those roots

`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

`SpectralDistances.spectralenergy`

— Method`spectralenergy(d::TimeEvolution, a::AbstractVector, b)`

Calculates the energy in the spectrum associated with transfer function `b/a`

`SpectralDistances.ContinuousRoots`

— Type`ContinuousRoots{T, V <: AbstractVector{T}} <: AbstractRoots{T}`

Represents roots in continuous time

`SpectralDistances.ContinuousRoots`

— Method`ContinuousRoots(r::DiscreteRoots) = begin`

Represents roots of a polynomial in continuous time

`SpectralDistances.DiscreteRoots`

— Type`DiscreteRoots{T, V <: AbstractVector{T}} <: AbstractRoots{T}`

Represent roots in discrete time

`SpectralDistances.DiscreteRoots`

— Method`DiscreteRoots(r::ContinuousRoots) = begin`

Represents roots of a polynomial in discrete time

`SpectralDistances.HungarianAssignement`

— Type`HungarianAssignement <: AbstractAssignmentMethod`

Sort roots using Hungarian method

`SpectralDistances.SortAssignement`

— Type`SortAssignement{F} <: AbstractAssignmentMethod`

Contains a single Function field that determines what to sort roots by.

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

`SpectralDistances.domain_transform`

— Method`domain_transform(d::Domain, e::AbstractRoots)`

Change the domain of the roots

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

`SpectralDistances.polar`

— Method`polar(e::Number)`

magnitude and angle of a complex number

`SpectralDistances.reflect`

— Method`reflect(r::AbstractRoots)`

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

`SpectralDistances.residueweight`

— Method`residueweight(e::AbstractRoots)`

Returns a vector where each entry is roughly corresponding to the amount of energy contributed to the spectrum be each pole. See also `simplex_residueweight`

for a normalized version.

`SpectralDistances.simplex_residueweight`

— Method`simplex_residueweight(x)`

Returns a vector where each entry is roughly corresponding to the amount of energy contributed to the spectrum be each pole, normalized to sum to 1. See `residueweight`

for a non-normalized version.

`SpectralDistances.unitweight`

— Method`unitweight(e)`

A weighting function that returns a vector of uniform weights that sum to 1.

`ControlSystemsBase.tf`

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

Convert model to a transfer function compatible with ControlSystemsBase.jl

`ControlSystemsBase.tf`

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

Convert model to a transfer function compatible with ControlSystemsBase.jl

`ControlSystemsBase.denvec`

— Method`ControlSystemsBase.denvec(::TimeEvolution, m::AbstractModel)`

Get the denominator polynomial vector