# [MRG] VIZ/DOC/ENH: add time frequency gfp example #4139

Merged
merged 20 commits into from Apr 1, 2017

## Conversation

Projects
None yet
5 participants
Member

### dengemann commented Mar 30, 2017 • edited Edited 1 time dengemann edited Mar 31, 2017 (most recent)

 This is using some old school methods in new school ways to explore spectral temporal evolution. The (latest) output looks like this: cc @Eric89GXL @agramfort @jona-sassenhagen
Member

### dengemann commented Mar 30, 2017

 Advanced example of a better way of doing #3123.
Member

### dengemann commented Mar 30, 2017

 image output updated.

### agramfort reviewed Mar 30, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 subtract the avarage and finally rectify the signals prior to averaging. Then the GFP is computed as described in [2], using the sum of the squares devided by the true degrees of freedom. Baselinging is

Member

Baselining

### jona-sassenhagen reviewed Mar 30, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 # Then we prepare a bootstrapping function to estimate confidence intervals def get_gfp_ci(average, rank=rank):

#### jona-sassenhagen Mar 30, 2017

Contributor

Nice, will have to use this one for plot_compare_evokeds :)

Member

:)

Contributor

:D

Contributor

### jona-sassenhagen commented Mar 30, 2017

 Not much time now, but it looks pretty cool at a glance!

### agramfort reviewed Mar 30, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 Then the GFP is computed as described in [2], using the sum of the squares devided by the true degrees of freedom. Baselinging is applied to make the GFPs comparabel between frequencies.

#### agramfort Mar 30, 2017

Member

comparable

examples/time_frequency/plot_time_frequency_global_field_power.py
 applied to make the GFPs comparabel between frequencies. The procedure is then repeated for each frequency band of interest and all GFPs are visualized. The non-parametric confidence intervals are computed as described in [3].

#### agramfort Mar 30, 2017

Member

computed using bootstrap as described

examples/time_frequency/plot_time_frequency_global_field_power.py
 The procedure is then repeated for each frequency band of interest and all GFPs are visualized. The non-parametric confidence intervals are computed as described in [3]. The advantage of this method over summarizing the the Space x Time x Frequency

#### agramfort Mar 30, 2017

Member

the the

examples/time_frequency/plot_time_frequency_global_field_power.py
 # set epoching parameters event_id, tmin, tmax = 1, -1., 3. baseline = None frequency_map = list()

#### agramfort Mar 30, 2017

Member

insert empty line

examples/time_frequency/plot_time_frequency_global_field_power.py
 baseline = None frequency_map = list() for band, fmin, fmax in iter_freqs:

#### agramfort Mar 30, 2017

Member

remove line

examples/time_frequency/plot_time_frequency_global_field_power.py
 frequency_map = list() for band, fmin, fmax in iter_freqs: raw = mne.io.read_raw_fif(raw_fname, preload=True)

#### agramfort Mar 30, 2017

Member

examples/time_frequency/plot_time_frequency_global_field_power.py
 # remove evoked response and get analytic signal (envelope) data = epochs.get_data() mean = np.mean(data, axis=0, keepdims=True) epochs._data = np.abs(data - mean)

#### agramfort Mar 30, 2017

Member

use epochs.subtract_evoked() method

examples/time_frequency/plot_time_frequency_global_field_power.py
 rank = raw.estimate_rank() print('The numerical rank is: %i' % rank) rng = np.random.RandomState(42) mne.utils.set_log_level('warning')

#### agramfort Mar 30, 2017

Member

don't put this in an example

# Codecov Report

Merging #4139 into master will increase coverage by 0.27%.
The diff coverage is n/a.

@@            Coverage Diff             @@
##           master    #4139      +/-   ##
==========================================
+ Coverage   86.18%   86.46%   +0.27%
==========================================
Files         354      276      -78
Lines       63739    60597    -3142
Branches     9709     9597     -112
==========================================
- Hits        54935    52395    -2540
+ Misses       6128     5741     -387
+ Partials     2676     2461     -215

Δ = absolute <relative> (impact), ø = not affected, ? = missing data

### jona-sassenhagen reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 return ci_low, ci_up # Now we can track the emergence of spatial patterns compared to baseline # for each freqyenct band

