[MRG+2] ENH: padding kwarg for psd_welch by mmagnuski · Pull Request #4003 · mne-tools/mne-python

mmagnuski · 2017-02-18T20:17:06Z

Fixes #3621
by adding padding kwarg to psd_welch.
I'll add an entry to whats_new once comments are resolved and CIs green.

codecov-io · 2017-02-18T21:59:50Z

Codecov Report

Merging #4003 into master will decrease coverage by -12.12%.
The diff coverage is 100%.

@@             Coverage Diff             @@
##           master    #4003       +/-   ##
===========================================
- Coverage   86.06%   73.95%   -12.12%     
===========================================
  Files         349      350        +1     
  Lines       62680    62999      +319     
  Branches     9591     9631       +40     
===========================================
- Hits        53945    46588     -7357     
- Misses       6040    13617     +7577     
- Partials     2695     2794       +99

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 1182668...fb5db5f. Read the comment docs.

agramfort · 2017-02-19T08:24:06Z

    """Aux function."""
    return func(epoch, fs=fs, nperseg=nfft, noverlap=noverlap,
-                nfft=nfft, window='hann')[2][..., freq_mask, :]
+                nfft=padding, window='hann')[2][..., freq_mask, :]


I don't like this. nfft and padding have different meanings.

This is the scipy.spectral.spectrogram API. It differs from mne's psd_welch: nperseg is the length of each welch window and nfft is the final length of each window. So if nfft > nperseg it results in zero padding nperseg segments to nfft. psd_welch already uses nfft to mean nperseg, so I added padding (or pad_to or a different kwarg name if you prefer). We could also change psd_welch a little to use the same logic as scipy spectrogram (nperseg and nfft) but that would require deprecation (although the behavior could stay backward compatible - the meaning of the parameters would change).

@mmagnuski can you summarize the new logic? You add a padding parameter which can override the n_fft if it is not None?

Previously both nperseg and nfft were populated with psd_welch nfft - which means that padding was never used. Now there is a separate kwarg for padding that is used only when not None and > nfft. I agree that naming is unfortunate - ideally nfft should control padding if > nperseg (just like in scipy). We can have such behavior in mne but this would probably require a deprecation...

mmagnuski · 2017-02-19T11:59:40Z

Or not - depending on the wording of the docs and your opinion. :)
We can pass nfft to scipys nfft and add nperseg kwarg.

dengemann · 2017-02-19T12:18:38Z

@mmagnuski I'm trying to better understand a) the motivation for this b) how it changes the exisiting behavior with regard to n_fft. Maybe you want to share your insights :)

mmagnuski · 2017-02-19T13:39:43Z

Sure, the motivation for this is adding control over padding - see #3621. It wasn't a long discussion but the general agreement was to add control over padding, likely through additional kwarg. Currently you can't have padding - only segment length is controlled.
Segment length is termed n_fft in mne and nperseg in scipy. Mne naming seems to stem from the fact that there was no ability to have zero-padding, so segment length was the same as n_fft.

Introducing control over padding requires either:

adding nperseg kwarg that defines segment length, and make n_fft control padding when n_fft > nperseg. But because previously n_fft was closer in meaning to nperseg it was adjusted down to len(inst.times) when n_fft > len(inst.times). So this kind of change would rather require deprecation.
another option is to add padding or pad_to, as done in this PR. A separate kwarg that controls padding, while n_fft controls segment length. Padding takes place only when padding > n_fft. The downside here is that n_fft works as nperseg which may be conterintuitive.

I also cover the logic proposed here in my previous comments - in case you didn't see all of them.
And I actually would prefer adding nperseg instead of padding - but this requires some kind of deprecation (n_fft >signal_length would do padding).

larsoner · 2017-02-19T13:45:33Z

I'm okay with deprecation here if it gets equivalent API with SciPy.

agramfort · 2017-02-19T17:12:25Z

what disturb me is that nfft = len(signal) + padding it was a semantic mismatch. Why would you want to do some padding in the first place? Have non power of 2 segments but fast ffts?

mmagnuski · 2017-02-19T18:07:34Z

@agramfort
Yeah, I understand the problem with semantic mismatch. What about solution 2?

It is not so much my need to add padding - but in #3621 you agreed it should be added (at least I read your "I agree" this way :) ). In general non power of 2 is one case but interpolation is another. I remember I was a bit surprised when I started using mne and found out that you can't set window length and padding - just as the person who submitted

mmagnuski · 2017-02-19T18:09:17Z

