## Testing Connectivity

Here we explore how different paddings of the MEG recordings affect the phase after applying the Hilbert transform.

Start Matlab session:

``````

In [1]:

``````
``````

Starting MATLAB on ZMQ socket ipc:///tmp/pymatbridge-48c7880e-8c96-4e6f-933b-7f755877a622
Send 'exit' command to kill the server
......MATLAB started and connected!

/usr/local/lib/python2.7/dist-packages/IPython/nbformat.py:13: ShimWarning: The `IPython.nbformat` package has been deprecated. You should import from nbformat instead.
"You should import from nbformat instead.", ShimWarning)

``````
``````

In [2]:

%%matlab

ft_defaults

``````

Plot something with Matlab:

``````

In [3]:

%%matlab
%%%%%%%%%%%%
%%% Plot Test to check that Matlab is loaded properly
%%%%%%%%%%%%

a = linspace(0.01,6*pi,100);
plot(sin(a))
grid on
hold on
plot(cos(a),'r')

``````
``````

``````

## First Part (Hilbert for a random signal)

#### Hilbert Transform

Here we look at how padding affects the phase when taking the Hilbert transform. This code is from Javier for a simple signal.

``````

In [7]:

%%matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Checks if padding a signal makes a difference when taking the
%%% Hilbert transform to find the phase.
%%%
%%% code from Javier (email, 3/07/2014) %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

signal = randn(1000,1);

hsig = hilbert(signal);

%% padding with the first and last values; ideally we would like to keep the trend too,
%% but this is just a quick example

%% the code above does something like this: for "signal" --> "ssssssssignallllllllll"

%plot(angle(signal)) % meaningless because this is not the instantaneous phase
hold on
title('Phase with angle() for simple and padded signal')
plot(angle(hsig),'b') %% this is the phase :)
plot(angle(hsigpad(101:end-100)),'r') %% the phase of the central part of the padded signal