#### jona-sassenhagen Mar 31, 2017

Contributor

frequency ... also spatial patterns?

#### dengemann Mar 31, 2017

Member

what is the problem? GFP indicates variance across sensors.

#### dengemann Mar 31, 2017

Member

but thanks for pointing me to a typo :)

### jona-sassenhagen reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 ############################################################################### # Now we can compute the Global Field Power # We first estimae the rank as this data is rank-reduced as SSS was applied.

Contributor

estimate :D

Member

### dengemann commented Mar 31, 2017

 @jona-sassenhagen you want to run it on your laptop to confirm that it is relatively fast / runs on a weakly performant architecture (I developed this on a 12" MB)? Final verdict?

### jona-sassenhagen reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 fig, axes = plt.subplots(4, 1, figsize=(10, 7), sharex=True, sharey=True) colors = [plt.cm.viridis(ii) for ii in (0.1, 0.35, 0.75, 0.95)]

#### jona-sassenhagen Mar 31, 2017

Contributor

colors = plt.cm.viridis((0.1, 0.35, 0.75, 0.95))?

Member

indeed ...

### jona-sassenhagen reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 gfps_bs = np.array(gfps_bs) gfps_bs = mne.baseline.rescale(gfps_bs, average.times, baseline=(None, 0)) ci_low = np.percentile(gfps_bs, 2.5, axis=0) ci_up = np.percentile(gfps_bs, 97.5, axis=0)

#### jona-sassenhagen Mar 31, 2017

Contributor

Can't you do ci = np.percentile(gfps_bs, (2.5, 97.5), axis=0)
and then later

ax.fill_between(times, *(gfp + ci), color=color, alpha=0.3)
or something like that

#### dengemann Mar 31, 2017

Member

maybe. But the code I present here is at least more explicit. Let me at least avoid two calls to np.percentile.

Contributor

### jona-sassenhagen commented Mar 31, 2017

 I don't have Somato here and the wifi is very bad .... I ran it on our server and it took a while. It would be good to optimize it further, if at all possible.
Member

### dengemann commented Mar 31, 2017

 I ran it on our server and it took a while. It would be good to optimize it further, if at all possible. It's difficult to believe this. It takes less then one minute on my 12" macbook (remember it's gold on top of it) that does not even have an intel i5 chip but something that's closer to an iPad ...
Contributor

### jona-sassenhagen commented Mar 31, 2017

 Well 'a while' < 1 min :)
Member

### dengemann commented Mar 31, 2017

 Well 'a while' < 1 min :) Ok but this is then indeed ok. Filtering four time takes a moment. It can only get faster on more serious computers (no gold).

### jona-sassenhagen reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 gfp_bs = np.sum(average.data[bs_indices] ** 2, 0) / rank gfps_bs.append(gfp_bs) gfps_bs = np.array(gfps_bs) gfps_bs = mne.baseline.rescale(gfps_bs, average.times, baseline=(None, 0))

#### jona-sassenhagen Mar 31, 2017

Contributor
def get_gfp_ci(average, rank=rank, n_draws=2000):
"""get confidence intervals from non-parametric bootstrap"""
indices = np.arange(len(average.ch_names), dtype=int)
gfps_bs = np.zeros((n_draws, average.data.shape[1]))
for iter_ in range(n_draws):
bs_indices = rng.choice(indices, replace=True, size=len(indices))
gfps_bs[iter_, :] = (average.data[bs_indices] ** 2).sum(0) / rank
gfps_bs = mne.baseline.rescale(gfps_bs, average.times, baseline=(None, 0))
ci_low, ci_up = np.percentile(gfps_bs, (2.5, 97.5), axis=0)
return ci_low, ci_up

?

#### dengemann Mar 31, 2017

Member

let's be a bit softer here. This function in the middle of the example is already enough I feel ... making a list and then an array is a stable pattern in research code. For a module level version I'd prefer your code.

#### dengemann Mar 31, 2017

Member

Ok, I tried to take into account your proposal.

#### jona-sassenhagen Mar 31, 2017 • edited Edited 1 time jona-sassenhagen edited Mar 31, 2017 (most recent)

Contributor

I'll probably steal it to turn it into a private function anyways for plot_compare evokeds, so then I can golf it all I want :D

