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

1. Model based
2. 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.

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(q::TimeVaryingAR, y::TimeVaryingAR, dist, rad::Int, args...; kwargs...)

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

• DTW-OT: Introduction
• DTW-OT: Detection

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.

Docstrings

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

TimeVaryingAR{T <: AbstractRoots} <: AbstractModel

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

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