[ENH] Add tutorial on time-frequency source estimation with STC viewer GUI

[alexrockhill]

Related to https://github.com/mne-tools/mne-python/issues/10721, GSoC.

So I realize there aren't really any tutorials on time-frequency source space resolution (there is DICS but that's just frequency in the tutorial) so I think this fills a very important gap in MNE functionality. Also, since it is the most complete use case (you can always average across time or trials) it covers a lot of other use cases so I'm pretty excited about adding this.

Basically, time-frequency resolved MEG data is filtered using an LCMV beamformer. This yields a list of lists of stcs where the outer layer is frequencies and the inner layer is trials. The above referenced issue is the discussion for why this is hard to keep track of and could use the introduction of a new object. There's a lot that could be done but the priorities in order are:

• [x] Make the user interface (in pyvistaqt first)
  • [ ] Once https://github.com/mne-tools/mne-python/pull/10913 is merged, use the abstraction
• [x] Make a TFRVolSourceEstimate class to hold the list of list stcs and be returned by apply_lcmv_epochs_tfr
• [x] Allow other functions to return corresponding TFR source objects as much as possible? This is a bit scope creep but it would be inconsistent to only have this functionality for LCMV. Making a new class for each existing source estimate class seems like a ton of overhead but potentially the more flexible source vertices+time+frequency+trials stc classes could supersede the source vertices+time source estimates and those could just be the flat case of the more general source estimates, that might be a huge refactoring though since it's so central to MNE.

Now is a great time to take five minutes to review the future direction if you have a minute. I will be out of town for the next few days so there won't be much movement.

cc @britta-wstnr @larsoner @drammock

[britta-wstnr]

I think there is some mixing of two concepts here, let me try and make this clearer:

Very generally speaking:

1. LCMV beamformers are for time-domain data (i.e., raw data, epochs, or evoked data).
2. DICS beamformers are for frequency-domain data (i.e., spectral power - in the supported case non-time-resolved power).

Some exceptions that come to mind (@alexrockhill and I have discussed some of these, so I am adding them for clarity):

1. The Hilbert beamformer or any type of bandpass-filtered signal. This is time-domain data, but (implicitly) also frequency-linked due to the filtering. But there is no power estimate (amplitude instead).
2. Our DICS code has e.g. apply_dics_evoked, which applies a DICS (i.e. frequency-resolved) beamformer to time-domain data. The use cases are not very typical I think - I would guess that people would maybe use a bandpass-filtered LCMV for this.

Now back to the proposal. I would suggest to have a DICS beamformer with this functionality. Meaning: We use make_dics to make a DICS beamformer (using a cross-spectral density matrix based on the whole epoch) and then apply this step-wise to a time-resolved power estimate (i.e., time-frequency representation, TFR). To me, it is not so ideal to use an LCMV beamformer for that. The output would be one stc per time point we estimated power in the TFR at.

Alternatively, we can go with the Hilbert approach, where we would use bandpass-filtered time-domain data. Here the output would be one stc per frequency-band we filtered the epochs at.

Hope this will add some clarity with respect to frequency- and time-domain beamforming.

[alexrockhill]

Ok, I'm totally on the same page with one thing that I don't quite understand. I will make it more of a priority to add the DICS support sooner rather than later. The one thing I don't understand:

To me, it is not so ideal to use an LCMV beamformer for that.

Why is it appropriate to use the Hilbert transform and not any other time-frequency decomposition? You can get amplitude and not power from Morlet wavelets or you can square the Hilbert amplitude and get power right? I'm not sure what's different about the Hilbert transform compared to other time-frequency methods. From our earlier conversation, I thought you could apply a LCMV or DICS filter with the tradeoff being that the DICS is discontinuous in time and the LCMV is discontinuous in frequencies so from that I would think it would be nice to be able to use either depending on your analysis.

[britta-wstnr]

Let me be a bit more specific, hope this will make it clearer. It is not wrong to apply a DICS filter to time-domain data or an LCMV filter to frequency-domain data. But I think it is not ideal and I am -1 on actually supporting this.