Member

### dengemann commented Mar 31, 2017

 @agramfort wanna have a look at this. I think we're good to go. … On Fri, Mar 31, 2017 at 12:31 AM jona-sassenhagen ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In examples/time_frequency/plot_time_frequency_global_field_power.py <#4139 (comment)>: > +print('The numerical rank is: %i' % rank) +rng = np.random.RandomState(42) + +# Then we prepare a bootstrapping function to estimate confidence intervals + + +def get_gfp_ci(average, rank=rank): + """get confidence intervals from non-parametric bootstrap""" + indices = np.arange(len(average.ch_names), dtype=int) + gfps_bs = list() + for bootstrap_iteration in range(2000): + bs_indices = rng.choice(indices, replace=True, size=len(indices)) + gfp_bs = np.sum(average.data[bs_indices] ** 2, 0) / rank + gfps_bs.append(gfp_bs) + gfps_bs = np.array(gfps_bs) + gfps_bs = mne.baseline.rescale(gfps_bs, average.times, baseline=(None, 0)) def get_gfp_ci(average, rank=rank, n_draws=2000): """get confidence intervals from non-parametric bootstrap""" indices = np.arange(len(average.ch_names), dtype=int) gfps_bs = np.zeros((n_draws, average.data.shape[1])) for iter_ in range(n_draws): bs_indices = rng.choice(indices, replace=True, size=len(indices)) gfps_bs[iter_, :] = (average.data[bs_indices] ** 2).sum(0) / rank gfps_bs = mne.baseline.rescale(gfps_bs, average.times, baseline=(None, 0)) ci_low, ci_up = np.percentile(gfps_bs, (2.5, 97.5), axis=0) return ci_low, ci_up ? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#4139 (review)>, or mute the thread .
Contributor

### jona-sassenhagen commented Mar 31, 2017

 Maybe a better title would be "Explore event-related synchronisation/desynchronisation for multiple bands" so people who do that kind of stuff can easily find it?
Member

### dengemann commented Mar 31, 2017

 I don't like that terminology as it is cultural / group specific. I'd keep it general. … On Fri, Mar 31, 2017 at 9:30 AM jona-sassenhagen ***@***.***> wrote: Maybe a better title would be "Explore event-related synchronisation/desynchronisation for multiple bands" so people who do that kind of stuff can easily find it? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#4139 (comment)>, or mute the thread .

Member

### dengemann commented Mar 31, 2017

 @jona-sassenhagen good for you? @Eric89GXL wanna have a look?

### larsoner reviewed Mar 31, 2017

Example:

https://4229-1301584-gh.circle-artifacts.com/0/home/ubuntu/mne-python/doc/_build/html/auto_examples/time_frequency/plot_time_frequency_global_field_power.html

Probably worth doing set(xlim=...) to avoid unused space to the right of each plot.

Otherwise LGTM

examples/time_frequency/plot_time_frequency_global_field_power.py
 The objective is to show you how to explore spectrally localized effects. For this purpose we adapt the method described in [1] and use it on the somato dataset. The idea is to track the band-limited temporal evolution of spatial patterns by using the Global Field Power (GFP).

#### larsoner Mar 31, 2017

Member

Need newline between paragraphs to render properly

examples/time_frequency/plot_time_frequency_global_field_power.py
 ============================================= The objective is to show you how to explore spectrally localized effects. For this purpose we adapt the method described in [1] and use it on

#### larsoner Mar 31, 2017

Member

[1]_

examples/time_frequency/plot_time_frequency_global_field_power.py
 We first bandpass filter the signals and then apply a Hilbert transform. To reveal oscillatory activity the evoked response is then subtracted from every single trial. Finally, we rectify the signals prior to averaging across trials. Then the GFP is computed as described in [2], using the

#### larsoner Mar 31, 2017

Member

[2]_

examples/time_frequency/plot_time_frequency_global_field_power.py
 # bandpass filter and compute Hilbert raw.filter(fmin, fmax, n_jobs=1) # use more jobs to speed up. raw.apply_hilbert(n_jobs=1, envelope=False)

#### larsoner Mar 31, 2017

Member

why envelope=False but epochs._data = abs later? For baseline purposes?

#### dengemann Mar 31, 2017

