Failed Attempt to use Time Shift PCA

test_tspca.py

# Author: Jean-Remi King
#
# Licence: BSD 3-clause
"""
This is a failed attempt to use a python implementation of Time Shift
PCA (https://github.com/pealco/python-meg-denoise) that aims as
denoising external MEG sources from the signals when we have reference
sensors available (e.g. in the KIT).

See http://audition.ens.fr/adc/NoiseTools/ for more info on TSPCA

The scripts fails for both Raw and Epochs data, due to some shape
issues.
"""

import os
import numpy as np

import mne
from mne import Epochs
from mne.io import read_raw_kit

from denoise import demean, fold, unfold
from tspca import tsr
from sns import sns
from dss import dss1


data_path = '/media/DATA/Pro/Projects/NewYork/audio_wm/data/'
subject_dir = 'R1022_Audio_Sequence_11.25.15'
fname_2tones = 'R1022_2Tones_11.25.15.sqd'
bad_sensors = [20, 40, 64, 112, 115, 152]

fname = os.path.join(data_path, subject_dir, fname_2tones)
raw = read_raw_kit(fname, preload=True)
events = mne.find_events(raw)
for bad in bad_sensors:
    raw.info['bads'] += ['MEG %.3i' % bad]

ch_types = {6001: 'grad', 6002: 'mag', 0: 'misc'}
ch_type = np.array([ch_types[ii['coil_type']] for ii in raw.info['chs']])
data_raw = raw._data.transpose()[:10000, :, None]  # time x chan x trial
epochs = Epochs(raw, events, event_id=[1, 2], tmin=0, tmax=.400,
                baseline=None, preload=True, verbose=False)
epochs.decimate(10)
data_epochs = epochs._data.transpose([2, 1, 0])

for data in [data_epochs, data_raw]:
    # data must be: time x chan x trial
    # data = np.concatenate([data, data], axis=2)
    meg_head = data[:, np.where(ch_type == 'grad')[0], :]
    meg_ref = data[:, np.where(ch_type == 'mag')[0], :]

    # remove means
    noisy_data = demean(meg_head)[0]
    noisy_ref = demean(meg_ref)[0]

    # shifts = np.arange(-50, 51) ???
    data_tspca, idx = tsr(noisy_data, noisy_ref)[0:2]
    data = data[idx, :, :]
    data_mean = np.mean(data, 2)
    data_tspca_mean = np.mean(data_tspca, 2)

    # apply SNS
    nneighbors = 10
    data_tspca_sns = sns(data_tspca, nneighbors)

    # apply DSS
    # --- Keep all PC components
    data_tspca_sns = demean(data_tspca_sns)[0]
    todss, fromdss, ratio, pwr = dss1(data_tspca_sns)

    # c3 = DSS components
    data_tspca_sns_dss = fold(np.dot(unfold(data_tspca_sns), todss),
                              data_tspca_sns.shape[0])