[ENH] - Add multitaper method for spectral density estimation #317
Conversation
There was a problem hiding this comment.
Thanks for making this! I left some general feedback/questions. I compared this to the mne's psd_array_multitaper and the spectra here seem noiser. I wonder if mne windows the signal and then compute tapers on a window basis, averages per window, and averages across windows? I'll have to take a look into mne more. Here's a comparison:
TomDonoghue
changed the title
|
The tests weren't running for some reason, so I made an edit to add the new func to the API list, to quickstart the tests, and that seemed to work. I'm on board with adding a basic multitaper option to ndsp. Since I think @mwprestonjr has a bit more experience with this method, I'm curious on his thoughts on this, and I can also do a quick review if useful. @SM-Figueroa - can you also add a test that runs the new multitaper function? You can see examples in the other test files - it doesn't have to be a rigorous check for accuracy (since we probably don't know what the exact numerical outputs should be), but at least include a "smoke test" that runs the function, checks the inputs & outputs are all as expected in shape, type, etc (this at least catches if the function stops running / returning properly). |
|
I'll have some time to review this in the next few days. |
add default args for compute_spectrum_multitaper() add spectral.checks.check_mt_settings. this function sets values for bandwidth and num_windows. update DPSS args. NW/bandwidth corrected
|
I added a check_mt_settings function which sets default parameter values for compute_spectrum_multitaper (analogous to check_spg_settings). This function also uses the bandwidth parameter to compute NW, which should be used as the second argument to scipy.DPSS. With this change, the results should be comparable to mne.time_frequency.psd_array_multitaper. It's worth noting that later tapers in the Slepian sequence have lower spectral concentration. Babadi and Brown (IEEE 2014) suggest limiting the number of tapers to significantly less than n_samples * bandwidth * dt * 0.5 for this reason; MNE weights the tapered spectra based on the concentration ratios, and provides an option (low_bias=True/False) for dropping tapers with ratios less than 0.9: in addition to this weighting method, MATLAB's multitaper method includes an 'adaptive' weighting method that is frequency dependent. In short, there doesn't seem to be a single correct way to weight and recombine tapers. But we should probably have an option for weighting/dropping tapers. I can work on implementing this. |
|
I added a simple output test comparing the shape of the returned frequency array to the power array as done for other methods in test_power.py. I'm unsure of how to test the actual output values of the method for accuracy. Any thoughts would be much appreciated. |
|
Hey @SM-Figueroa & @mwprestonjr - sorry for dropping off here! In terms of status - it seems we have a basic working implementation, and a basic smoke test (checks the function runs), right - with no particular known issues? In terms of testing accuracy - that's actually quite hard to do to any detail. In To @mwprestonjr's point about weighting tapers - this seems like an important extension. Not knowing too much about this - do you think we can / should finish and merge this, and then try to extend it later, or is the weighting so important that it should really be added here? |
drop tapers with concentratio ratios less than 0.9, by default.
weight spectral estimates by concentratio ratio prior to combining. this argument can be set to false for simple averaging over estimates.
|
Hey @TomDonoghue! That's right. I think adding a basic test of accuracy as described would be good. I went ahead and implemented the weighting procedure. There's now a 'low_bias' argument that's True by default, which drops tapers with low concentration ratios. I also added an 'eigenvalue_weighting' option that's True by default, which weights each spectral estimate by the concentration ratio of its respective taper. The procedure and results are similar to MNE. |
|
This all looks great - thanks for the updates @mwprestonjr! Everything here looks good to me, so I'm going to go ahead and merge this in! Thanks for the work on this! |
This is an initial attempt at implementing a multi-taper method for estimating power spectral density. Refer to the issue raised here for more details.