Member

I thought it mattered for the subtraction of the evoked response to first subtract and then rectify.

#### agramfort Mar 31, 2017

Member

to be able to subtract the evoked

examples/time_frequency/plot_time_frequency_global_field_power.py
 effects. For this purpose we adapt the method described in [1] and use it on the somato dataset. The idea is to track the band-limited temporal evolution of spatial patterns by using the Global Field Power (GFP). We first bandpass filter the signals and then apply a Hilbert transform. To

#### larsoner Mar 31, 2017

Member

... and take the magnitude

examples/time_frequency/plot_time_frequency_global_field_power.py
 frequencies. The procedure is then repeated for each frequency band of interest and all GFPs are visualized. To estimate uncertainty, non-parametric confidence intervals are computed as described in [3] across channels.

#### larsoner Mar 31, 2017

Member

[3]_

examples/time_frequency/plot_time_frequency_global_field_power.py
 raw.pick_types(meg='grad', eog=True) # we just look at gradiometers # bandpass filter and compute Hilbert raw.filter(fmin, fmax, n_jobs=1) # use more jobs to speed up.

#### larsoner Mar 31, 2017

Member

these filters will have different transition bandwidths, because l_trans_bandwidth='auto' mode makes the bandwidth depend on the freq. To make them more comparable, set all these to the same value (maybe 1?)

#### dengemann Mar 31, 2017

Member

good call, checked it out and it looks good.

examples/time_frequency/plot_time_frequency_global_field_power.py
 sharex=True, sharey=True) colors = plt.cm.viridis((0.1, 0.35, 0.75, 0.95)) for ((freq_name, fmin, fmax), average), color, ax in zip( frequency_map, colors, reversed(axes.ravel())):

#### larsoner Mar 31, 2017

Member

why not just axes[::-1, 0]?

#### dengemann Mar 31, 2017

Member

reversed is more explicit.

### agramfort reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 reject=dict(grad=4000e-13, eog=350e-6), preload=True) # remove evoked response and get analytic signal (envelope) epochs.subtract_evoked() epochs._data = np.abs(epochs.get_data())

#### agramfort Mar 31, 2017

Member

hacking a private attribute is problematic :(

### agramfort reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 for iteration in range(n_bootstraps): bs_indices = rng.choice(indices, replace=True, size=len(indices)) gfps_bs[iteration] = np.sum(average.data[bs_indices] ** 2, 0) / rank gfps_bs = np.array(gfps_bs)

Member

### agramfort reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 def get_gfp_ci(average, rank=rank, n_bootstraps=2000): """get confidence intervals from non-parametric bootstrap""" indices = np.arange(len(average.ch_names), dtype=int) gfps_bs = np.zeros((n_bootstraps, len(average.times)))

Member

np.empty

### agramfort reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 # Now we can track the emergence of spatial patterns compared to baseline # for each frequency band

#### agramfort Mar 31, 2017

Member

too many blank lines

### agramfort reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 axes.ravel()[-1].set_xlabel('Time [ms]') # We see dominant responses in the Alpha and Beta bands.

#### agramfort Mar 31, 2017

Member

I would put this in a header of the last cell

Member

Member

### agramfort reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 # Now we can compute the Global Field Power # We first estimate the rank as this data is rank-reduced as SSS was applied. # Therefore the degrees of freedom are less then the number of sensors.

#### agramfort Mar 31, 2017 • edited Edited 1 time agramfort edited Mar 31, 2017 (most recent)

Member

there is no rank anymore

### agramfort reviewed Mar 31, 2017

 single trial. Finally, we rectify the signals prior to averaging across trials by taking the magniude of the Hilbert. Then the GFP is computed as described in [2]_, using the sum of the squares but without normalization by the rank.

Member

no rank anymore

#### dengemann Mar 31, 2017

Member

yes this is what I write here :) see "without"

### agramfort reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 for ((freq_name, fmin, fmax), average), color, ax in zip( frequency_map, colors, reversed(axes.ravel())): times = average.times * 1e3 gfp = np.sum(average.data ** 2, 0)

#### agramfort Mar 31, 2017

Member

gfp = np.sum(average.data ** 2, axis=0)

