search

mne-tools / mne-python

MNE: Magnetoencephalography (MEG) and Electroencephalography (EEG) in Python
https://mne.tools
BSD 3-Clause "New" or "Revised" License
2.66k stars 1.31k forks source link

ENH: Power line noise filtering - spectrum interpolation and zapline #7170

Open DominiqueMakowski opened 4 years ago

DominiqueMakowski commented 4 years ago

Two "new" methods have been recently published that are supposedly superior to notch filtering for power line noise removal.

These might interesting additions for MNE. As my matlab skills are not enough for the task, I'm leaving this here for reference ☺️

References:

nbara commented 4 years ago

Hi Dominique,

I've implemented ZapLine into my MEEGkit toolbox for python :

https://github.com/nbara/python-meegkit/blob/master/meegkit/dss.py#L126

Let me know if you have any comments or need some help with it.

agramfort commented 4 years ago

@nbara do you have rendered doc of this package? I discover it now... should we mutualize a bit the efforts between this code and mne/preprocessing ? do you have examples using real data? why not using sphinx-gallery for the examples?

nbara commented 4 years ago

@agramfort

do you have rendered doc of this package? I discover it now... [...] why not using sphinx-gallery for the examples?

I've never taken the time to properly set up sphinx-gallery, or upload the doc to readthedocs...

should we mutualize a bit the efforts between this code and mne/preprocessing ?

Well the reasons I created the package alone in the first place were a) as an practical exercice to learn python b) because the DSS code had been stuck in MNE-sandbox for several years, so I assumed there was no interest.

There are actually quite a lot of things on MEEGkit that could be integrated into MNE :

All in all, I'd be happy to help integrate some of this into MNE, potentially during a sprint. But it depends on what people actually find useful.

do you have examples using real data?

I use it all the time on my data (EEG now, MEG in the past). Pretty sure it could easily fit the example datasets in MNE (e.g. the audiovisual dataset).

agramfort commented 4 years ago

@agramfort https://github.com/agramfort

do you have rendered doc of this package? I discover it now... [...] why not using sphinx-gallery for the examples?

I've never taken the time to properly set up sphinx-gallery, or upload the doc to readthedocs...

sphinx-gallery and read the docs is not necessarily easy to work together

having circle ci be able to push to the gh-pages branch is easier.

let us know if we can help

should we mutualize a bit the efforts between this code and mne/preprocessing ?

Well the reasons I created the package alone in the first place were a) as an practical exercice to learn python b) because the DSS code had been stuck in MNE-sandbox for several years, so I assumed there was no interest.

reading https://github.com/mne-tools/mne-python/issues/2112 again I see that the conclusion when looking at empirical data was not super clear. If we stick to the rule that code should enter mne when it comes with a convincing example on real data that explains why it was moved to sandbox. Do we need to have another dataset in mne?

There are actually quite a lot of things on MEEGkit that could be integrated into MNE :

  • TSPCA (I use it less now as it's mostly useful for MEG)
  • SNS (not that useful for me in practice)
  • STAR : to remove sparse time-artifacts (may be a bit redundant with autoreject)

maybe it can be used to spot and fix flux jumps? you have flux jumps in MNE sample data so maybe it can be a convincing example?

we have line noise in many mne datasets so maybe it can be considered?

I've heard good things about this too. How would you illustrate the benefit in an example?

one thing at a time :)

All in all, I'd be happy to help integrate some of this into MNE, potentially during a sprint. But it depends on what people actually find useful.

do you have examples using real data?

I use it all the time on my data (EEG now, MEG in the past). Pretty sure it could easily fit the example datasets in MNE (e.g. the audiovisual dataset).

I would like to push you one more time to provide examples that use real data and not simulations

nbara commented 4 years ago

reading #2112 again I see that the conclusion when looking at empirical data was not super clear. If we stick to the rule that code should enter mne when it comes with a convincing example on real data that explains why it was moved to sandbox. Do we need to have another dataset in mne?

