Thanks for the contribution! Could you comment on how it differs or improves upon the mt_pgram multitaper spectrogram function already in DSP.jl?

Ah yes, it's a spectrogram. Nevermind!

Made a bit of progress here today: I reworked things to have an MTConfig object that stores the settings and temporary variables for mt_pgram, and added an in-place version that mt_pgram! that uses it. Then changed the multitaper spectrogram config to use it too. This way we don't have a separate implementation for mt_pgram vs multitaper_spectrogram! and that you can do both in-place, and even reuse config object s between them. We also have a multitaper coherence implementation that I'd like to add, which can also re-use the same config object.

I also changed the order to be freq x time to match spectrogram and added some docstrings. Still needs tests though.

I'll try to take a look this weekend, sorry for the delay.

No worries! It wasn't really ready yet before anyway. Now is a good time to take a look though, if you get a chance!

I'll try to address review comments and add tests on Monday (assuming other priorities don't come up first...). The cross-spectra stuff is probably best as a separate PR, so I'll move that out then too (right now it's just commented out and unfinished).

I'm going to force push to clean up the git history and properly credit some of the code's origins with co-author commits. Hope you don't mind!

Ok, I think this is ready for review. Thanks for the git help, @palday!

I don't have permission to formally ask for review on github, but I'd appreciate it if @kolia and @ararslan could give this a review too, since it's a big bunch of code. I will start making some PRs to Beacon private repos to try out this branch as well.

Choices I made that possibly should be questioned:

• configuration objects named MTConfig{T}, MTSpectrogramConfig{T}, MTCrossSpectraConfig{T}, MTCoherenceConfig{T} hold settings and workspaces. These can be re-used to perform the same operations re-using workspace variables on different signals of the same size (with the same settings). (If you're using multiple threads, you should use a separate workspace for each thread, but note FFTW can be threaded too, so don't oversubscribe!).
• I did not make a default numeric type, e.g. MTConfig(...) = MTConfig{Float64}(...), since it's not clear what the default should be. One option is to make e.g. MTConfig(signal) methods to extract size and eltype info from given signal data.
• allocate_output(config) produces an output object of the correct size given the configuration
• I did not name the "cross spectra" things as "csd" since there is not a normalization to make it a true "density". Is that annoying / too pendantic? I did not really know what to call those functions.
• I ended up putting the cross spectra and coherence stuff in this PR instead of splitting it out, since I realized some of my decisions were guided by needing that code to work too. However, if it's too bothersome to review this way please let me know and I'll split that out into a separate PR.
• mt_pgram now allocates a factor of ntaper more memory, because I save the results of each taper separately in tapered_spectra! before combining them in mt_pgram!. However, that memory is a workspace in the MTConfig object, so if you re-use the config you of course save a lot of memory. Why save each taper separately? Because then the same codepath and inner config objects can be reused for mt_cross_spectral!. I thought it was simpler and less error-prone this way, but if that's too big of an issue then I can make a separate code path for each case. This is fixed in 4e5c6a5
• I thought this one was kinda nice: we pass FFTW.MEASURE as the default flag when you create a config object, but FFTW.ESTIMATE when you use the out-of-place methods that create the config on the fly. The idea is if you're bothering to create a config object, then you're probably going to re-use it, so it's worth spending a bit of time optimizing the plan. If you're not going to re-use it (since you don't even get access to it with the out-of-place methods), then let's just get a quick plan as usual.
• AFAIK, the cross spectra and coherence code is the first multi-channel code in DSP.jl, and thus a convention of which dimension is the channel and which is time is needed. I chose the ordering channels x samples since that's what we use in Onda.jl, but of course the other ordering is equally valid.
• I tried to make as "simple" as possible by layering the config objects, so for example MTConfig settings are used basically whenever you're doing an FFT in this code, and then if you're doing a bunch of them as a spectrogram, MTSpectrogramConfig uses an MTConfig to govern the individual ones and just holds the configuration needed for how you are combining them.
• I tried to match the existing API, by naming things mt_* to match mt_pgram, and using fs, nfft, nw, etc, which are existing parameter names. I also tried to match which things are keyword arguments vs positional, especially between creating the config objects and the out-of-place functions. Reviewers, please take a look and see if I missed anything!

Quick benchmark

Using the first test from test/periodogram.jl as example data, we have

julia> using BenchmarkTools

julia> @btime mt_spectrogram!($out, $x0, $spec_config);
  14.186 μs (1 allocation: 2.12 KiB)

julia> @btime mt_spectrogram(x0, 256, 128; fs=10);
  638.610 μs (64 allocations: 118.73 KiB)

julia> @btime spectrogram(x0, 256, 128; fs=10);
  32.185 μs (47 allocations: 10.73 KiB)

You'll note that the in-place one still has one pesky allocation. Not sure where that is, but if someone spots it, let me know! And also the in-place multitaper spectrogram is faster than the existing out-of-place spectrogram, which is kind of fun (considering the multitaper one does like 7x more work), but almost certainly that advantage will go away as the data gets bigger.

Just seeing where this stands, were you going to add an option to weight each taper's estimation, and add a similar option to multitaper estimation beyond spectrograms?

Just seeing where this stands, were you going to add an option to weight each taper's estimation, and add a similar option to multitaper estimation beyond spectrograms?

Yep! Just pushed up a commit. Some updates:

• I am pretty sure the cross_spectral code is comupting power not energy (i.e. there's a normalization by the length of the sequence), and that's built into the tapering. I wrote up my understanding of this here: https://mne.discourse.group/t/understanding-mne-time-frequency-csd-multitaper-implementation/3043/2. Therefore I've adjusted the names to reflect that.
• The low_bias and eigenvalue weighting stuff seems a bit ad-hoc so I've pulled it out of the main API. At the moment, it's exposed by a dpss_config function to create a config with those parameters, which we need for the reference tests. I'm not totally sure how it should be exposed.
• The parameter r already is used for (inverse) normalization, so I've expanded it to allow per-taper normalization involving any desired taper_weights, which you can pass to the MTConfig constructor.
• I've added methods so we have all four of the methods discussed in add in-place mt_pgram!, and multitaper spectrogram, multitaper csd, multitaper coherence  #401 (comment)

I am guessing the codecov will take a hit with the new methods I added, so I'll add tests for them (...tomorrrow). I pushed up what I had so far so you could take a look already in case you had any comments.

I've added tests for the new methods and also added an error if the user does mt_cross_spectra (or mt_coherence, which uses that method) with complex input or a two-sided FFT, because I realized the current code assumes a one-sided FFT. I don't think it should be that hard to generalize it, but I'd prefer to leave that for later work in the interest of getting this in quicker.

I think I've addressed all your great feedback, please let me know if I've missed anything or as new things come up :).

Would be amazing to merge this :) Is there anything we can do to help?
@ararslan, @kolia, could you take another look at this, now that it's ready to get merged?

Can you add the new functions to the docs? Unless I missed it in the long and impressive file changes.

Can you add the new functions to the docs?

Good point, done!

Bump 🙂

Would really love a re-review @galenlynch! Or maybe another DSP maintainer can help out, e.g. @ararslan?

I've got a pile of branches depending on this branch and they've started to propagate further downstream...

Hey sorry that I disappeared. I'm happy with this now! I'll approve the PR and merge it tomorrow unless someone objects.

Whoo! Very exciting :). Could we tag a release?