My main reason: If you take a time-domain signal and convert it from broadband to frequency-bandpassed, you use a (frequency) filter that has different characteristics (e.g. passband) than if you take this signal and transform it into the frequency-domain (e.g. by using a Fourier transform). Ultimately, you will end up with slightly (or even really) different frequencies being resolved in the signal. So why would you want to use e.g. an LCMV beamformer created on a bandpass-filtered signal resolving 14-17 Hz (toy example), to a power estimate that is 15 +/- 1 Hz, if you can instead apply a DICS beamformer that matches this frequency resolution perfectly? I am not convinced this is a good idea.

Thus, for the Hilbert-case (which operates on a bandpass-filtered time-domain signal), we use an LCMV (which uses a covariance matrix based on the same bandpass-filtered signal), while for a TFR reconstruction, we should rather use a DICS beamformer that is based on a cross-spectral density matrix that uses power estimates of the same frequency resolution.

Furthermore, this way we do not have to worry about some of the not-shared features between LCMV and DICS (e.g. complex valued filters). Plus, not supporting this might lead to less confusion in users :)

From our earlier conversation, I thought you could apply a LCMV or DICS filter with the tradeoff being that the DICS is discontinuous in time and the LCMV is discontinuous in frequencies so from that I would think it would be nice to be able to use either depending on your analysis.

Remember that the filter as such is never time-resolved (not to be mixed up with time-domain!), the time-resolution is only recovered in the application of the filter to the data (where, theoretically, it is possible to apply both LCMV and DICS filters to a time-resolved signal).

[alexrockhill]

So why would you want to use e.g. an LCMV beamformer created on a bandpass-filtered signal resolving 14-17 Hz (toy example), to a power estimate that is 15 +/- 1 Hz, if you can instead apply a DICS beamformer that matches this frequency resolution perfectly?

Because the DICS cannot be computed on a single time point and instead requires a window, so if you needed time resolution at your sampling rate but didn't mind if your frequencies were a bit smeared, I think it would be reasonable to use LCMV.

Thus, for the Hilbert-case (which operates on a bandpass-filtered time-domain signal), we use an LCMV (which uses a covariance matrix based on the same bandpass-filtered signal), while for a TFR reconstruction, we should rather use a DICS beamformer that is based on a cross-spectral density matrix that uses power estimates of the same frequency resolution.

I don't follow the rest of the argument about why a collection of bandpassed Hilbert transformed time-frequency data is any different than any other time-frequency decomposition. I understand everything you said and I think that's all reasonable but I don't understand how it applies to this point. You could take the covariance and not cross-spectral density of the TFR epochs and I'm not convinced that's any different than doing that for the Hilbert method.

Thanks for writing back, sorry for a few more questions :)

[britta-wstnr]

Because the DICS cannot be computed on a single time point and instead requires a window, so if you needed time resolution at your sampling rate but didn't mind if your frequencies were a bit smeared, I think it would be reasonable to use LCMV.

The DICS and the LCMV spatial filter are not different from each other in terms of how much data they need to go into the CSD or covariance matrix (i.e. the time interval the filter is built on). The signal you apply the spatial filter to is of course bound in its time resolution. However, if you want the high time resolution of a Hilbert signal, why not use an LCMV beamformer for that? Why are you so focused on mixing Hilbert input with a DICS spatial filter? If you instead use an LCMV, your frequencies in the spatial filter and the input signal you apply it to will match up perfectly.

Basically, it boils down to: why would you mix-and-match if staying within the time-domain or the frequency-domain will give you a methodologically more sound and cleaner result?

[alexrockhill]

Basically, it boils down to: why would you mix-and-match if staying within the time-domain or the frequency-domain will give you a methodologically more sound and cleaner result?

I'm still not following, you have to mix and match because you want time-frequency resolved data so you can either 1) convert to time-frequency with many bandpass filtered signals and apply Hilbert or use Morlet wavelets and then compute the covariance matrix and apply the LCMV beamformer or 2) compute the DICS beamformer on sliding time-windows using cross-spectral density, right?

To me the bandpass filter and sliding time windows are unideal because they inexactly resolve in frequency and time respectively. So wouldn't you want to select the more reasonable one for your use-case or compare results?

