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

news/2024/11/23 5:37:02/

使用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