%%%  Get analytical phase (as per http://www.scholarpedia.org/article/Hilbert_transform_for_brain_waves)

%% plot phase from atan for the initial signal and the padded one
%figure; hold on;
phase = atan(imag(hsig)./signal);
%plot(unwrap(phase), 'b')

figure; hold on;
title('Analytical phase (atan) for simple and padded signal');
plot(phase, 'b')

figure; hold on;
title('Unwrapped analytical phase atan2 for simple and padded signal');
plot(unwrap(atan2(imag(hsig), signal)), 'b')
%% these do the same as the two lines above
%plot(unwrap(angle(hsig)), 'b')

``````
``````

``````

Apparently for a complex number z:

angle(z) - is equal to theta = atan2(imag(z),real(z)). Quick check below and comparison to unwrapped phase:

``````

In [8]:

%%matlab
signal = randn(10,1);
hilbert(signal);
%phase = atan(signal./imag(hilbert(signal)))
phase = atan2(imag(hilbert(signal)), signal);
x = imag(hilbert(signal));
hold on

% the first 2 below are exactly the same, so they are indeed the same
plot(angle(hilbert(signal)), 'g')
plot(phase, 'b')
plot(unwrap(phase), 'r')

``````
``````

``````

Looking at how padding affects the phase of cleaned MEG signals when taking the Hilbert transform:

``````

In [9]:

%%matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Padds the signal with 200 columns at the beginning with the first column of the epoch
%%%     and 200 columns at the end with the last column of the epoch.
%%%
%%% OBS: In this cell, the code takes the Hilbert transform of
%%%   the "channels (rows) X time samples (columns)" matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% add data folder to path

signal = epoch.meg_no_ecg;
% size(signal) %% this is 148 (channels) x 1695 (samples)

%%% subtract mean before doing anything else?

%% plot the MEG signals
%plot(signal)

hsignal = hilbert(signal);

% size(signalPadded) %% this is 148 x 1895 --> correct as we padded with 200 columns in total
%signalPadded(:,98:102) %% extra check --> OK
%signalPadded(:, 1793:1799) %% extra check --> OK

``````
``````

In [10]:

%%matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Plots the phase for the unpadded MEG signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure; hold on;
plot(angle(hsignal),'b')

``````
``````

``````
``````

In [11]:

%%matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Plots the phase for the padded MEG signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure; hold on;

``````
``````

``````
``````

In [12]:

%%matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure; hold on;
plot(unwrap(angle(hsignal)),'b')
%% these two lines below generate the same thing w.r.t the plot above??

``````
``````

``````

Here we take the Hilbert transform of the transposed MEG matrix which yields a time samples x channels matrix. In Matlab, the hilbert() function operates columnwise. I guess this is the correct way to compute the phase.

• appending 200 columns identical to the first column of the epoch at the beginning of the signal
• appending 200 columns identical to the last column of the epoch at the end of the signal
``````

In [13]:

%%matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Load MEG signals for epochs specified in the "data" array.
%%% Plots the previous graphs.
%%%
%%% OBS: Code takes Hilbert transform of
%%%  the "time samples (rows) X channels (columns)" matrix (using transpose()).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% add data folders to path

data = {'MEGnoECG_07_segm9.mat', ...
'MEGnoECG_07_segm20.mat', ...
'MEGnoECG_AL41D_segm6.mat', ...
'MEGnoECG_AL41D_segm29'};

for i=1:length(data)
signal = epoch.meg_no_ecg;
hsignal = hilbert( transpose(signal) );

%% pad signal (200 columns at beginning and end) and take Hilbert

figure; hold on;
title(strcat('Phase for unpadded signal ', data{i}))
plot(angle(hsignal), 'b')

figure; hold on;
plot(angle(hsignalPadded(201:end-200,:)),'r') %% need to get rid of rows as I've transposed

figure; hold on;
plot(unwrap(angle(hsignal)),'b')

figure; hold on;
plot(angle(hsignal(:, end-5:end)),'b')

figure; hold on;
title(strcat('Unwrapped Phase for unpadded and padded signal (only last 5 channels)', data{i}))
plot(unwrap(angle(hsignal(:, end-5:end))),'b')
end

``````
``````

``````

## Fourth Part

Here we also compare the phases of unpadded and padded signal as in Part Three, with the difference that the padding is done using symmetric extension (also called reflective padding in the past). Padding is done using wextend().

``````

In [14]:

%%matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Load MEG signals for epochs specified in the "data" array.
%%%     Padding is done using symmetric extension.
%%% Plots the previous graphs.
%%%
%%% OBS: Code takes Hilbert transform of
%%%  the "time samples (rows) X channels (columns)" matrix (using transpose()).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% add data folders to path

data = {'MEGnoECG_07_segm9.mat', ...
'MEGnoECG_07_segm20.mat', ...
'MEGnoECG_AL41D_segm6.mat', ...
'MEGnoECG_AL41D_segm29'};

for i=1:length(data)
signal = epoch.meg_no_ecg;
hsignal = hilbert( transpose(signal) );

%% pad signal symmetrically (200 columns at beginning and end) and take Hilbert
L = [0, 200]; %% add 0 rows, 200 columns (on each side)
%% wextend(TYPE,MODE,X,L,LOC)

figure; hold on;
title(strcat('Phase for unpadded signal ', data{i}))
plot(angle(hsignal), 'b')

figure; hold on;
plot(angle(hsignalPadded(201:end-200,:)),'r') %% need to get rid of rows as I've transposed

figure; hold on;
plot(unwrap(angle(hsignal)),'b')

figure; hold on;
plot(angle(hsignal(:, end-5:end)),'b')

figure; hold on;
title(strcat('Unwrapped Phase for unpadded and padded signal (only last 5 channels)', data{i}))
plot(unwrap(angle(hsignal(:, end-5:end))),'b')
end

``````
``````

``````

This section looks at estimating instantaneous phase with FieldTrip.

We filter the our data and then apply the Hilbert transform. This can be done in _ftpreprocessing() like so:

``````

In [18]:

%%matlab

filterOrder = 750;
Fs = 169.54;

%% add data folders to path

megFiles = {'MEGnoECG_07_segm9.mat', ...
'MEGnoECG_07_segm20.mat', ...
'MEGnoECG_AL41D_segm6.mat', ...
'MEGnoECG_AL41D_segm29'};

for i=1:length(megFiles)

% configuration structure
cfg = [];
cfg.bpfilter = 'yes'; %% bandpass filter
cfg.bpfreq = [0.5 4]; %% bandpass frequency range as [low high] in Hz (delta here)
cfg.bpfiltord = filterOrder; %% bandpass filter order
cfg.bpfilttype = 'fir'; %% band pass filter type - FIR
cfg.bpfiltdir = 'twopass'; %% filter direction - two pass, like filtfilt
cfg.demean = 'yes'; %% apply baseline correction (remove DC offset)

cfg.hilbert = 'angle'; % this gives you just the phase, you can
%specify 'complex' to get both phase and amplitude

% data structure
data = [];
data.fsample = Fs;
data.trial = {epoch.meg_no_ecg};

data.time = {timeFor10s'};
startsample = 1;
endsample = size(epoch.meg_no_ecg, 2);
data.sampleinfo = [startsample endsample];

[processedData] = ft_preprocessing(cfg, data);

figure;
plot(1:1695,unwrap(processedData.trial{1}));

end

``````
``````

preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and required the additional allocation of an estimated 3 MB
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and required the additional allocation of an estimated 0 MB
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and required the additional allocation of an estimated 0 MB
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and required the additional allocation of an estimated 0 MB
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and required the additional allocation of an estimated 0 MB
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and required the additional allocation of an estimated 0 MB
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and required the additional allocation of an estimated 0 MB
preprocessing
preprocessing trial 1 from 1
the call to "ft_preprocessing" took 0 seconds and required the additional allocation of an estimated 0 MB

``````

Looking at how padding affects the phase of RAW MEG signals when taking the Hilbert transform. Have to merge raw data first...

``````

In [19]:

%%matlab

%% add RAW data folder to path

``````
``````

In [21]:

``````
``````

MATLAB closed

``````