ENH: Power line noise filtering - spectrum interpolation and zapline

[DominiqueMakowski]

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

• spectrum interpolation, implemented in fieldtrip, seems to have a rather straightforward definition.
• ZapLine is also implemented in matlab.

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

References:

• Leske, S., & Dalal, S. S. (2019). Reducing power line noise in EEG and MEG data via spectrum interpolation. NeuroImage, 189, 763-776.
• de Cheveigné, A. (2019). ZapLine: a simple and effective method to remove power line artifacts. NeuroImage, 116356.

[nbara]

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]

@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]

@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 :

• 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)
• Zapline to remove 50Hz : useful and easy + minimal overhead since DSS was already implemented by @larsoner here
• Robust detrending : very useful, maked high-pass filtering obsolete for a lot of applications
• CCA and multiway CCA, as a (more performant, according to the authors) alternative to ReceptiveField : very useful, but would require more work to be coded to MNE-s standards :)

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]

@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 #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?

- Zapline to remove 50Hz : useful and easy + minimal overhead since DSS was already implemented by @larsoner <https://github.com/larsoner> here <mne-tools/mne-incubator#30> we have line noise in many mne datasets so maybe it can be considered? - Robust detrending <https://www.sciencedirect.com/science/article/pii/S1053811918300351#bib32> : very useful, maked high-pass filtering obsolete for a lot of applications I've heard good things about this too. How would you illustrate the

benefit in an example?

- CCA <https://doi.org/10.1016/j.neuroimage.2018.01.033> and multiway CCA <https://doi.org/10.1016/j.neuroimage.2018.11.026>, as a (more performant, according to the authors) alternative to ReceptiveField : very useful, but would require more work to be coded to MNE-s standards :) 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]

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 :

• DSS (+ Zapline): could be easily applied like in this tutorial
• Robust detrending (high-pass filtering is a big issue for ERP research, so this is important); there's at least one independent paper that I know of showing benefits of detrending on MVPA; maybe @kingjr will want to chime in on this because it would potentially be helpful for his temporal generalization analyses
• STAR (can be useful)
• SNS (maybe)
• TSPCA (probably not)

[larsoner]

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]

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]

Your plan sounds good, except that:

• 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.
• Instead of having one that works on plain covariances I would just go to Epochs. For example that is how XDawn works: https://mne.tools/stable/generated/mne.preprocessing.Xdawn.html

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]

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]

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]

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]

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]

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

[britta-wstnr]

@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]

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

[britta-wstnr]

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]

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

[Lubell]

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

[agramfort]

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

[Lubell]

I did not write the original code

[larsoner]

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]

Ahh we need to talk to the original author then

[sarangnemo]

@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]

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

[britta-wstnr]

@larsoner would that work?

[larsoner]

Seems to me like that should be sufficient

[britta-wstnr]

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

[larsoner]

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]

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

[hoechenberger]

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

[jnvandermeer]

I second that!

[larsoner]

@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]

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

[cbrnr]

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

[chapochn]

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).