使用ICA去除32通道脑电数据中的眼电

% 去除听音乐1的脑电数据眼电

% Method: ICA

% 1.导入数据 32通道 1000hz采样率

% 2. 数据预处理(50hz陷波+0.3hz-250hz带通，可以自己改参数调整)

% 3. ICA(下载的FastICA_25)

% 4. 计算ICA的各ICs分量与原水平、垂直眼电信号相关性

% 5. 将眼电噪声的ICs置零后重建成没有眼电的脑电信号

% 注意代码中设置断点

% copyright：seu-ww

% 2019-05-07

clc;clear all;close all;

MUSIC1PATH=‘D:\Share\music_breath_2018\1-EEG\MUSIC2’;

SAVEPATH1 = ‘D:\Share\music_breath_2018\1-EEG\CLEAN_DATA\MUSIC2’;

data_list = dir(MUSIC1PATH);

eeg_number = length(data_list)-2;

fs = 1000;

for i =3:48 %45名被试的数据

tic;

filename = data_list(i).name;

datapath = [MUSIC1PATH,filename];

num_med = regexp(data_list(i).name,‘D’,‘split’);

subjectID = str2num(cell2mat(num_med(1)));

load(datapath);

data_org = data_music2([1:18,20:23,25:36]?;% M1、M2

%% Preprocessing(Notch + bandpass)

Wo = 50/(fs/2); BW = Wo/35;

[b,a] = iirnotch(Wo,BW);

data_notch= filtfilt(b,a,(double(data_org(:?))’)’;

[b,a] = butter(3,[0.3/(fs/2) 250/(fs/2)],‘bandpass’);

data_bp= filtfilt(b,a,data_notch’)’;%使用零相移数字滤波器(IIR)

L=length(data_bp);t=(0:(L-1))/fs;

%% ICA

level = 32; %ICA分层数，最高设为脑电的通道数

[Zica,A,W] = fastica(data_bp,‘numOfIC’, level);

plot_ICAs(Zica,fs,10);

%% Correlation with VEOG and HEOG

for k=1:(length(Zica(:,1)))

s=corrcoef(data_bp(33,:)’,Zica(k,:)’);

ss1(k) = s(2,1);

s=corrcoef(data_bp(34,:)’,Zica(k,:)’);

ss2(k) = s(2,1);

end

figure;plot(abs(ss1));title(‘各ICs与原水平眼电信号相关性’);

figure;plot(abs(ss2));title(‘各ICs与原垂直眼电信号相关性’);

%% Reconstruction

rec_no = [3，7] %这里输入的是与原垂直眼电信号相关性最高的IC分量

Zica_del = Zica; Zica_eye = zeros(size(Zica,1),size(Zica,2));

Zica_del(rec_no,:) = 0; Zica_eye(rec_no,:) = Zica(rec_no,:);

data_rec = AZica;

data_clean = AZica_del;

data_eog = A*Zica_eye;

%去眼电前后对比(以FP1通道的数据为例)

figure;plot(t,data_bp(1,:),‘k’);title('32通道ICA去眼电前EEG信号 ');

figure;plot(t,data_clean(1,:),‘k’);title('32通道ICA去眼电后EEG信号 ');

data_music2=data_clean;

save([SAVEPATH1,filename],‘data_music2’); %保存数据

end

%画ICA分解后ICs分量的函数plot_ICAs

% 画出ICA 独立分量

% Zica : ICs

% fs : fs

% d : offset

function plot_ICAs(Zica,fs,d)

level = length(Zica(:,1));

L = length(Zica(1,:));

t=(0:(L-1))/fs;

figure;

offs= 1:1:level;

d= 10;

for ics = 1:level

plot(t,Zica(ics,:) + offs(ics)d,‘k’);hold on;

end

ylim([0 d(level+1)]); title(‘Independent Components’);

IC = regexp(sprintf('IC-%02d ‘,[1:level]),’ ',‘split’);

set(gca,‘YTick’,[d:d:d*(level)],‘YTickLabel’,IC(1:end));

end

ICA分解后的独立分量

3明显为垂直眼电，7明显为水平眼电。

去眼电前后的脑电信号对比：