# Time-Frequency distances

For non-stationary signals, it is important to consider how the spectrum changes with time. This package has some, so far, basic support for time-frequency representations of non-stationary signals.

## Overview

There are two main approaches for time-frequency representations supported

- Model based
- Spectrogram based

### Model based

We define a custom fit method for fitting time varying spectra, `TimeWindow`

. It takes as arguments an inner fitmethod, the number of points that form a time window, and the number of points that overlap between two consecutive time windows:

`julia> fitmethod = TimeWindow(LS(na=2), 1000, 500)`

`TimeWindow{LS}(LS(2, 0.01), 1000, 500)`

`julia> y = sin.(0:0.1:1000);`

`julia> model = fitmethod(y)`

`Model estimation 11%|███▋ | ETA: 0:00:18 Model estimation 100%|███████████████████████████████████| Time: 0:00:02 TimeVaryingAR{AR{Vector{Float64}, DiscreteRoots{ComplexF64, Vector{ComplexF64}}, ContinuousRoots{ComplexF64, Vector{ComplexF64}}}, Vector{AR{Vector{Float64}, DiscreteRoots{ComplexF64, Vector{ComplexF64}}, ContinuousRoots{ComplexF64, Vector{ComplexF64}}}}}(AR{Vector{Float64}, DiscreteRoots{ComplexF64, Vector{ComplexF64}}, ContinuousRoots{ComplexF64, Vector{ComplexF64}}}[AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.997, 0.997), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.997, 0.997), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99702, 0.99702), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99702, 0.99702), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99703, 0.99703), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99703, 0.99703), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99703, 0.99703), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99702, 0.99702), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99701, 0.99701), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.997, 0.997), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99699, 0.99699), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99699, 0.99699), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.997, 0.997), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99701, 0.99701), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99702, 0.99702), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99703, 0.99703), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99703, 0.99703), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99703, 0.99703), AR{coeff: Float64, root: ComplexF64} Abs disc ∈: (0.99702, 0.99702)])`

This produces a custom model type, `TimeVaryingAR`

that internally stores a vector of `ContinuousRoots`

.

`TimeWindow`

currently only supports fitmethod `LS`

.

Accompanying this time-varying model is a time-aware distance, `TimeDistance`

. It contains an inner distance (currently only `OptimalTransportRootDistance`

supported), and some parameters that are specific to the time dimension, example:

`julia> dist = TimeDistance( inner = OptimalTransportRootDistance( domain = Continuous(), p = 2, weight = simplex_residueweight, ), tp = 2, c = 0.1, )`

`TimeDistance{OptimalTransportRootDistance{Continuous, typeof(identity), typeof(simplex_residueweight), Nothing}, Float64}(OptimalTransportRootDistance{Continuous, typeof(identity), typeof(simplex_residueweight), Nothing}(Continuous(), identity, SpectralDistances.simplex_residueweight, 0.01, 2, nothing), 2, 0.1)`

`tp`

is the same as `p`

but for the time dimension, and `c`

trades off the distance along the time axis with the distance along the frequency axis. A smaller `c`

makes it cheaper to transport mass across time. The frequency axis spans `[-π,π]`

and the time axis is the non-negative integers, which should give you an idea for how to make this trade-off.

Below, we construct a signal that changes frequency after half the time, and measure the distance between two different such signals. If the time penalty is large, it is cheaper to transport mass along the frequency axis, but if we make `c`

smaller, after a while the transport is cheaper along the time dimension

```
# Construct a signal that changes freq after half
signal(f1,f2) = [sin.((0:0.1:49.9).*f1);sin.((50:0.1:99.9).*f2)]
fm = TimeWindow(LS(na=2), 500, 0)
m = signal(1,2) |> fm # Signal 1
m2 = signal(2,1) |> fm # Signal 2 has the reverse frequencies
# First it is expensive to transport along time
dist = TimeDistance(
inner = OptimalTransportRootDistance(
domain = Continuous(),
p = 1,
weight = simplex_residueweight,
),
tp = 1,
c = 1.0, # c is large
)
evaluate(dist, m, m2)
```

`0.10010035006916945`

Then we make it cheaper

```
dist = TimeDistance(
inner = OptimalTransportRootDistance(
domain = Continuous(),
p = 1,
weight = simplex_residueweight,
),
tp = 1,
c = 0.01, # c is small
)
evaluate(dist, m, m2)
```

`0.010030949189558745`

#### Chirp example

Here we consider the estimation of the distance between two signals containing chirps, where the onset of the chirp differs. We start by creating some signals with different chirp onsets:

```
fs = 100000
T = 3
t = 0:1/fs:T
N = length(t)
f = range(1000, stop=10_000, length=N)
chirp0 = sin.(f.*t)
function chirp(onset)
y = 0.1sin.(20_000 .* t)
inds = max(round(Int,fs*onset), 1):N
y[inds] .+= chirp0[1:length(inds)]
y
end
using DSP
plot(spectrogram(chirp(0), window=hanning), layout=2)
plot!(spectrogram(chirp(1), window=hanning), sp=2)
```

We then define the fit method and the distance, similar to previous examples

```
fm = TimeWindow(LS(na=4, λ=1e-4), 20000, 0)
m = chirp(1) |> fm # This is the signal we'll measure the distance to
onsets = LinRange(0, 2, 21) # A vector of onset times
cv = exp10.(LinRange(-3, -0.5, 6)); # a vector of `c` values for the time-transport cost
```

```
6-element Vector{Float64}:
0.001
0.003162277660168376
0.010000000000000005
0.03162277660168378
0.10000000000000002
0.31622776601683794
```

We now calculate the distance to the base signal for varying onsets and varying time-transport costs.