EDIT: After re-reading, it sounds like you associate the Hilbert transform with the time domain in some way that is different than Morlet wavelets. I would question this assumption: they're different methods but fundamentally they are trying to compute the same thing; the time-frequency amplitude or power. I've read and understand what you're saying but I don't agree that the Hilbert transform is somehow more closely associated with the time domain and therefore more cleanly associated with using an LCMV beamformer. Fundamentally, I think it's an empirical question of which works better based on that method's tradeoffs with estimating time-frequencies (i.e. choosing things like the bandwidth of the filter for Hilbert and the number of cycles for Morlet wavelets). Undoubtably, the practical experience of someone who uses these methods is very valuable and should be used to inform which method gets used or highlighted in the MNE tutorial, but I think it's important not to conflate what is mathematically inconsistent with what empirically just doesn't work as well. As I mentioned in the meeting, I think something that would really improve the methods in M/EEG in general, not just in MNE, is a really good set of decoding datasets to test methodological choices on. Maybe some of the mne-bids pipeline data that @agramfort is working on would be good for that but it's a bit outside of GSoC. Yes, if some choice of method works well for one dataset that's not super helpful in all cases, but if some methodological choice is consistently better across several datasets, I think that's much better evidence than intuition based on prior experience, again, although that is super helpful too. I don't want to dwell too much on beamforming choices but this is super relevant to the API; should the LCMV beamformer apply method take in an EpochsTFR object or Epochs and a set of frequencies and bandwidths to apply the Hilbert? I would argue that the first option is both more elegant/fits within existing MNE structures better and, until I see data or a convincing argument about why it's incoherent to use Morlet wavelets/Welch/Stockwell/whatever with LCMV, it is agnostic to the way the EpochsTFR was constructed and so is maximally abstract and generalizable. This is also on theme with trying to make a TFRVolSourceEstimate that can be returned both by the LCMV apply method, the DICS method and the apply_inverse method; although I don't want to put a lot of effort into developing all these methods, the API won't be very good if it isn't ready for these methods to be introduced in the future. That way we can highlight what is empirically shown to be the best method in the tutorial but still allow for other methods in a way that is really good for the code base because it leaves the time-frequency stuff to be done in a different, separate module of the code.

[alexrockhill]

Ok @britta-wstnr and @agramfort, if you want to try running through this tutorial, I think now would give a good indication of the direction I think we should go (and that we've agreed upon).

The main reason I bring this up is that the LCMV (of Hilbert data) and DICS apply functions returns a list of list of stcs but the LCMV is nested frequencies first and then epochs and then the stcs are shape (n_verts, n_samples) whereas the DICS is nested sliding time windows first then epochs and then the stcs are shape (n_verts, n_freqs). I think it's pretty difficult to rearrange the dimensions in lists of list so it really seems to be crying out for a class abstraction that will arrange them correctly for you. That's the next step that I'll try and get done today which is the main thing holding up making an actual GUI.

As for whether it's okay to use Morlet wavelets instead of Hilbert and why you can't use power as input to a minimum norm estimation, I think we should discuss at the next meeting and I'm looking forward to learning about how these work :)

It's not crucial that everything be resolved with the source estimation before moving forward (it really needs to be optimized for parallelization also), I just wanted to make sure the class worked so that it could be used as input to the GUI without having to be changed.

I have no idea what the data looks like and I haven't invented any new methods, just used the old ones, so I assume that it will be reasonable but a time-frequency viewer is sorely needed with this high of a dimension of data; it's just very hard to look at. Will have a reasonable screenshot soon!

EDIT: Also very likely the final minimum norm source estimation method is wrong since it's done with wavelets only in mne.minimum_norm.source_induced_power but it would be good to know why that is so. I'll switch it over to that code actually now that I think about it.

EDIT2: Actually I tried source_induced_power and waited five minutes but it didn't finish computing so I'll wait to discuss the best options.

[alexrockhill]

Notes per method:

1a) The LCMV on the Hilbert transformed data works great but is quite slow because the bandpass filter has to be computed on the raw data before epoching. 1b) The LCMV on Morlet wavelets seems to be efficient and has the fewest drawbacks computationally but we'll see about the results 2) The DICS is quite slow as well and there are issues with the sliding window size for lower frequencies since it must be larger. I had to go up to 90 Hz in the tests with the testing sampling frequency of 300 Hz because I kept getting warnings that the filter was longer than the data. This seems like a big drawback for this method in making a single TFR source estimate with the same sliding windows. If you allowed for bigger and smaller windows the code would be very complicated as it would no longer be matrix based... this is a tough one. 3) I imagine just applying the minimum norm to time-frequency epochs data may be suboptimal and incorrect but it computes in a reasonable time everything works well with window sizes since the time-frequencies are computed beforehand.

[agramfort]