sorry, I'm on my phone, pushed wrong button by accident.
...the person who submitted #3621.

agramfort · 2017-02-19T20:21:27Z

please show me how high level API would look like

mmagnuski · 2017-02-19T20:38:45Z

with solution 2 that would be:

psds, freqs = psd_welch(raw, nperseg=128, n_fft=256)

with solution 1:

psds, freqs = psd_welch(raw, n_fft=128, pad_to=256)

agramfort · 2017-02-20T09:08:25Z

ok fine with solution 2 with nperseg defaulting to None meaning use n_fft so we have no deprecation.

mmagnuski · 2017-02-20T09:24:03Z

n_fft > signal length would do zero-padding, while previously it would be set back to signal length. So there should be some deprecation info if nperseg=None and n_fft > signal length, right? Or do you suggest that n_fft >

mmagnuski · 2017-02-20T09:25:35Z

(eh, same problem, sorry)
n_fft > signal length should not do zero padding when nperseg=None?

agramfort · 2017-02-20T09:55:32Z

well calling welch on a signal that is less than nperseg is not really doing welch and should probably be not done. You should use mulit-tapers for such short signals

larsoner · 2017-02-20T16:01:41Z

well calling welch on a signal that is less than nperseg is not really doing welch and should probably be not done.

Agreed. I think it's okay to raise an error if someone tries to do this. If it breaks some existing code and forces someone who really wants to do it to now be explicit about it, I think that's okay

mmagnuski · 2017-02-20T16:27:39Z

Ok, so if nperseg=None and n_fft > signal_length -> raise an error stating that if one really wants this nperseg has to specified.

larsoner · 2017-02-20T16:29:03Z

Sounds reasonable to me. This allows us to proceed with solution 2 (SciPy API), right?

agramfort · 2017-02-20T16:43:26Z

yes

mmagnuski · 2017-02-21T15:15:38Z

I've added nperseg kwarg, but everywhere in mne there is rather n_something than nsomething. For example - psd_welch has n_fft, not nfft. So n_per_seg might be more consistent with the rest of naming (although I don't like two underscores there). I've also changed the docstring a little, so take a look if you like the changes.

larsoner

I guess since we already have different underscores we might as well have n_per_seg. At least the meaning immediately translates to scipy

mmagnuski · 2017-02-21T16:10:41Z

Just in case another option could be segment_len (or maybe n_segment) or window_len (etc.)

larsoner · 2017-02-21T16:16:34Z

I think it's better to stay close to SciPy. n_per_seg seems like the right compromise between equivalence with SciPy (which we can't get fully anyway anymore) and our internal style practices

mmagnuski · 2017-02-21T16:21:45Z

Ok, I'll change it to n_per_seg then.

larsoner · 2017-02-21T16:45:42Z

Ready for review/merge from your end?

mmagnuski · 2017-02-21T18:05:34Z

Yes, it should be good for review/merge now.

larsoner

Otherwise LGTM

larsoner · 2017-02-21T20:17:36Z

+    """Helper to make sure n_fft, n_per_seg and n_overlap make sense."""
+    n_per_seg = n_fft if n_per_seg is None or n_per_seg > n_fft else n_per_seg
+    n_per_seg = n if n_per_seg > n else n_per_seg
+    n_overlap = n_per_seg - 1 if n_overlap >= n_per_seg else n_overlap


Maybe we covered this already, but this should probably just be an error, no?

that's the previous behavior, I can change it to an error if you prefer.

I'd rather have people be explicit if they want this behavior. This looks like it's silent-errant-behavior prone to me

larsoner · 2017-02-21T20:18:44Z

-        it is the smoother are the PSDs. The default value is 256.
-        If ``n_fft > len(inst.times)``, it will be adjusted down to
-        ``len(inst.times)``.
+        The length of FFT used. If n_per_seg is None n_fft sets the length of the


This is good enough:

The length of FFT used, must be ``>= n_per_seg`` (default: 256). The segments will be zero-padded if ``n_fft > n_per_seg``.

Yes, that's better, thanks!

larsoner · 2017-02-21T20:20:18Z

+    n_per_seg : int | None
+        Length of each Welch segment. The smaller it is with respect to the
+        signal length the smoother are the PSDs. Defaults to None, which sets
+        n_per_seg equal to n_fft. If n_per_seg is smaller than n_fft, each window


If you add the bits about padding above, they don't need to be here.

larsoner · 2017-02-21T20:22:06Z

