In [1]:
import numpy as np
import sys

# comment out or adjust as appropriate:
sys.path.insert(0, '/home/ctw/src/ptsa_new/build/lib.linux-x86_64-2.7/')
import ptsa

from ptsa.data.filters.MorletWaveletFilterCpp import MorletWaveletFilterCpp
from ptsa.data.readers import EEGReader #EDFRawReader
from ptsa.data.filters import MonopolarToBipolarMapper
from ptsa.data.readers import JsonIndexReader
from ptsa.data.readers import BaseEventReader


# when running on rhino set rhino_mount to '', otherwise adjust as appropriate:
rhino_mount = '/home/ctw/fusemounts/rhino'

In [2]:
ev = BaseEventReader(filename=rhino_mount+
                     '/data/eeg/scalp/ltp/ltpFR2/behavioral/events/events_all_LTP093.mat',
                     use_reref_eeg=True,common_root='data').read()
ev = ev[ev['type'] == 'WORD']
ev = ev[ev['session'] == 1]

channels = np.array(['{:03}'.format(x) for x in range(1,130)])
eeg = EEGReader(events=ev, channels=channels, start_time=0.0, end_time=1.0, buffer_time=2.0).read()

freqs=np.logspace(np.log10(2), np.log10(200), 15)

In [3]:
wf = MorletWaveletFilterCpp(time_series=eeg,
                            freqs=freqs,
                            output='power',
                            frequency_dim_pos=0,
                            verbose=True,
                            cpus=1)
%timeit pow_wavelet, phase_wavelet = wf.filter()


('CPP total time wavelet loop: ', 102.09135293960571)
('CPP total time wavelet loop: ', 106.2283661365509)
('CPP total time wavelet loop: ', 105.33881711959839)
('CPP total time wavelet loop: ', 105.72558093070984)
1 loop, best of 3: 1min 46s per loop

In [4]:
wf = MorletWaveletFilterCpp(time_series=eeg,
                            freqs=freqs,
                            output='power',
                            frequency_dim_pos=0,
                            verbose=True,
                            cpus=10)
%timeit pow_wavelet, phase_wavelet = wf.filter()


('CPP total time wavelet loop: ', 11.442231893539429)
('CPP total time wavelet loop: ', 11.50298810005188)
('CPP total time wavelet loop: ', 11.33631181716919)
('CPP total time wavelet loop: ', 11.525317192077637)
1 loop, best of 3: 12.7 s per loop

In [5]:
import mne
info = mne.create_info(129, 500)
eparray = np.transpose(np.array(eeg), (1,0,2))
dat_mne = mne.EpochsArray(eparray, info)


576 matching events found
0 projection items activated
0 bad epochs dropped

In [6]:
%timeit power = mne.time_frequency.tfr_array_morlet(epoch_data=dat_mne, sfreq=500, freqs=freqs, n_cycles=5, output='power', n_jobs=1)


1 loop, best of 3: 5min 12s per loop

In [7]:
%timeit power = mne.time_frequency.tfr_array_morlet(epoch_data=dat_mne, sfreq=500, freqs=freqs, n_cycles=5, output='power', n_jobs=10)


[Parallel(n_jobs=10)]: Done  52 tasks      | elapsed:   51.5s
[Parallel(n_jobs=10)]: Done 129 out of 129 | elapsed:  2.1min finished
[Parallel(n_jobs=10)]: Done  52 tasks      | elapsed:   52.1s
[Parallel(n_jobs=10)]: Done 129 out of 129 | elapsed:  2.1min finished
[Parallel(n_jobs=10)]: Done  52 tasks      | elapsed:   51.7s
[Parallel(n_jobs=10)]: Done 129 out of 129 | elapsed:  2.1min finished
[Parallel(n_jobs=10)]: Done  52 tasks      | elapsed:   51.9s
[Parallel(n_jobs=10)]: Done 129 out of 129 | elapsed:  2.1min finished
1 loop, best of 3: 2min 24s per loop

In [ ]: