自然情景范式下脑电与音乐特征的相关性分析：EEG张量分解，音乐特征提取MIRToolbox和topoplot脑地形图绘制与响应脑区

---

未经允许，禁止转载。

---

摘要

本问详细讲解了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？

1. 大连理工大学丛丰裕和王小宇开发了名为ASAP_ICA_EEG_ArtifactRecjection的工具箱，该工具箱可以执行ICA算法，提取并识别脑电成分，且能剔除伪迹，
2. 将该工具箱引入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。

---

参考

1. Introduction to Event-Related Potential Technique, second edition. Steven J. Luck. MIT Press.
2. 脑电EEG/ERP+神经脑科学、认知脑科学研究推荐必看数目(附带所有书籍高清pdf)
3. 独立成分分析（ICA）
4. 浅析张量分解（Tensor Decomposition）
5. 不同意识状态对音乐感知的差异性:基于音乐特征与脑电张量分解的研究

附录

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