文章目录
未经允许,禁止转载。
摘要
本问详细讲解了ongoing脑电特征与音乐特征相关分析的整个数据处理流程。对于初入MIR,自然情景范式领域的人来说,十分有必要了解一下。流程包括了6个部分:EEG预处理:带通滤波;ICA消除脑电干扰成分:去除眼电、肌电等噪声;音乐特征提取;EEG张量构建和分解;脑电成分与音乐特征的相关分析;可视化分析:脑地形图。本文就其中内容和细节展开了详细讨论,希望给各位读者一定启发和帮助。
本人也是MIR领域的一个新人,在经历一个月的实验采集和导师的指导后,感觉在数据处理和分析方面收获颇多,便写下了此文,仅供大家参考,如有偏颇之处,翘首以待各位佬的指正。
一、介绍
1. 我介绍了哪些东西?
根据实验范式的不同,EEG信号可以分为自发脑活动EEG,事件相关电位EEG(ERP)和ongoing EEG信号,分别对应于无刺激,快速重复短刺激,自然情景刺激三个场景。 本文针对的信号则是对应于自然情景刺激下的ongoing EEG。本文主要介绍相关分析技术的6个流程,并就其中细节展开深入讨论。
1)EEG预处理:带通滤波
它包括了提取指定频带的信号,滤波后噪声和被滤波频带的影响。
2)ICA消除脑电干扰成分:去除眼电、肌电等噪声
它包括了ICA简单介绍和如何使用工具箱消除脑电伪迹,该工具箱名为“ASAP_ICA_EEG_ArtifactRejection”,可以有效识别出眼电,肌电等伪迹成分,并指导消除。
3)音乐特征提取
它包括了如何使用工具箱提取音乐特征,包含了5个方面:dynamics(动态),rhythm(音律),timbre(音色),pitch(音高)和tonality(音调),每个特征方面包含了若干指标。该工具箱名为“MIRToolbox”。
5个特征术语可见该工具箱的手册,中文翻译可能略有出入。
4)EEG张量构建和张量分解
它包括了使用短时傅里叶变换将时域信号转为频域信号,得到一个张量,形状为[channel, frequency, frame],以及张量分解方法介绍与工具箱使用。同时,介绍了帧数的计算方法,即frame值。该工具箱名为“Ongoing_EEG_Tensor_Decomposition-main”。
5)脑电成分与音乐特征的相关分析
它包括了如何使用相关函数计算两段序列的相关系数。
6)可视化分析:脑地形图
它包括了脑地形的含义,如何使用MATLAB绘制脑地形图,如何查看脑地形图的响应脑区。
2. 你应该学会哪些?
1)基于内容,你应该学会独立完成图1的所有流程。
图1的流程图有两个路线,上面的路线是脑电张量构建和张量分解的流程,下面的路线是音乐特征提取的流程。两图路线合并后即使相关分析。其中,张量分解(此处使用CP算法)后得到了若干个[空间,频域,时间]的张量,这些其实是张量分解后的脑电成分。
图1。张量分解技术在脑电与音乐特征相关分析的应用。
2)基于流程,你应该学会套用该流程去执行或完善自己的实验。
这一点也是我最推崇的。理解流程的深刻内涵,并能够举一反三,才是我们学习的终极目标。
二、方法
1.EEG预处理:带通滤波
1)在线滤波和离线滤波
在线滤波和离线滤波,即online filter和offline filter。当我们阅读采集脑电信号相关的文献时,它们应该是经过会提及的术语。其本质都是滤波,但区别在于前者是采集时滤波,后者是采集后滤波。
Steven J. Luck是脑电领域顶级的教授,他在他的书《Introduction to Event-Related Potential Technique》确定了它们的地位:尽可能少的在线滤波,尽可能多的离线滤波。用Steven的话解释就是,离线滤波还能“unfilter”。翻译翻译就是,在线滤波不可撤销,离线滤波还能撤销。
因此,在采集数据的时候可以不滤波,采集后我们再徐徐图之。
Tips:数据处理前,可以先明确我们的频谱范围,例如1~30Hz,这样就可以删除无关频带,加快计算速度。
2)滤波对波形的影响
模拟滤波器和数字滤波器在消除给定噪声频带后对波形的影响。
2.ICA消除脑电干扰成分:去除眼电、肌电等噪声
1)什么是ICA?
ICA,即independent component analysis,或者独立成分分析(ICA)。这里不会详细讲述技术原理,有需要的可以点击该链接。我会应用层面介绍ICA是个什么样子。
首先,我们需要明白脑电数据的数学模型——线性变换模型。如下图,图中记录到了三个脑电数据:E1,E2,E3。我们可以清晰地看见这三个脑电数据是源波形(source waveform)C1、C2和C3的加权和。由于容积传导效应,权重的绝对值小于1。这就是我们熟知的脑电数据的数学模型——任意EEG波形都是各个源波形的加权和。权重与响应脑区的响应程度和传播距离有关。 脑电数据处理与分析的目的是分离目标源波形,消除干扰波形,ICA恰恰符合这一需求。
图2。脑电成分的线性混合模型示例
2)如何做ICA?
- 大连理工大学丛丰裕和王小宇开发了名为ASAP_ICA_EEG_ArtifactRecjection的工具箱,该工具箱可以执行ICA算法,提取并识别脑电成分,且能剔除伪迹,
- 将该工具箱引入MATLAB路径目录后,就可以通过命令行来使用该工具箱。其实,就是运行脚本,没有专用的GUI供用户使用。
图3呈现了该工具箱的文件列表,标注了我们需要的两个代码文件。文件1和文件2分别完了ICA成分提取和伪迹成分剔除。具体操作既可以参照下面,也可以参考手册“Manual”。
图3。ASAP_ICA_EEG_ArtifactReject工具箱的文件列表
当我们运行Run_Lcasso脚本的时候,我们会看见图4弹出的窗口。此时,我们需要按照提示,设置ICA的参数。第一个参数“Number of Selected Independent Components”表示ICA提取出的成分的数量,第二参数“Number of running times”表示该脚本迭代拟合的次数,第三个参数,你可以懂的。如果不太明白,我们可以按照默认参数运行脚本,然后再修改。反正,我们的原始数据不会改变。当我们做完ICA后,数据会以Spatial_filtered_ICA.set重新保存,不用担心修改原始数据,当然这是两个文件都运行后才有的。
图4。ASAP_ICA_EEG_ArtifactRejection执行ICA的步骤示意图
当我们完成上面的步骤,它会弹出一系列的窗口,其中最有价值的是图4所示的窗口,它显示了ICA提取了哪些成分,并给出了该成分的概率,十分直观。其中,显示的编号十分重要,因为它将是我们剔除该成分的操作对象。
图5。ICA成分提取后的成分类别与其概率
当我们得到图5,说明我们已完全掌握文件2的信息。此时,我们需要运行文件1。如果两次运行选择的.set文件都是一样的,我们就应该得到下面图5的弹窗。它显示了所有伪迹的编号。实际上,我不建议(包括文档说明)直接点击OK,这可能在结果上出现很多的问题,例如突兀的波形。我建议(包括文档说明)只选择最明显的伪迹,然后重复上述步骤,直到没有伪迹,即每次只选择一个伪迹,重复循环知道没有伪迹。
图6。脚本Remove_ICs_Artifact.m运行后的界面
上述应该是使用该工具箱执行ICA和剔除伪迹的一个完整流程。两个步骤都会显示很多信息或者指标,这些信息或指标不是该工具箱发明的,如果我们想要去深究它们,只需要更加专业的知识。
3. 音乐特征提取(包括帧数的计算)
1)MIR
MIR, 即music information retrieval,音乐特征提取。Petri Toiviainen教授在他的MIRToolbox工具箱手册里给出了5个方面的音乐特征指标:dynamics(动态),rhythm(音律),timbre(音色),pitch(音高)和tonality(音调),具体可以参考MIRToolbox的用户手册(manual)。
2)MIRToolbox简要说明
实际上,MIRToolbox比上面的工具箱更加友好,文档内容十分详尽。这里我给出了建议是:选定某个指标后,直接查看对应函数如何调用,参数的设置参照“先默认,后精进”,不需要花里胡哨的设置诸多参数也能达到预期效果。一切从简!
3)帧数的计算
经过对模拟信号采样,数字信号由有限个样本点组成,但是我们一般不会选择整段信号来做数据处理。一般而言,我们会选择一个窗(window,例如时频分析常提及的矩形窗)。**一个窗截取部分样本点,进行频域转换,得到一组频域信号以及一个时间点(一般为窗中心点),我们通常把该时间点称为帧(frame),它代表了时间信息。**这样,把窗口从左往右滑动,就得到了整个样本长度的信息序列。整个信息序列由每次窗口得到的帧组成,窗口的宽度和移动距离决定了帧的数量,或者说序列的长度。
图7。MIRToolbox帧数的计算公式。
其中,numFrames 是帧数,Noverlap是窗口滑动重叠的样本点数。
4. EEG张量构建和分解(包括帧数的计算)
1)数据结构形式介绍
首先,我们需要明白什么是张量。通常来说,张量是一种数据结构形式。对于单个数据,例如11.11,我们称之为scalar,即标量;对于一维数据,例如[6.18, 11.11],我们称之为vector, 即向量;对于二维数据,例如[6.18, 11.11; 5.1, 10.1],我们称之为matrix,即矩阵;对于三维及以上,例如[[6.18, 11.11; 5.1, 10.1], [6.18, 11.11; 5.1, 10.1]],我们称之为tensor,即张量。
2)构建EEG张量
如果我们了解过深度学习,就会发现它比较喜欢以一定形式去组织张量结构,例如[sample, feature1,feature2,label]。这里我们需要构建的张量形式是[channel, frequency, frame]。frame是我们说的帧。
这其实就是一个时频分析技术,对应于图1,我们应该走到了STFT(短时傅里叶变换)步骤。在MATLAB中,spectrogram可以完成STFT。
3) 帧数的计算
理解帧数的计算十分重要,它关乎你能否使用此数据进行后续的相关分析的矩阵操作。公式如下图。尽管下面的公式和上面的公式不一样,但是它们其实是等效。如果把Noverlap转换为windowLength的百分比系数。那么,两个公式的系数之和等于1。
证明过程可以查看附录1。
图8。帧数计算公式
【注意】张量形式中frame代表的是time维度,之所以不用time而用frame,是因为frame更能表现时频分析技术的特点,即使用窗函数的中心点时刻位置,作为该段时间的频域信号的时间估计。
4)张量分解
这里使用的工具箱来自王德庆主页,也可以从GitHub上下载该工具箱。该工具箱相比于其他张量分解工具箱,具有运行速度快的特点(10倍)。下载后,我们应该包含下面的文件和文件夹。
图9。王德庆张量分解工具箱的文件列表
此外,我们可能需要安装tensor_toolbox来支撑该工具箱。
图10。MATLAB添加Tensor Toolbox
以次执行Example_iAPG_OngoingEEG_Step1和Example_iAPG_OngoingEEG_Step2文件后,我们应该得到至少一张这样的图:
图11。高相关的脑成分可视化图
图中描述了某成分的脑地形图,频谱图以及该成分与高相关音乐特征的时域图,CC给出了相关系数。
5. 脑电成分与音乐特征的相关分析:脑电张量成分的空间分布分量
在提取到有效脑电成分序列和音乐特征序列后,它们的序列长度应该一致,且在时刻上一一对应。这样,把两段序列做相关后,就可以查看它们在时间上的相关性,即如果该成分由某音乐特征诱发,那么该成分就应该伴随该音乐特征变化,随机性显著降低,表现在高相关系数值。
计算公式如下:
图12。相关系数公式
6. 可视化分析:脑地形图+频谱图+时序图
脑地形可以直观地呈现响应脑区,这对于快速定位脑成分的源区来说十分有帮助。也许用不着我们来绘图,因为张量分解的工具箱会绘制结果图,显示了高相关性的脑地形图+频谱图+成分与特征的时序图。
下面这是一个绘制32导的脑地形的范例。第一个参数是30个通道的值,第二参数是.locs文件的路径。这里为什么是30个通道而不是32?这是因为32导系统中包含了一个GND电极和一个参考电极。正是因为有了它们,我们才能得到经过参考电极重参考后的电压差。
topo = topoplot(corr_channel(:),locs_fullpath,'plotchans',[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30]);
总结
从数据处理的角度来说,自然情景范式的相关性分析处理并不复杂——提取成分序列和特征序列,执行相关操作,借助脑地形图+频谱图+时序图进行可视化分析。附录给出了构建张量和提取音乐特征的代码,有需要的道友可以尝试调试一下。
代码中选择的音乐特征是中心化波动值(fluctuation centroid),波动熵(fluctuation entropy),音阶清晰值(key clarity),脉冲清晰值(pulse clarity)和调值(mode),可在MIRToolbox工具箱中查看提取方式。关于原因,详见参考5。
参考
- Introduction to Event-Related Potential Technique, second edition. Steven J. Luck. MIT Press.
- 脑电EEG/ERP+神经脑科学、认知脑科学研究推荐必看数目(附带所有书籍高清pdf)
- 独立成分分析(ICA)
- 浅析张量分解(Tensor Decomposition)
- 不同意识状态对音乐感知的差异性:基于音乐特征与脑电张量分解的研究
附录
1. MIRToolbox和STFT帧数计算公式等效性证明
2. ICA工具箱使用说明
3. 代码
third_order_tensor_constrction.m
clc
clear
%close all
tic
%%
sampling = 1000;
sp = 3000;
window_length = sp;
n_overlap = 2*window_length/3;
NFFT=1000;
fs=msgbox('please select .mat dataset','selecting');
set(fs, 'position', [200 200 500 80]);
ft = findobj(fs, 'Type','text');
set(ft, 'FontSize', 30, 'Unit', 'normal');
uiwait(fs)
[filename,pathname]= uigetfile('*.mat','Select Datasets',...
'MultiSelect','on');
filename = cellstr(filename);
temp=load(fullfile(pathname,filename{1}));
matname=char(fieldnames(temp));
EEG_data=getfield(temp,matname);
fprintf('Size of EEG_data: %d*%d \n',size(EEG_data));
D=double(EEG_data);
d =D(16,:);% choose channel
figure
spectrogram(d,window_length,n_overlap,NFFT,sampling);
title('channel1 spectrum');
hold on
for channel = 1:size(EEG_data,1)
d =D(channel,:);
[S,F,T,P]=spectrogram(d,window_length,n_overlap,NFFT,sampling);
% S:[frequency, time]
temp1=abs(S);
third_order_tensor(channel,:,:)= temp1;
end
third_order_tensor = third_order_tensor(:,1:30,:);
fprintf('Size of tensor:%d*%d*%d\n',size(third_order_tensor));
save('third_order_tensor','third_order_tensor');
music_features_extraction_5.m
clc;
clear;
%close all
tic
%%
[music_filename,music_path]= uigetfile('*.mp3','Select a MUSIC',...
'MultiSelect','on');
sampling = 1000;
sp = 3000;
%% Music: Load a music
audio=miraudio(fullfile(music_path,music_filename));
audio_resample = miraudio(audio,'Sampling',sampling);
audio_data = get(audio,'Data');
audio_data = audio_data{1}{1};
audio_resample_data = get(audio_resample,'Data');
audio_resample_data = audio_resample_data{1}{1};
fprintf('Size of audio with fs=44100Hz: %d*%d\n',size(audio_data));
fprintf('Size of audio_rsample with fs=1000Hz: %d*%d\n',size(audio_resample_data));
%% Music: Extracting a given period of data
start_point = 1;
end_point =291236;
%end_point = size(EEG_data,2);
audio_cut = miraudio(audio_resample,'Extract', start_point, end_point,'sp','Start');
fprintf('Finished to extracte data from %d to %d \n',start_point, end_point);
audio_cut_data = get(audio_cut,'Data');
audio_cut_data = audio_cut_data{1}{1};
fprintf('Size of audio_cut after period cutting: %d*%d\n',size(audio_cut_data));
%% Music:frame decomposition
hop_rate = 100/3;
frames=mirframe(audio_cut,'Length',sp,'sp','Hop',hop_rate,'%')
fprintf('Frame decomposition with sp=%d and hop_rate=%d% \n',sp,hop_rate);
frames_data = get(frames,'Data');
frames_data = frames_data{1}{1};
framenumber=size(frames_data,2);
fprintf('Size of frames::%d*%d\n',size(frames_data));
%% feature: Tonal key clarity
% Key clarity: measure of the tonal clarity
% spetrum_keys = mirspectrum(audio_cut,'Frame',sp,'sp',hop_rate,'%')
spectrum_keys = mirspectrum(audio_cut,'Frame',sp,'sp',hop_rate,'%')
chromagram_keys = mirchromagram(spectrum_keys);
keystrength_keys = mirkeystrength(chromagram_keys);
[keys, clarity, keystrength] = mirkey(keystrength_keys);
keyclarity = get(clarity,'Data');
keyclarity = keyclarity{1}{1};
keyclarity = keyclarity.';
fprintf('Size of clarity: %d*%d \n',size(keyclarity));
%% feature:Tonal mode
spectrum_mode = mirspectrum(audio_cut,'Frame',sp,'sp',hop_rate,'%');
mode_data = mirmode(spectrum_mode);
mode =get(mode_data,'Data');
mode =mode{1}{1};
mode= mode.';
fprintf('Size of mode: %d*%d \n',size(mode));
%% feature: fluctuation centroid;fluctration entropy;pulse clarity
Fluctuation_Centroid=zeros(framenumber,1);
FluctuationEntropy=zeros(framenumber,1);
Pulse_Clarity=zeros(framenumber,1);
for i=0:framenumber-1
temp_audiocut=miraudio(audio_resample,'Extract',i,i+3);
temp_frame= mirframe(temp_audiocut,'Length',sp,'sp','Hop',hop_rate,'%');
flucdata=mirfluctuation(temp_frame,'summary');
flucdata=mirgetdata(flucdata);
%fluctuation centroid
temp_framenumber=size(flucdata,2);
temp_windowlength=size(flucdata,1);
m=((10/(temp_windowlength-1))*(0:temp_windowlength-1))';
fctemp1=flucdata(:,1);
fctemp1=sum(m.*fctemp1)/(sum(fctemp1)+eps);
Fluctuation_Centroid(i+1,1)=fctemp1;
%fluctration entropy
fetemp2=abs(flucdata);
summat=repmat(sum(fetemp2)+repmat(eps,...
[1 size(fetemp2,2) size(fetemp2,3) size(fetemp2,4)]),...
[size(fetemp2,1) 1 1 1]);
fetemp2=fetemp2./summat;
value=-sum(fetemp2.*log(fetemp2+eps))./log(size(fetemp2,1));
Fluctuation_Entropy(i+1,1)=value;
%pulse clarity
Pulse_Clarity(i+1,1)=mirgetdata(mirpulseclarity(temp_frame));
end
%% plot
figure;plot(1:framenumber,keyclarity);xlabel('Frame Number');ylabel('Value');title('Tonal Key Clarity');
figure;plot(1:framenumber,mode);xlabel('Frame Number');ylabel('Value');title('Tonal Mode');
figure;plot(1:framenumber,Fluctuation_Centroid);xlabel('Frame Number');ylabel('Value');title('Fluctuation Centroid');
figure;plot(1:framenumber,Fluctuation_Entropy);xlabel('Frame Number');ylabel('Value');title('Fluctuation Entropy');
figure;plot(1:framenumber,Pulse_Clarity);xlabel('Frame Number');ylabel('Value');title('Pulse Clarity');
%%
long_term_features=zeros(framenumber,5);
long_term_features(:,1)=keyclarity;
long_term_features(:,2)=mode;
long_term_features(:,3)=Fluctuation_Centroid;
long_term_features(:,4)=Fluctuation_Entropy;
long_term_features(:,5)=Pulse_Clarity;
music_feature_names={'keyclarity','mode','Fluctuation_Centroid','Fluctuation_Entropy','Pulse_Clarity'};
%%
save mirLongTermFeatures long_term_features music_feature_names
toc