That issue you refer to mostly deals with this SSS method (not familiar with it), and AFAICT @kingjr wasn't able to get the TSPCA code to run ?

Either way, regarding TSCPA, I don't use it at all as I don't do MEG anymore.

On the other hand I think DSS is super versatile (can be used to remove line noise / ZapLine, enhance evoked activity, maximise difference between two experimental conditions, etc.) and easy to implement. It has been used in numerous publications (including mine, bot on EEG and MEG), so finding working examples would be easy.

So if I had to rank these methods in order of priority :

larsoner commented 4 years ago

Agreed, the first 2-4 seem reasonable to me.

It looks like my PR was for tsDSS, do we need that or just plain DSS? In any case, want to make a PR to add DSS to start or, if you're not sure about API, open an issue for API discussion for it?

nbara commented 4 years ago

It looks like my PR was for tsDSS, do we need that or just plain DSS?

just plain DSS.

In terms of API, I think it would make sense to have a base DSS estimator, and children classes like: 

class DSS(BaseEstimator):  # base class, works on arbitrary covariances  
    ...

class EvokedDSS(DSS):  # enhance evoked activity
    ...

class tsDSS(DSS):  # time-shift DSS
    ...

class ZapLine(DSS):  # removes line noise
    ...

? And yes we can open an issue to discuss API.

Robust detrending is a separate issue IMO.

larsoner commented 4 years ago

Your plan sounds good, except that:

