Open ericphanson opened 3 years ago
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.
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 https://github.com/beacon-biosignals/TimeSpans.jl/issues/2 😁
That sounds useful
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).
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?
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.
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.
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 agetindex
function.Onda.jl does something smilar with its
Onda.Samples
which wrapschannel x time
signal data, and also allow indexing byTimeSpan
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 intoSpectrogram
s 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
?