Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add getindex methods to DSP.Spectrogram? #402

Open
ericphanson opened this issue Feb 11, 2021 · 7 comments
Open

Add getindex methods to DSP.Spectrogram? #402

ericphanson opened this issue Feb 11, 2021 · 7 comments

Comments

@ericphanson
Copy link
Contributor

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
Copy link
Contributor

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
Copy link

jrevels commented Feb 11, 2021

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
Copy link
Member

That sounds useful

@ericphanson
Copy link
Contributor Author

ericphanson commented Apr 27, 2021

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
Copy link
Member

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
Copy link
Member

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
Copy link
Contributor Author

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

5 participants