### agramfort reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 fig, axes = plt.subplots(4, 1, figsize=(10, 7), sharex=True, sharey=True) colors = plt.cm.viridis((0.1, 0.35, 0.75, 0.95)) for ((freq_name, fmin, fmax), average), color, ax in zip( frequency_map, colors, reversed(axes.ravel())):

#### agramfort Mar 31, 2017

Member

axes.ravel()[::-1]

would keep it an array

### agramfort reviewed Mar 31, 2017

examples/time_frequency/plot_time_frequency_global_field_power.py
 ############################################################################### # Now we can compute the Global Field Power

#### agramfort Mar 31, 2017

Member

I would make it a full block of text and don't have the inline comments that won't pop out

Contributor

### jona-sassenhagen commented Mar 31, 2017

 Since the gfp_bs = np.array(gfp_bs) line is gone, I'm happy :) Would prefer to have a title that advertises this to ERD/ERS people though.
Member

### dengemann commented Mar 31, 2017

 Are you happy if we say it in the text? By putting it in the title we take position. This is not neutral. … On Fri, Mar 31, 2017 at 5:25 PM jona-sassenhagen ***@***.***> wrote: Since the gfp_bs = np.array(gfp_bs) line is gone, I'm happy :) Would prefer to have a title that advertises this to ERD/ERS people though. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#4139 (comment)>, or mute the thread .
Contributor

### jona-sassenhagen commented Mar 31, 2017

 Well what'd be most user-friendly? How do you actually reach people? I think ERD/ERS is pretty common, and that's what people who need this code are probably looking for. (Not specifically these words - but something that's more descriptive of what's actually going on wrt. analysing experimental data, e.g. 'Event-Related Spectral Perturbations' as Makeig et al. would say.
Contributor

### jona-sassenhagen commented Mar 31, 2017

 I mean, if you want to reach different people, that's a different thing.
Member

### dengemann commented Mar 31, 2017

 I mean, if you want to reach different people, that's a different thing. I would veto this. ERD/ERS is too local. Only certain people call it like this and this goes to far.
Contributor

### jona-sassenhagen commented Mar 31, 2017

 I don't mean those specific words - just something that provides the required information. Explore oscillatory activity in sensor space doesn't really say much about what's specific about this - that it's 1. a smart way of getting specific bands, and 2. about their event-related dynamics.
Member

### dengemann commented Mar 31, 2017

 @jona-sassenhagen propse a name that is informative but does not explicitly buy into certain theoretical camps
Contributor

### jona-sassenhagen commented Mar 31, 2017

 How about "Explore event-related dynamics for specific frequency bands [in sensor space]"? Or "Unbiased isolation of frequency band specific time courses" ... something like this. I mean, the current title isn't bad, I just think it could sell the thing to more people :)
Member

### dengemann commented Mar 31, 2017

 I mean, the current title isn't bad, I just think it could sell the thing to more people :) ok I'm totally for this :) "Explore event-related dynamics for specific frequency bands ..." sounds good.
Contributor

### jona-sassenhagen commented Mar 31, 2017

 I actually really prefer this approach to the currently much more popular broad-band ERSP approach. As Lopes Da Silva and Pfurtscheller write: "the term ERD is only meaningful if the baseline measured some seconds before the event represents a rhythmicity seen as a clear peak in the power spectrum. Similarly, the term ERS only has a meaning if the event results in the appearance of a rhythmic component and therewith in a spectral peak that was initially not detect- able (Lopes da Silva and Pfurtscheller, 1999)." I guess people just love 2D arrays in Jet too much though.
Member

### dengemann commented Apr 1, 2017

 Does anyone wants to have another look before merging? It seems comments by 3 reviewers are addressed.
Member

### agramfort commented Apr 1, 2017

 refs don't see to render fine in circle: https://4269-1301584-gh.circle-artifacts.com/0/home/ubuntu/mne-python/doc/_build/html/auto_examples/time_frequency/plot_time_frequency_global_field_power.html#id4
Member

### agramfort merged commit b6bb7bf into mne-tools:master Apr 1, 2017 2 of 3 checks passed

#### 2 of 3 checks passed

continuous-integration/travis-ci/pr The Travis CI build is in progress
Details
ci/circleci Your tests passed on CircleCI!
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details