An example showing the basic usage of wICA algorithm

This code is for illustration of the method described in: N.P. Castellanos, and V.A. Makarov (2006). "Recovering EEG brain signals: Artifact suppression with wavelet enhanced independent component analysis" J. Neurosci. Methods, 158, 300–312.

Requirements: runica from EEGLAB toolbox, and rwt - Rice Wavelet Toolbox (both freely available in internet).

This code is copyright © by the authors, and we hope you acknowledge our work. We distribute it in the hope that it will be useful, but without any warranty.

Author: Valeri A. Makarov e-mail: vmakarov@opt.ucm.es

2006

Contents

Reset environment

close all
clear all

Read test data

Fs = 256;  % Sampling frequency in Hz
[Data, ChanTitels, FileName] = ReadTruScan_ascii('./Data/test1.tdt');

Conventional notch filter (you may add more)

This is an optional step (suppress 50Hz noise)

Fnyq  = Fs/2;
F_notch = 50; % Notch at 50 Hz
[b,a] = iirnotch(F_notch/Fnyq, F_notch/Fnyq/20);
Data  = filtfilt(b,a, Data);

Conventional High Pass Filter

This is an optional step (suppress low <4Hz frequency noise)

F_cut = 4;
[b,a] = ellip(1, 0.5, 20, F_cut/Fnyq, 'high');
Data  = filtfilt(b,a, Data);

Remove mean values from the channels and plot raw data

Data  = detrend(Data,'constant');
% Transpose the data matrix to get (channel x time) orientation
Data = Data';
figure;
PlotEEG(Data, Fs, ChanTitels, 200, 'Raw data (19 channels 10-20 system)');

Find independent components

EEGLAB is required!!! You can also use other algorithms, e.g. fICA. Note, the use of long (in time) data sets reduces the algorithm performance see for details the abovementioned paper.

[weight, sphere] = runica(Data, 'verbose', 'off');
W = weight*sphere;    % EEGLAB --> W unmixing matrix
icaEEG = W*Data;      % EEGLAB --> U = W.X activations
figure;
PlotEEG(icaEEG, Fs, [], [], 'Independent Components');
xlabel('Time (s)')

Inverting negative activations: 1 2 -3 4 5 6 -7 8 9 -10 11 12 -13 -14 15 -16 17 18 -19 

Common ICA artifact suppression method

Remove artifact components and rebuild signals

An input dialog will be open, where you should enter the ICs to be removed (those with artifacts)

answer  = inputdlg({'Artifact ICs'},'Components to remove',1);
ArtICs = str2num(answer{1});
disp(['Suppressed components: ' answer{1}]);
icaEEG1 = icaEEG;
icaEEG1(ArtICs, :) = 0;          % suppress artifacts
Data_ICA = inv(W)*icaEEG1;       % rebuild data
figure;
PlotEEG(Data_ICA, Fs, ChanTitels, 100, 'ICA cleanned EEG');

Suppressed components: 1

wICA artifact rejection

Rice Wavelet Toolbox (rwt) is required!!!

You can also tune the third argument, or/and provide (in the second argument) numbers of only selected components to be processed by wavelet e.g. [1:4,7] (in the example we use all 19 components).

[icaEEG2, opt]= RemoveStrongArtifacts(icaEEG, (1:19), 1.25, Fs);
Data_wICA = inv(W)*icaEEG2;
figure;
PlotEEG(icaEEG2, Fs, [], [], 'wavelet filtered ICs');
figure;
PlotEEG(Data_wICA, Fs, ChanTitels, 100, 'wICA cleanned EEG');

The component #1 has been filtered
The component #2 has been filtered
The component #3 has been filtered
The component #4 has been filtered
The component #5 has passed unchaneged
The component #6 has been filtered
The component #7 has passed unchaneged
The component #8 has passed unchaneged
The component #9 has been filtered
The component #10 has been filtered
The component #11 has been filtered
The component #12 has been filtered
The component #13 has passed unchaneged
The component #14 has been filtered
The component #15 has passed unchaneged
The component #16 has been filtered
The component #17 has been filtered
The component #18 has been filtered
The component #19 has been filtered

Make a comparative plot

We plot a segment with artifacts and a segment outside of (ocular) artifacts. Note that outside of artifacts an artifact suppression method should conserve (not perturb) the original EEG signal.

figure
subplot(2,1,1);
T1 = 2; T2 = 6;
nT1 = T1*Fs; nT2 = T2*Fs;
plot((nT1:nT2)/Fs, Data(1,nT1:nT2),'k',(nT1:nT2)/Fs, Data_ICA(1,nT1:nT2),'b',...
    (nT1:nT2)/Fs, Data_wICA(1,nT1:nT2),'r');
legend({'Raw EEG','ICA cleaned EEG','wICA cleaned EEG'},'Location','NorthWest');
title('Segment with artifacts (FP1)')
subplot(2,1,2);
T1 = 2.5; T2 = 3;
nT1 = T1*Fs; nT2 = T2*Fs;
plot((nT1:nT2)/Fs, Data(1,nT1:nT2),'k',(nT1:nT2)/Fs, Data_ICA(1,nT1:nT2),'b',...
    (nT1:nT2)/Fs, Data_wICA(1,nT1:nT2),'r');
title('Amplified, outside of artifacts')
xlabel('Time (s)')
figure
T1 = 8.2; T2 = 9.2;
nT1 = T1*Fs; nT2 = T2*Fs;
plot((nT1:nT2)/Fs, Data(1,nT1:nT2),'k',(nT1:nT2)/Fs, Data_ICA(1,nT1:nT2),'b',...
    (nT1:nT2)/Fs, Data_wICA(1,nT1:nT2),'r');
legend({'Raw EEG','ICA cleaned EEG','wICA cleaned EEG'},'Location','NorthWest');
title('Segment with another artifact')
xlabel('Time (s)')

Warning: Integer operands are required for colon operator when used as index
Warning: Integer operands are required for colon operator when used as index
Warning: Integer operands are required for colon operator when used as index

Published with MATLAB® 7.2