@alexrockhill it would help to go over a rendered doc page / tutorial to see the user perspective on this. Is there any? Sorry if I let you babysit me here ...

[alexrockhill]

@alexrockhill it would help to go over a rendered doc page / tutorial to see the user perspective on this. Is there any? Sorry if I let you babysit me here ...

There will be soon since #10940 merged, it was too much memory, causing a segfault in circle but it runs locally. Let me see about now, otherwise I'll comment out 2/3 methods (LCMV, DICS, MNE).

[britta-wstnr]

Okay, I can only reiterate my points from before:

1. Hilbert outputs the same time-resolution as you put in, thus I sort this into a time-domain signal. The Morlet wavelet approach mostly does not do this - you traditionally have not the same time steps as in your input signal (but fewer).
2. And more importantly: The Hilbert uses a time-domain filter. The covariance matrix of the LCMV beamformer can utilize a time domain filter. The Morlet wavelet uses a wavelet. The CSD matrix of the DICS beamformer uses a wavelet. Thus, staying "within" domain makes the beamformer filter and the signal you apply the filter to match way better.

I am still opposed to applying an LCMV beamformer to a Morlet wavelet resolved TFR. Of course, you could show this actually works way better than the DICS beamformer on a Morlet TFR - but that would be a whole project (you cannot just show it on one data set and be done but need carefully simulated data etc).

Until then: let's not confuse our users.

[alexrockhill]

1. Hilbert outputs the same time-resolution as you put in, thus I sort this into a time-domain signal. The Morlet wavelet approach mostly does not do this - you traditionally have not the same time steps as in your input signal (but fewer).

That's not my understanding, here is the API, https://mne.tools/stable/generated/mne.time_frequency.tfr_array_morlet.html, you can decimate but by default there are the same number of time steps. And you can decimate the Hilbert in the same way.

2. And more importantly: The Hilbert uses a time-domain filter. The covariance matrix of the LCMV beamformer can utilize a time domain filter. The Morlet wavelet uses a wavelet. The CSD matrix of the DICS beamformer uses a wavelet. Thus, staying "within" domain makes the beamformer filter and the signal you apply the filter to match way better.

I agree that the CSD matrix uses a wavelet and thus matching CSD methods with your filter and and time-series object that you apply it to seems like the only logically consistent approach but I don't know what this has to do with LCMV, that seems like a misleading analogy. The Hilbert transform itself doesn't use a time-domain filter, you first filter your data before taking the analytic representation of the data which is very similar to a Fourier transform. Both Hilbert and Morlet methods generally window data in frequency (via the filter before Hilbert) and time-frequency respectively and do a Fourier transform or Fourier-transform-like operation to compute the amplitude or power so I don't see the big difference.

Happy to talk more in person, I do think by iterating we've clarified where each of our thought processes are coming from but it seems like I haven't done a good job convincing you so maybe we should see if other people weighing in helps us get on the same page.

EDIT: Maybe @wmvanvliet might have a minute to chime in here since you've done a lot of the DICS development. The question is: Is it okay to use Morlet wavelets instead of narrow bandpass filtered epoch data with a Hilbert transform as input to an LCMV beamformer? @britta-wstnr 's and my points of view is summarized above, if you have a few minutes.

[alexrockhill]

Here's the artifact @agramfort https://output.circle-artifacts.com/output/job/10323476-dfd6-4c9a-86c0-58228c1aa59f/artifacts/0/dev/auto_tutorials/inverse/55_tfr_stc.html?highlight=55_tfr, sorry I know it's late your time. It doesn't have many visuals because that has to be done when the GUI is built but basically:

1. You narrowband bandpass filter the data and apply the Hilbert transform to make a numpy array and then use that to construct an EpochsTFR object from scratch, then you compute LCMV filters and apply them which are now returned as an stc object with dimensions (n_epochs, n_verts, n_orient, n_freqs, n_times).

2. You compute cross-spectral density on sliding time windows to get the same stc object with DICS but this time the times are the center of the windows so there are fewer of them but the frequency resolution is really good.

3. You do the same thing as LCMV but treat the time-frequency data just like time series data and apply standard MNE looping over frequencies.

[alexrockhill]

Fixes https://github.com/mne-tools/mne-python/issues/10721.

I only mentioned passingly in a comment to a review from @britta-wstnr that I would like to keep the epochs dimension for my own analysis at single-trial resolution if that's not too much of a headache for everyone else.

