Add getindex methods to DSP.Spectrogram?

[ericphanson]

I think it would be convienent to be able to index into a Spectrogram with times or frequencies, and obtain the power in the corresponding bin. Spectrogram already has all the necessary information (the times, frequencies, and powers), so it just needs a getindex function.

Onda.jl does something smilar with its Onda.Samples which wraps channel x time signal data, and also allow indexing by TimeSpan objects to return a view of the corresponding samples. I think it would be nice to use something like Invenia's Intervals.jl to allow indexing into Spectrograms with intervals of times or frequencies and obtain a corresponding view of the power from the spectrogram.

Does this seem like a good idea for DSP.Spectrogram?

[kleinschmidt]

Also relevant is that https://github.com/beacon-biosignals/TimeSpans.jl was spun out of Onda.jl so it's a pretty lightweight dependency which provides the convenience functions that Onda uses to work with time spans.

[jrevels]

Also relevant is that https://github.com/beacon-biosignals/TimeSpans.jl was spun out of Onda.jl so it's a pretty lightweight dependency which provides the convenience functions that Onda uses to work with time spans.

should also ref beacon-biosignals/TimeSpans.jl#2 😁

[galenlynch]

That sounds useful

[ericphanson]

I've started using AxisKeys.jl for this in the hopefully soon-to-be open sourced OndaDSP, which connects DSP.jl to Onda.jl. For example, for a multichannel periodogram of an Onda samples object (which represents a regularly sampled multichannel signal), there's the method

@views function DSP.mt_pgram!(output::AbstractMatrix, samples::Samples, config::MTConfig{T}) where {T}
    if samples.info.sample_rate != config.fs
        throw(ArgumentError("`samples` does not have the same `sample_rate` as `config.fs`; got `samples.info.sample_rate`=$(samples.info.sample_rate) and `config.fs` = $(config.fs)."))
    end
    require_decoded(samples)
    for c in 1:channel_count(samples)
        mt_pgram!(output[c, :], samples.data[c, :], config)
    end
    unit = sample_unit(samples)
    return KeyedArray(with_unit(output, unit^2); channel=samples.info.channels, freq = config.freq * Hz)
end

(and out-of-place variants etc). Then for example, with spectrograms, you can do

julia> waves = sine_waves([60Hz, 45Hz, 30Hz])
Samples (00:10:00.000000000):
  info.kind: "synthetic"
  info.channels: ["c1", "c2", "c3"]
  info.sample_unit: "microvolt"
  info.sample_resolution_in_unit: 0.25
  info.sample_offset_in_unit: 0
  info.sample_type: Int16
  info.sample_rate: 200 Hz
  encoded: false
  data:
3×120000 Matrix{Float64}:
 0.0  0.951052  -0.587811  -0.587747   0.951076  -7.85405e-5  -0.951027   0.587874  …   7.85405e-5  -0.951076   0.587747   0.587811  -0.951052   2.32385e-12
 0.0  0.98769    0.308995  -0.891023  -0.587747   0.707148     0.808975  -0.454064     -0.707148     0.587747   0.891023  -0.308995  -0.98769   -5.53307e-12
 0.0  0.809022   0.951052   0.308995  -0.587811  -1.0         -0.587747   0.309069      1.0          0.587811  -0.308995  -0.951052  -0.809022   1.16193e-12

julia> spec = mt_spectrogram(waves)
3-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   channel ∈ 3-element Vector{String}
→   freq ∈ 257-element Vector{Quantity{Float64, 𝐓^-1, Unitful.FreeUnits{(Hz,), 𝐓^-1, nothing}}}
□   time ∈ 598-element StepRange{Nanosecond,...}
And data, 3×257×598 Array{Quantity{Float64, 𝐋^4 𝐌^2 𝐈^-2 𝐓^-6, Unitful.FreeUnits{(μV^2,), 𝐋^4 𝐌^2 𝐈^-2 𝐓^-6, nothing}}, 3}:
[showing 3 of 598 slices]
[:, :, 1] ~ (:, :, Nanosecond(1000000000)):
                 (0.0 Hz)    (0.390625 Hz)     (0.78125 Hz)     (1.17188 Hz)  …     (98.8281 Hz)     (99.2188 Hz)     (99.6094 Hz)       (100.0 Hz)
  ("c1")  1.36997e-7 μV^2  9.51958e-7 μV^2  1.89614e-6 μV^2  1.58865e-6 μV^2     5.78691e-6 μV^2  6.89311e-6 μV^2  3.45635e-6 μV^2  4.96658e-7 μV^2
  ("c2")  3.56939e-7 μV^2  2.48233e-6 μV^2  4.94787e-6 μV^2  4.15016e-6 μV^2     2.22651e-6 μV^2  2.65622e-6 μV^2  1.33314e-6 μV^2  1.91763e-7 μV^2
  ("c3")  1.00183e-6 μV^2  6.98199e-6 μV^2  1.39443e-5 μV^2  1.17342e-5 μV^2     7.92392e-7 μV^2  9.45911e-7 μV^2   4.7493e-7 μV^2  6.83439e-8 μV^2