```
dists = map(Iterators.product(cv, onsets)) do (c, onset)
m2 = chirp(onset) |> fm
dist = TimeDistance(
inner = OptimalTransportRootDistance(
domain = Continuous(),
p = 1,
weight = simplex_residueweight,
),
tp = 1,
c = c,
)
evaluate(
dist,
change_precision(Float64, m), # we reduce the precision for faster computations
change_precision(Float64, m2),
iters = 10000,
tol = 1e-2, # A coarse tolerance is okay for this example
)
end
plot(onsets, dists',
lab = cv',
line_z = log10.(cv)',
color = :inferno,
legend = false,
colorbar = true,
xlabel = "Onset [s]",
title = "Distance as function of onset and time cost"
)
```

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

The results are shown below. The figure indicates the cost `log10(c)`

using the color scale. We can see that the distance between the signals is smallest at `onset=1`

, which was the onset for the base signal. We also see that for small values of `c`

, it's cheap to transport along time. After increasing `c`

for a while it stops getting more expensive, indicating that it's now cheaper to transport along the frequency axis instead.

In this example, we chose the weight function `simplex_residueweight`

, which ensures that each time step has the same amount of spectral mass. Individual poles will still have different masses within each timestep, as determined by the pole's residue.

### Spectrogram based

The distance `ConvOptimalTransportDistance`

can operate on 2D measures on a grid, for example spectrograms. This distance supports calculating barycenters and barycentric coordinates. See Barycenters between spectrograms and Barycentric coordiantes.

This distance has a special option

`invariant_axis::Int = 0`

If this is set to 1 or 2, the distance will be approximately invariant to translations of the entire measure along the invariant axis. As an example, to be invariant to a spectrogram being shifted slightly in time, set `invariant_axis = 2`

. If this option is used, the distance is not differentiable and this option will be ignored when calculating barycenters etc.

## Dynamic Time Warping

This package interfaces with DynamicAxisWarping.jl and provides optimized methods for `DynamicAxisWarping.dtwnn`

. Below is an example of how to search for a query pattern `Qm`

in a much longer pattern `Ym`

```
searchresult = dtwnn(
Qm,
Ym,
OptimalTransportRootDistance(p = 1, weight = simplex_residueweight),
rad,
saveall = true,
tol = 1e-3,
)
```

Both `Qm`

and `Ym`

are expected to be of type `TimeVaryingAR`

.

`DynamicAxisWarping.dtwnn`

— Method`DynamicAxisWarping.dtwnn(q::TimeVaryingAR, y::TimeVaryingAR, dist, rad::Int, args...; kwargs...)`

Wrapper for `dtwnn`

. To save allocations between multiple calls to `dtwnn`

, you may manually create a helper object and call

```
using DynamicAxisWarping
h = DTWHelper(dist, q, y)
dtwnn(h, rad::Int, args...; kwargs...)
```

For examples of the combination of DTW and OT, see the following notebooks

## Distance profile

A distance profile between a query pattern `Qm`

and a much longer pattern `Ym`

can be computed efficiently with (example)

```
dist = TimeDistance(
inner = OptimalTransportRootDistance(p = 1, β = 0.5, weight = simplex_residueweight),
tp = 1,
c = 0.1,
)
res_tt = distance_profile(dist, Qm, Ym, tol=1e-3)
```

Both `Qm`

and `Ym`

are expected to be of type `TimeVaryingAR`

.

Missing docstring for `SlidingDistancesBase.distance_profile(od::TimeDistance, q::TimeVaryingAR, y::TimeVaryingAR)`

. Check Documenter's build log for details.

Missing docstring for `SlidingDistancesBase.distance_profile(od::ConvOptimalTransportDistance, q::DSP.Periodograms.TFR, y::DSP.Periodograms.TFR)`

. Check Documenter's build log for details.

## Docstrings

`SpectralDistances.TimeDistance`

— TypeThis is a time-aware distance. It contains an inner distance (currently only `OptimalTransportRootDistance`

and `KernelWassersteinRootDistance`

supported), and some parameters that are specific to the time dimension, example:

```
dist = TimeDistance(
inner = OptimalTransportRootDistance(
domain = Continuous(),
p = 2,
weight = simplex_residueweight,
),
tp = 2,
c = 0.1,
)
```

`tp`

is the same as `p`

but for the time dimension, and `c`

trades off the distance along the time axis with the distance along the frequency axis. A smaller `c`

makes it cheaper to transport mass across time. The frequency axis spans `[-π,π]`

and the time axis is the non-negative integers, which should give you an idea for how to make this trade-off.

`SpectralDistances.TimeVaryingAR`

— Type`TimeVaryingAR{T <: AbstractRoots} <: AbstractModel`

This model represents a rational spectrogram, i.e., a rational spectrum that changes with time. See `TimeWindow`

for the corresponding fit method and `TimeDistance`

for a time-aware distance.

Internally, this model stores a vector of `AR`

.

`SpectralDistances.TimeWindow`

— TypeWe define a custom fit method for fitting time varying spectra, `TimeWindow`

. It takes as arguments an inner fitmethod, the number of points that form a time window, and the number of points that overlap between two consecutive time windows:

**Example**

```
fitmethod = TimeWindow(TLS(na=2), 1000, 500)
y = sin.(0:0.1:100);
model = fitmethod(y)
```

`SlidingDistancesBase.distance_profile`

— Method`distance_profile(od::TimeDistance, q::TimeVaryingAR, y::TimeVaryingAR; normalize_each_timestep = false, kwargs...)`

Optimized method to compute the distance profile corresponding to sliding the short query `q`

over the longer `y`

.