[Mne_analysis] Hilbert Transform for Time Frequency Analysis (Aditya P Singh)

Aditya P Singh adityasingh at utexas.edu
Fri Jul 7 03:21:02 EDT 2017
Search archives:

Hi all,

Thank you Philip for the help, and I used the loop example but I am not
exactly sure how to save the amplitude values for each frequency and time
window in the code. I have attached some of my code as a .py file but I
think just looking at the example:

tf_scores[freq, t] = np.mean(cross_val_score(estimator=clf, X=X, y=y,
                                             scoring='roc_auc', cv=cv,
n_jobs=1), axis=0)

And it is initialized as:

tf_scores = np.zeros((n_freqs - 1, n_windows))

which I changed the tf_scores to the amplitude by applying the hilbert
transform then drawing out the amplitude

 h = raw_filter.apply_hilbert(picks,envelope = True)
 a = raw_filter.apply_function(np.abs, picks)

So instead of tf_scores, it would be

 amplitude[freq, t] = [a._data[picks[0]]]

But here is where I am having the trouble since I am not not sure how
to save the sequence of amplitudes over time for each frequency band
in the loop? I am sure it is a simple fix, but I am just not sure how
to do it. Also what should I be setting my vmin as in the tfr plot? I
have attached my code for further detail.


Thanks,

Aditya
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.nmr.mgh.harvard.edu/pipermail/mne_analysis/attachments/20170707/9d6a045a/attachment-0001.html 
-------------- next part --------------
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Fri Jul  7 00:33:58 2017

@author: aditya
"""
import numpy as np
import os.path as op
import mne
import matplotlib.pyplot as plt
from mne import Epochs, find_events, create_info
from mne.time_frequency import AverageTFR

raw = mne.io.read_raw_fif(
        "/Users/aditya/desktop/MEG_TrainingData/SelfPaced_ButtonPress/mm_selfpaced_index_rt_raw.fif",
        preload = True)
event_id, tmin, tmax = 1, -2, 2 #specify times and event id 
events = mne.find_events(raw, stim_channel='STI101')


# Extract information from the raw file
sfreq = raw.info['sfreq']

min_freq = 5.
max_freq = 30.
n_freqs = 10  # how many frequency bins to use

# Assemble list of frequency range tuples
freqs = np.linspace(min_freq, max_freq, n_freqs)  # assemble frequencies
freq_ranges = list(zip(freqs[:-1], freqs[1:]))  # make freqs list of tuples

# Infer window spacing from the max freq and number of cycles to avoid gaps
window_spacing = (n_cycles / np.max(freqs) / 2.)
centered_w_times = np.arange(tmin, tmax, window_spacing)[1:]
n_windows = len(centered_w_times)


# init 
r = np.zeros((n_freqs - 1, n_windows))

# Loop through each frequency range of interest
for freq, (fmin, fmax) in enumerate(freq_ranges):

    # Infer window size based on the frequency being used
    w_size = n_cycles / ((fmax + fmin) / 2.)  # in seconds

    # Apply band-pass filter to isolate the specified frequencies
    raw_filter = raw.copy().filter(fmin, fmax)
    picks = mne.pick_types(raw_filter.info, meg= True, eeg=True, stim=True, eog=False,
                       exclude='bads')
    h = raw_filter.apply_hilbert(picks,envelope = True)
    a = raw_filter.apply_function(np.abs, picks)

#     Extract epochs from filtered data, padded by window size
#    epochs = Epochs(raw_filter, events, event_id, tmin - w_size, tmax + w_size,
#                    proj=False, baseline=None, preload=True)
#    epochs.drop_bad()
#    y = le.fit_transform(epochs.events[:, 2])

    # Roll over time
    for t, w_time in enumerate(centered_w_times):

        # Center the min and max of the window
        w_tmin = w_time - w_size / 2.
        w_tmax = w_time + w_size / 2.

        # Crop data into time-window of interest
#        X = epochs.copy().crop(w_tmin, w_tmax).get_data()

        # Save mean scores over folds for each frequency and time window
        amplitude[freq, t] = [a._data[picks[0]],t]
        
        
# Set up time frequency object
av_tfr = AverageTFR(create_info(['freq'], sfreq), amplitude[np.newaxis, :],
                    centered_w_times, freqs[1:], 1)

chance = np.mean(y)  # set chance level to white in the plot
av_tfr.plot([0], vmin=0, title="Time-Frequency Decoding Scores",
            cmap=plt.cm.Reds)


More information about the Mne_analysis mailing list