+    if n_per_seg is None and n_fft > n_times:
+        raise ValueError('If n_per_seg is None n_fft is not allowed to be >'
+                         ' n_times. If you want zero-padding, you have to set'
+                         ' n_per_seg to relevant length.')


this check should go in _check_nfft

ok, I'll move it there

mmagnuski · 2017-02-22T12:56:46Z

Codecov is confused, it says "Absolute coverage decreased by -12.11% ..." and lists many tests as not tested...

mmagnuski · 2017-02-22T12:57:25Z

and Travis did not run

larsoner · 2017-02-22T14:56:06Z

@agramfort +1 for merge from you? looks like comments have been addressed

@mmagnuski sometimes there is a delay, everything looks good now. +1 for merge from me

larsoner · 2017-02-22T14:57:24Z

+    n_per_seg = n if n_per_seg > n else n_per_seg
+    if n_overlap >= n_per_seg:
+        msg = ('n_overlap cannot be greater than n_per_seg (or n_fft). Got '
+               'n_overlap of {} while n_per_seg is {}.')


FYI for future reference by convention we 1. don't split these across multiple vars/commands and 2. use/prefer % (even though there is talk of deprecation)

Thanks for stating what the preference is - I think I saw this using format and a separate var with text) a few times in Jona's code so I thought it is ok. I'll remember now. :) I can change it here too.

Sometimes .format is a lot more convenient (e.g., mne/datasets/utils.py I think?), and we do use it in those cases. But it's not a big deal, so don't bother unless someone else has changes they want you to make

agramfort · 2017-02-22T15:30:10Z

please update what's new and then +1 for MRG

agramfort · 2017-02-23T10:41:34Z

yes ok for me

mmagnuski · 2017-02-24T09:37:04Z

@Eric89GXL the errors we added make a test fail:

======================================================================
ERROR: Test creating raw from array.
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/travis/miniconda/envs/test/lib/python2.7/site-packages/nose/case.py", line 197, in runTest
    self.test(*self.arg)
  File "/home/travis/miniconda/envs/test/lib/python2.7/site-packages/numpy/testing/decorators.py", line 147, in skipper_func
    return f(*args, **kwargs)
  File "/home/travis/build/mne-tools/mne-python/mne/io/array/tests/test_array.py", line 95, in test_array_raw
    raw2.plot_psd(tmax=np.inf, average=True, spatial_colors=False)
  File "<string>", line 2, in plot_psd
  File "/home/travis/build/mne-tools/mne-python/mne/utils.py", line 709, in verbose
    return function(*args, **kwargs)
  File "/home/travis/build/mne-tools/mne-python/mne/io/base.py", line 1638, in plot_psd
    spatial_colors=spatial_colors)
  File "<string>", line 2, in plot_raw_psd
  File "/home/travis/build/mne-tools/mne-python/mne/utils.py", line 709, in verbose
    return function(*args, **kwargs)
  File "/home/travis/build/mne-tools/mne-python/mne/viz/raw.py", line 598, in plot_raw_psd
    n_overlap=n_overlap, n_jobs=n_jobs)
  File "<string>", line 2, in psd_welch
  File "/home/travis/build/mne-tools/mne-python/mne/utils.py", line 709, in verbose
    return function(*args, **kwargs)
  File "/home/travis/build/mne-tools/mne-python/mne/time_frequency/psd.py", line 200, in psd_welch
    n_jobs=n_jobs, verbose=verbose)
  File "<string>", line 2, in psd_array_welch
  File "/home/travis/build/mne-tools/mne-python/mne/utils.py", line 709, in verbose
    return function(*args, **kwargs)
  File "/home/travis/build/mne-tools/mne-python/mne/time_frequency/psd.py", line 114, in psd_array_welch
    n_overlap)
  File "/home/travis/build/mne-tools/mne-python/mne/time_frequency/psd.py", line 26, in _check_nfft
    ' signal length is %d.') % (n_fft, n))
ValueError: If n_per_seg is None n_fft is not allowed to be > n_times. If you want zero-padding, you have to set n_per_seg to relevant length. Got n_fft of 2048 while signal length is 1803.

Previous behavior was to set n_fft equal to signal length if n_fft > signal length so I'll change this test by setting n_fft to signal length.

mmagnuski · 2017-02-24T13:35:06Z

All builds are green except one:

======================================================================
FAIL: Test object size estimation.
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/travis/miniconda/envs/test/lib/python3.5/site-packages/nose/case.py", line 198, in runTest
    self.test(*self.arg)
  File "/home/travis/miniconda/envs/test/lib/python3.5/site-packages/mne-0.14.dev0-py3.5.egg/mne/tests/test_utils.py", line 111, in test_object_size
    msg='%s < %s < %s:\n%s' % (lower, size, upper, obj))