Want to take my code (or yours if it's better) from #30 and give it a shot?

As for the ZapLine class I'm not sure, it depends on what ZapLine actually needs to do to operate. I was naively hoping we could just add a method='zapline' to mne.filter.notch_filter. But let's do that after (ts)DSS is in.

nbara commented 4 years ago

I wouldn't separate tsDSS from DSS since the latter is just a case of the former. Since the more common name is DSS I would probably call the class DSS.

Sure that's fine.

Instead of having one that works on plain covariances I would just go to Epochs. For example that is how XDawn works

Well in principle DSS doesn't have to be used on epoched data... It depends on your objective function.

For me the most basic use of DSS (and one that makes it really versatile) is to provide a base covariance (c0 in the original paper nomenclature), as well as a biased covariance c1.

c1 can be anything, and we can't provide an estimator that takes data in the temporal form (n_channels, n_times) to fit all the use cases.

For instance, a nice use case is to compute c1 as the covariance of the difference between two experimental condition, then DSS finds the components which are maximally different between those two conditions.

As for the ZapLine class I'm not sure, it depends on what ZapLine actually needs to do to operate. I was naively hoping we could just add a method='zapline' to mne.filter.notch_filter. But let's do that after (ts)DSS is in.

ZapLine is simply a special case of DSS, where c1 is the covariance of the data filtered with a comb-shaped filter @ 50Hz + harmonics. This finds the most line-dominated components.

A final step consists in regressing out these line-dominated components from the original data. It's nice as it doesn't create any 'holes' in the spectrum, and doesn't cause long edge artifacts.

larsoner commented 4 years ago

For me the most basic use of DSS (and one that makes it really versatile) is to provide a base covariance (c0 in the original paper nomenclature), as well as a biased covariance c1.

How about the c0 gets passed to the constructor, and you pass c1 to the fit method then? And apply operates on instances like Raw / Evoked / Epochs?

class DSS(object)
    def __init__(self, c0, max_delay=0)
        ...

    def fit(self, c1)
        ...

    def apply(self, inst)
        ...

?

agramfort commented 4 years ago

I would start by robust detrending as I think this can be well illustrated with examples.

If I can suggest something it is to aim for small incremental contributions. Nothing too big in one go.

nbara commented 4 years ago

How about the c0 gets passed to the constructor, and you pass c1 to the fit method then? And apply operates on instances like Raw / Evoked / Epochs?

class DSS(object)
    def __init__(self, c0, max_delay=0)
        ...

    def fit(self, c1)
        ...

    def apply(self, inst)
        ...

?

Why would you provide c0 in the constructor ? It's a bit weird no? I would provide c0 and c1 in fit(X=c0, y=c1) ...

I would start by robust detrending as I think this can be well illustrated with examples. If I can suggest something it is to aim for small incremental contributions. Nothing too big in one go.

Yes I agree. I will try to put up a PR for robust detrending soon.

DominiqueMakowski commented 4 years ago

I'm glad I could connect you guys ☺️ really looking forward to using these new methods!

britta-wstnr commented 4 years ago

@agramfort @nbara Do your previous efforts include the spectrum interpolation (Leske & Dalal 2019)? @lubell who is a PhD student in @sarangnemo's lab has a Python implementation of this and we thought about looking into contributing it if this is still relevant!

Here is the code: https://github.com/Lubell/python_dftfilter

agramfort commented 4 years ago

I would only suggest to compare. I don't have any experience myself.

britta-wstnr commented 4 years ago

Our implementation is the spectrum interpolation as implemented in FieldTrip that @DominiqueMakowski mentioned in his initial post, developed by Sabine Leske and @sarangnemo, tackling line noise. Is anything like this already implemented now? What should we compare it to?

larsoner commented 4 years ago

That code is GPL so if we do want to add it we need permission to relicense under BSD

Lubell commented 4 years ago

@larsoner Happy to remove the GPL, just part of my routine to include it. Would that be enough to remove re-license obstacles?

agramfort commented 4 years ago

did you write the original code in fieldtrip? is the corresponding fieldtrip GPL?

Lubell commented 4 years ago

I did not write the original code

larsoner commented 4 years ago

I think that some explicit permission from the code author here in our conversation should be enough. Do we have permission?

(Also from what I recall, just removing the license would not fix the problem as code with no explicit license by default gets a pretty restrictive one. You'd have to instead replace with a permissive BSD-compatible license like BSD, MIT, Apache, etc.)

larsoner commented 4 years ago

Ahh we need to talk to the original author then

sarangnemo commented 4 years ago

@sabineleske and I wrote the original FieldTrip code, we didn't specify a special license so I suppose it's GPL like FieldTrip in general. You have my permission, @sabineleske are you there?

sabineleske commented 4 years ago

Hi! yes, of course this is fine with me. You have my permission as well.

britta-wstnr commented 4 years ago

@larsoner would that work?

larsoner commented 4 years ago

Seems to me like that should be sufficient

britta-wstnr commented 4 years ago

@Lubell Let's get started? @larsoner Should we create a new issue to discuss the implementation details before Jamie starts on a PR?

larsoner commented 4 years ago

My first thought would be to add it as a new method in mne.filter.notch_filter. If it's not clear that this is a good choice or how to do it then we can discuss in another issue

britta-wstnr commented 4 years ago

@larsoner, that sounds reasonable. I will get together with @Lubell and we'll get started on this!

hoechenberger commented 4 years ago

Just wanted to chime in and and say that I believe ZapLine would be an extremely very useful addition to MNE.

jnvandermeer commented 4 years ago

I second that!

larsoner commented 4 years ago

@britta-wstnr maybe a good one for when #7085 is merged. Or if you have other things to work on maybe @hoechenberger is interested in working on it with @Lubell ?

adswa commented 2 years ago

Also just wanted to chime in and say that I would find it very useful if a ZAPline implementation exists in MNE :)

cbrnr commented 2 years ago

I think mne.filter.notch_filter() with method='spectrum_fit' should already do something like that.

chapochn commented 1 year ago

In my understanding method='spectrum_fit' and ZapLine are very different. method='spectrum_fit' subtracts cos/sin functions from each individual electrode, whereas ZapLine does DSS on all the electrodes to find the subspace that contains the powerline and then removes it from all the electrodes that have power in that subspace. Zapline works on all electrodes at the same time and finds the linear combination of electrodes that maximizes the 60Hz power (and harmonics).