Also, I didn't inherit TFRSourceEstimate from _BaseSourceEstimate yet because I thought it might be nice to use the newer info structure and inherit the mixins such as those being refactored here https://github.com/mne-tools/mne-python/pull/10945. I'm happy to inherit if that's too iconoclastic but I thought it might be a good time to start that kind of refactoring and build everything around mne.viz.Brain and the GUI abstraction https://github.com/mne-tools/mne-python/pull/10913. I won't be able to refactor all the other MNE SourceEstimates this summer but I think it's not too redundant to maintain during a transition period.

[alexrockhill]

@agramfort and @britta-wstnr , I get the following stack trace when I try to use complex phase-amplitude and not power:

ValueError                                Traceback (most recent call last)
<ipython-input-2-a143c40a3035> in <cell line: 11>()
     14     data_rank.update({ch_type: ch_rank for ch_type, ch_rank in
     15                       noise_rank.items() if ch_rank < data_rank[ch_type]})
---> 16     filters.append(mne.beamformer.make_lcmv(
     17         epochs_tfr.info, forward, data_cov, reg=0.05,
     18         noise_cov=noise_cov, pick_ori='vector',

<decorator-gen-497> in make_lcmv(info, forward, data_cov, reg, noise_cov, label, pick_ori, rank, weight_norm, reduce_rank, depth, inversion, verbose)

~/projects/mne-python/mne/beamformer/_lcmv.py in make_lcmv(***failed resolving arguments***)
    174     # compute spatial filter
    175     n_orient = 3 if is_free_ori else 1
--> 176     W, max_power_ori = _compute_beamformer(
    177         G, Cm, reg, n_orient, weight_norm, pick_ori, reduce_rank, rank_int,
    178         inversion=inversion, nn=nn, orient_std=orient_std,

~/projects/mne-python/mne/beamformer/_compute_beamformer.py in _compute_beamformer(***failed resolving arguments***)
    310     # Here W is W_ug, i.e.:
    311     # G.T @ Cm_inv / (G.T @ Cm_inv @ G)
--> 312     bf_denom_inv = _sym_inv_sm(bf_denom, reduce_rank, inversion, sk)
    313     assert bf_denom_inv.shape == (n_sources, n_orient, n_orient)
    314     W = np.matmul(bf_denom_inv, bf_numer)

~/projects/mne-python/mne/beamformer/_compute_beamformer.py in _sym_inv_sm(x, reduce_rank, inversion, sk)
    120         assert x.shape[1:] == (3, 3)
    121         if inversion == 'matrix':
--> 122             x_inv = _sym_mat_pow(x, -1, reduce_rank=reduce_rank)
    123             # Reapply source covariance after inversion
    124             x_inv *= sk[:, :, np.newaxis]

~/projects/mne-python/mne/utils/linalg.py in _sym_mat_pow(A, power, rcond, reduce_rank, return_s)
    170     limit = s[..., -1:] * rcond
    171     if not (s >= -limit).all():  # allow some tiny small negative ones
--> 172         raise ValueError('Matrix is not positive semi-definite')
    173     s[s <= limit] = np.inf if power < 0 else 0
    174     if reduce_rank:

ValueError: Matrix is not positive semi-definite

Any thoughts on how to fix this?

EDIT: I'm pushing a hacked version for now that uses single dipole inversion and suppressed an error so that you can see the changes but it definitely needs to be fixed.

The minimum norm source estimation worked with no errors also by the way.

[alexrockhill]

Ok, I couldn't figure out why the eigenvalues were not all positive but it has to do with the data covariance matrix which is positive definite for the CSD and is complex as well in that case so complex numbers are supported, it's just an issue with the conditioning of the matrix. The CSD matrix had a much larger diagonal than off diagonal entries for the real part compared to the covariance of the Hilbert which had larger off the diagonal than on the diagonal. I'm not exactly sure what the issue is with the conditioning though, I'm a bit stuck.

As I'm writing this, I'm thinking that real part of the covariance of the Hilbert should be positive since it's amplitude and I didn't check that this was true, I'll look into that next.

[agramfort]

hum hard to say without debugging. I think the code uses eigh which should also work with a complex valued matrix. Did you see how negative the eigenvalues are? Can you check at least that the input matrix is hermitian is A == A.H ?

[wmvanvliet]

Maybe @wmvanvliet might have a minute to chime in here since you've done a lot of the DICS development. The question is: Is it okay to use Morlet wavelets instead of narrow bandpass filtered epoch data with a Hilbert transform as input to an LCMV beamformer?

DICS == LCMV applied to a TFR (wavelets, multitaper, hilbert, whatever). A "tfr covariance" == CSD. So if you do LCMV using a "covariance matrix" derived from some sort of time-frequency decomposition, you are doing DICS. Literally, DICS and LCMV use the same code path (_compute_beamformer) with the only difference being s/cov/csd. So, it makes no sense to have an example that does this "manually" without using the DICS code.

What you are doing here is event-related DICS (Laaksonen, Kujala and Salmelin 2007) which we never got around to implementing. It looks like you are doing this with the apply_dics_evoked_tfr function, which is greatly appreciated.

[alexrockhill]

That was my understanding but now that you say that maybe it doesn't make sense to have the sliding windows but rather to only use the compute_covariance function on the epochs_tfr.

[alexrockhill]

hum hard to say without debugging. I think the code uses eigh which should also work with a complex valued matrix. Did you see how negative the eigenvalues are? Can you check at least that the input matrix is hermitian is A == A.H ?

They are large and negative, the largest was -13000 or so. I'll have to check on Hermitian.

[wmvanvliet]

now that you say that maybe it doesn't make sense to have the sliding windows

I think it is best to follow published methods for the main MNE-Python package. So I recommend following whatever approach they used in the paper. A sliding window will always have to be applied somewhere. When computing the TFR, when computing the CSD, or when applying the beamformer... I guess in the end it doesn't matter where we apply it, or even if we apply multiple sliding windows.

[wmvanvliet]

the compute_covariance function on the epochs_tfr.

Could you please modify the compute_csd function to be able to operate on TFRs instead? Computing "covariance" on a TFR object means computing a CSD.

[britta-wstnr]

I agree that the CSD matrix uses a wavelet and thus matching CSD methods with your filter and and time-series object that you apply it to seems like the only logically consistent approach but I don't know what this has to do with LCMV, that seems like a misleading analogy.

The LCMV beamformer uses a covariance matrix based on time-domain data. The Hilbert beamformer, as it was defined by my former supervisor @sarangnemo, computes a covariance matrix based on filtered time-domain data and then uses this in the LCMV beamformer. If you would now use a DICS for the same operation instead, you would need a CSD matrix as input, e.g. obtained via Morlet wavelets (as implemented in MNE), and the exact frequencies resolved would not match anymore. Especially so, since a Hilbert beamformer normally uses a broader frequency band (to achieve higher time resolution/less smearing in the time domain per the time-frequency tradeoff), which is hard to achieve with wavelets or Fourier (you'd need multitaper for that). That mismatch in frequencies resolved is what I am suspecting to be a really bad idea in theory and practice, what has not been done before in practice, and what I see as misleading for the user (as many will not be aware of the mismatch in frequency resolution).

I would vote for not making the focus to invent a new beamformer, but keep with what is used in the literature already - and for your GSOC just use synthetic input as suggested by @agramfort. The details of the computation of the beamformer should not influence the layout of the source space object. As mentioned before, in cases where it would make more sense to compute the spatial filter iteratively (Hilbert case, because you also have to do the filtering etc, which might make more sense to do on the fly in a loop instead of first creating all those epoch objects in different frequencies and with and without Hilbert applied), the output could still be transformed into a multi-dimensional source space object after the fact.

[britta-wstnr]

DICS == LCMV applied to a TFR (wavelets, multitaper, hilbert, whatever). A "tfr covariance" == CSD. So if you do LCMV using a "covariance matrix" derived from some sort of time-frequency decomposition, you are doing DICS. Literally, DICS and LCMV use the same code path (_compute_beamformer) with the only difference being s/cov/csd. So, it makes no sense to have an example that does this "manually" without using the DICS code.

What I would do is computed on CSD matrix across power for the whole time period - and then apply this time-step by time-step --- as you would with an LCMV beamformer that you apply to an evoked. For everything else we get very close to the 5-D beamformer by @sarangnemo, which we have just deprecated (as it is just not the best way of doing things anymore: your spatial filter suffers under not having many samples for your covariance/CSD matrix if you "slide" the make_filter operation across time.

[wmvanvliet]

Ok, so if I understand correctly, the Hilbert transform approach remains in the time domain (no complex values). So in that case you are doing LCMV. Clear!

[alexrockhill]

Ok, so if I understand correctly, the Hilbert transform approach remains in the time domain (no complex values). So in that case you are doing LCMV. Clear!

I think the idea was to use complex values and then take the power in source space, that's what's in the tutorial now. Since it's all using the same backend code in _beamformer.py is this just an issue of semantics? Or are you bringing up some technical issue with complex valued LCMV filters @wmvanvliet? because they did fail (causing Circle to fail in the CIs) because the cross-spectral density matrix was not positive semi-definite.

If you would now use a DICS for the same operation instead, you would need a CSD matrix as input, e.g. obtained via Morlet wavelets (as implemented in MNE), and the exact frequencies resolved would not match anymore. Especially so, since a Hilbert beamformer normally uses a broader frequency band (to achieve higher time resolution/less smearing in the time domain per the time-frequency tradeoff), which is hard to achieve with wavelets or Fourier (you'd need multitaper for that). That mismatch in frequencies resolved is what I am suspecting to be a really bad idea in theory and practice, what has not been done before in practice, and what I see as misleading for the user (as many will not be aware of the mismatch in frequency resolution).

We had a conversation about Morlet wavelets vs Hilbert as input to LCMV last Office Hours and the consensus seemed to be that they are different parameterizations of the same thing (i.e. by choosing the number of cycles you can get a very close approximation of a Hilbert transform of with a given bandwidth to the filter). Per @agramfort with @drammock and I agreeing. Again, I'm more than happy to use Hilbert as the example in the tutorial, but I haven't heard any persuasive argument why wavelets or Fourier transforms (multitaper) are a "bad idea" compared to using the Hilbert transform. I agree that there is more smearing in time with wavelets because they are complex-valued but the Hilbert transform has, in my opinion, the somewhat equivalent drawback that the bandwidth of the time-domain filter is can sub-optimally combine signals of interest by being too wide since you can't make a filter that perfectly selects one frequency without artifact. So, I think it's a recommendation-level issue and I agree with a slight recommendation for Hilbert, especially since it has already been published but we would have to go out of our way to not allow wavelets and I don't think we should do that.

[drammock]

We had a conversation about Morlet wavelets vs Hilbert as input to LCMV last Office Hours and the consensus seemed to be that they are different parameterizations of the same thing (i.e. by choosing the number of cycles you can get a very close approximation of a Hilbert transform of with a given bandwidth to the filter). Per @agramfort with @drammock and I agreeing.

For the record, I was present for that conversation but I was working on other stuff and not really paying close attention. I don't recall agreeing to anything.

[alexrockhill]

For the record, I was present for that conversation but I was working on other stuff and not really paying close attention. I don't recall agreeing to anything.

That's reasonable, it was mostly a conversation between Alex and I. You are welcome to weigh in one way or the other.

[alexrockhill]

hum hard to say without debugging. I think the code uses eigh which should also work with a complex valued matrix. Did you see how negative the eigenvalues are? Can you check at least that the input matrix is hermitian is A == A.H ?

Here are the first ten values of s returned from `eigh':

array([-816.53435783,  -63.88160429,  -37.18014888,  -34.082529  ,
        -27.26136952,  -21.52558745,  -17.72465881,  -13.89060314,
        -11.40168105,   -9.08428419])

And (np.matrix(Cm).H == Cm).all() returns True so it is Hermitian.

These lines convert it to Hermitian if there was rounding error:

# Whiten the data covariance
Cm = whitener @ Cm @ whitener.T.conj()
# Restore to properly Hermitian as large whitening coefs can have bad
# rounding error
Cm[:] = (Cm + Cm.T.conj()) / 2.

Interestingly, if I don't apply the baseline, I get a warning from computing the covariance matrix that it might be inaccurate and then I get another warning that the matrix is not positive semi-definite but I don't get the error that it isn't positive semi-definite. Not exactly sure what to make of that. Also, I tried it with Morlet wavelets and it gives the same error. I'm really not sure what the root cause is...

[agramfort]

I am far from my keyboard these days. I cannot look. Just for me Hilbert gives you the analytic signal which is complex values. When you take the abs of it you get the envelope. Now to me calling lcmv on the envelope is wrong due to nonlinear abs function