AssertionError: False is not true : 0 < 68 < 60:
foo
    'False is not true : 0 < 68 < 60:\nfoo' = self._formatMessage('False is not true : 0 < 68 < 60:\nfoo', "%s is not true" % safe_repr(False))
>>  raise self.failureException('False is not true : 0 < 68 < 60:\nfoo')

larsoner · 2017-02-24T15:14:34Z

Yep, failure is unrelated. Thanks @mmagnuski !

mmagnuski · 2017-02-24T15:32:57Z

Thanks @Eric89GXL @agramfort 🚀

larsoner · 2017-02-24T15:35:24Z

Onto #4000 :)

agramfort reviewed Feb 19, 2017

View reviewed changes

mmagnuski closed this Feb 19, 2017

mmagnuski reopened this Feb 19, 2017

mmagnuski closed this Feb 20, 2017

mmagnuski reopened this Feb 20, 2017

larsoner reviewed Feb 21, 2017

View reviewed changes

larsoner changed the title ~~[MRG] ENH: padding kwarg for psd_welch~~ [MRG+1] ENH: padding kwarg for psd_welch Feb 22, 2017

larsoner reviewed Feb 22, 2017

View reviewed changes

larsoner changed the title ~~[MRG+1] ENH: padding kwarg for psd_welch~~ [MRG+2] ENH: padding kwarg for psd_welch Feb 23, 2017

mmagnuski added 9 commits February 24, 2017 00:35

ENH: add padding

d6e9410

TST: add test for psd_welch padding param

acf1961

FIX: python 2 compatibility

7dd6e4a

STY: fix overindent

30a6e3f

ENH: nperseg kwarg instead of padding

d744236

rename nperseg -> n_per_seg

93c33d5

adress comments

2f38d23

change error messages to conform with mne style

880a20a

add whats new

48ac066

mmagnuski force-pushed the welch_padding branch from fb5db5f to 48ac066 Compare February 24, 2017 00:26

FIX: leftover from rebase conflict

c32c45c

FIX: fix failing test due to changes in how n_fft works

15cac5e

mmagnuski mentioned this pull request Feb 24, 2017

WIP: add reduction keyword to pwelch #4000

Closed

2 tasks

larsoner merged commit 9e4ac8a into mne-tools:master Feb 24, 2017

mmagnuski deleted the welch_padding branch February 24, 2017 15:38

CedricConday mentioned this pull request Jun 30, 2026

Fix psd_array_welch for good data spans shorter than n_per_seg (#13039) #14003

Open

Uh oh!

Uh oh!

Conversation

mmagnuski commented Feb 18, 2017 • edited by larsoner Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codecov-io commented Feb 18, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mmagnuski commented Feb 19, 2017

Uh oh!

dengemann commented Feb 19, 2017

Uh oh!

mmagnuski commented Feb 19, 2017

Uh oh!

larsoner commented Feb 19, 2017 via email

Uh oh!

agramfort commented Feb 19, 2017 via email

Uh oh!

mmagnuski commented Feb 19, 2017

Uh oh!

mmagnuski commented Feb 19, 2017

Uh oh!

agramfort commented Feb 19, 2017 via email

Uh oh!

mmagnuski commented Feb 19, 2017

Uh oh!

agramfort commented Feb 20, 2017 via email

Uh oh!

mmagnuski commented Feb 20, 2017

Uh oh!

mmagnuski commented Feb 20, 2017

Uh oh!

agramfort commented Feb 20, 2017 via email

Uh oh!

larsoner commented Feb 20, 2017

Uh oh!

mmagnuski commented Feb 20, 2017

Uh oh!

larsoner commented Feb 20, 2017

Uh oh!

agramfort commented Feb 20, 2017 via email

Uh oh!

mmagnuski commented Feb 21, 2017

Uh oh!

larsoner left a comment

Choose a reason for hiding this comment

Uh oh!

mmagnuski commented Feb 21, 2017

Uh oh!

larsoner commented Feb 21, 2017

Uh oh!

mmagnuski commented Feb 21, 2017

Uh oh!

larsoner commented Feb 21, 2017

Uh oh!

mmagnuski commented Feb 21, 2017

Uh oh!

larsoner left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mmagnuski commented Feb 18, 2017 •

edited by larsoner

Loading

codecov-io commented Feb 18, 2017 •

edited

Loading