[:, :, 300] ~ (:, :, Nanosecond(298505000000)):
                 (0.0 Hz)    (0.390625 Hz)     (0.78125 Hz)     (1.17188 Hz)  …     (98.8281 Hz)     (99.2188 Hz)     (99.6094 Hz)       (100.0 Hz)
  ("c1")  1.17957e-6 μV^2  1.82058e-6 μV^2  1.04293e-6 μV^2  1.23734e-6 μV^2     2.57473e-6 μV^2  2.73426e-6 μV^2  2.56513e-6 μV^2  1.20436e-6 μV^2
  ("c2")  9.02885e-7 μV^2  2.35708e-6 μV^2  3.09894e-6 μV^2  2.80178e-6 μV^2     1.71676e-6 μV^2   1.7806e-6 μV^2  1.71387e-6 μV^2  8.25888e-7 μV^2
  ("c3")  1.08946e-6 μV^2  5.93824e-6 μV^2  1.11714e-5 μV^2  9.45796e-6 μV^2     8.51161e-7 μV^2  8.74268e-7 μV^2  8.50162e-7 μV^2  4.13851e-7 μV^2

[:, :, 598] ~ (:, :, Nanosecond(595015000000)):
                 (0.0 Hz)    (0.390625 Hz)     (0.78125 Hz)     (1.17188 Hz)  …     (98.8281 Hz)     (99.2188 Hz)     (99.6094 Hz)       (100.0 Hz)
  ("c1")   1.5836e-6 μV^2  2.15721e-6 μV^2  7.12283e-7 μV^2  1.10119e-6 μV^2     1.32989e-6 μV^2  1.12256e-6 μV^2  2.21976e-6 μV^2  1.47862e-6 μV^2
  ("c2")  1.61773e-6 μV^2  2.19307e-6 μV^2  6.78029e-7 μV^2  1.03627e-6 μV^2     1.04932e-6 μV^2  6.34086e-7 μV^2  2.21238e-6 μV^2  1.65619e-6 μV^2
  ("c3")  1.27885e-6 μV^2  3.68262e-6 μV^2  5.17916e-6 μV^2  4.53877e-6 μV^2     9.78165e-7 μV^2  7.19442e-7 μV^2  1.66107e-6 μV^2  1.16052e-6 μV^2

julia> dropdims(mean(spec(channel=in(("c1", "c2")), freq=Interval(20Hz, 30Hz)); dims=(:channel,:freq)); dims=(:channel,:freq))
1-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   time ∈ 599-element StepRange{Dates.Nanosecond,...}
And data, 599-element Vector{Unitful.Quantity{Float64, 𝐋^4 𝐌^2 𝐈^-2 𝐓^-6, Unitful.FreeUnits{(μV^2,), 𝐋^4 𝐌^2 𝐈^-2 𝐓^-6, nothing}}}:
  Dates.Nanosecond(1000000000)     4.579166986065568e-6 μV^2
  Dates.Nanosecond(2000000000)     4.594530756686665e-6 μV^2
  Dates.Nanosecond(3000000000)     4.609899610640189e-6 μV^2
  Dates.Nanosecond(4000000000)     4.625273098048024e-6 μV^2
  ⋮                                
  Dates.Nanosecond(597000000000)   4.609899610671954e-6 μV^2
  Dates.Nanosecond(598000000000)   4.594530756819026e-6 μV^2
  Dates.Nanosecond(599000000000)   4.579166985782873e-6 μV^2

That's a lot of text, but the cool thing is the output of mt_spectrogram is an object you can treat like an array, but also knows the names of its channels, the time dimension, freq dimension etc. So you can easily answer questions like: what's the average power over time in the 20-30Hz frequency band over channels "c1" and "c2"?

I'm not sure if DSP wants to depend on AxisKeys but I think it's pretty nice! (I'm also using Unitful for units here but that's kinda a separate thing).

[rob-luke]

I just looked up AxisKeys and it is super cool. I will definitely start using it in my local projects.

However, for DSP.jl it would create additional dependencies and thus increased maintenance for this project (particularly if the dep is by a single author or small org). Could this be implemented without adding more dependencies?

[galenlynch]

Yeah I'm also enthusiastic about using AxisKeys for my own projects, but don't really want to add a dependency unless it can be avoided.

[ericphanson]

don't really want to add a dependency

Sounds good, I think it makes sense to keep DSP minimal and light.

Could this be implemented without adding more dependencies?

Yeah, maybe a minimal getindex method would be helpful already-- that's what I was thinking when I filed this issue. But now I also worry it could easily scope-creep (e.g. I can index into a DSP.Spectrogram, so why can't I broadcast over it?) and am starting to think maybe it's better to leave that kind of functionality to downstream packages where it probably can get more development attention than a partial reimplementation here.