本文旨在結合小波分析與經驗模態分解(EMD), 充分提取音樂干預下的腦電(EEG)信號特征參數, 提高情緒狀態評估的分類準確率與可靠性, 以期為輔助音樂治療提供支持與幫助。采用音樂誘發情緒的多通道標準情感數據庫(DEAP)中的數據, 利用小波變換提取出額區(F3, F4)、顳區(T7, T8)和中央(C3, C4)通道的α波、β波以及θ波節律; 對提取的腦電節律進行EMD以獲得固有模態函數(IMF)分量, 再進一步提取腦電節律波的IMF分量平均能量和幅度差特征值, 即每種節律波中包含3個平均能量特征和2個幅度差特征值, 以達到充分提取EEG特征信息的目的; 最后基于支持向量機分類器實現情感狀態評估。結果表明, 利用該算法可以使無情緒、積極情緒、消極情緒之間分類最優正確率達到100%, 使得積極與消極情緒之間的識別率提升10%左右, 可以實現無情緒與積極、無情緒與消極情緒等情感狀態的有效評估。處于不同情感狀態下, 音樂治療效果差異較大, 提高情感狀態評估的分類正確率, 將幫助提高音樂治療的效果, 更好地為音樂治療提供支持。
引用本文： 李昕, 田彥秀, 侯永捷, 齊曉英, 孫小棋, 范夢頔, 蔡二娟. 小波變換結合經驗模態分解在音樂干預腦電分析中的應用. 生物醫學工程學雜志, 2016, 33(4): 762-769. doi: 10.7507/1001-5515.20160124 復制
引言
音樂治療是一門融合音樂、醫學、心理學等多領域的新興、邊緣學科，以無痛苦、無不良反應、易管理等特點在音樂胎教、兒童自閉癥、抑郁癥及腦損傷、阿爾茨海默癥等神經精神性疾病的心理康復治療等方面取得了一定的研究成果[1-3]。
腦電(electroencephalogram, EEG)信號能夠綜合反映人體生理與心理狀態，結合腦電信號的變化分析情感狀態，進而評估音樂治療的效果是一種有效的評估方法。Karthick等[4]通過采集安靜閉眼休息、聆聽Indian Carnatic古典音樂和搖滾音樂這三種狀態下的額葉腦電信號，基于去趨勢波動分析(detrended fluctuation analysis, DFA)算法和多尺度熵(multiscale entropy, MSE)，表明有、無音樂時的腦電信號有顯著差異，大腦處理音樂時腦波的復雜性增大。Mikutta等[5]通過受試者閉眼聆聽貝多芬第五交響曲時記錄腦電信號、不記錄腦電信號而進行主觀喚醒度評估及問卷調查，基于譜功率協方差圖，發現處于高興奮期時，右額葉區域抑制低頻α波活動。Kabuto等[6]提取音樂腦電信號基本節律的頻域特征，發現愉快情緒狀態下，腦電θ波的能量升高，而α波的能量卻明顯下降。王小甜[7]基于時域幅度直方圖、頻域功率譜密度、時頻希爾伯特譜熵，定量、定性分析不同音樂情感下的腦電信號變化，發現額顳區θ節律與抑制性情緒成正比，而β節律與激勵性情緒成正比。賴永秀等[8]發現額區正性情緒的α波功率相對大于負性情緒。陳曾等[9]對情感腦電信號進行小波包分解重構，將β節律用于情感狀態識別。不同的情緒所引起的某種節律波活躍程度不同，因此，利用小波變換將原始信號分解重構獲得基本節律波，提取節律波的腦電特征將能更好地分析音樂情緒腦電的變化。
基于腦電評估音樂治療效果過程中，腦電特征提取是關鍵問題之一。腦電信號特征提取分析方法主要包括：①時頻分析法，如小波分析、短時傅里葉變換(short-time Fourier transform, STFT)、Cabor變換、Wigner-Ville分布(Wigner-Ville distribution, WVD)、Cohen類雙線性時頻分布等；②能量分析法，如幅度能量分析、平均能量分析、小波系數能量特征分析等；③基于非線性分析方法，如分維數、DFA、希爾伯特-黃變換、近似熵、LZ(Lempel-Ziv, LZ)復雜度、經驗模態分解(empirical mode decomposition, EMD)等。EMD是一種特別適用于處理非平穩、非線性信號的分析方法，該方法能根據信號本身的尺度特征對信號進行分解，獲得一系列的固有模態函數(intrinsic mode function, IMF)分量，信號的非線性和非平穩特征能在各階IMF分量中顯示出來[10]。
1998年，Huang等[11]提出EMD算法，應用于腦電信號分析中得到較好的效果[12]。隨后EMD算法迅速得到發展。例如，謝平等[13]將EMD與多尺度熵結合；余煒等[14]利用EMD分解提取IMF分量特征；楊航與郭曉金將EMD與能量特征結合均用于運動想象腦電數據的特征提取中，并且得到了較好的分類效果；Diez等[15]對有關心理任務誘發腦電進行EMD，對每個IMF分量提取出均方根、方差、Shannon熵、LZ復雜度、中心頻率和最大頻率6個特征，得到180組特征向量，然后利用Wilk的λ參數進行特征降維，將得到的最相關的特征向量作為腦電特征。Parvez等[16]將發作期和發作間期的癲癇腦電進行如下處理：①離散余弦變換(discrete cosine transform, DCT)后分離出高頻組分，提取出高頻部分的熵和能量特征；②EMD分解，計算出每個IMF分量的瞬時頻率和幅值，將其兩者結合產生帶寬特征。最后將DCT提取的特征和EMD特征結合用于腦電數據分類，達到了較好的分類結果。Arasteh等[17]把EMD靈敏度的信噪比作為腦電特征。楊鵬圓[18]利用小波包和希爾伯特-黃變換提取情感腦電特征進行了愉悅度識別。EMD在腦電信號的應用主要集中于運動想象腦電、疲勞腦電等，結合小波變換在情感腦電中的應用研究則較少，小波變換結合EMD進行信號去噪有很好的效果，因此，先用小波變換提取所需節律波，再結合EMD方法從情感腦電節律波中提取腦電特征用于情緒分類識別具有一定的研究意義。
本文以不同音樂刺激下的腦電狀態為分析對象，利用小波變換分解重構出與情緒相關程度較大的頻段——α波、β波和θ波，基于EMD的自適應性特色，充分利用所提取的節律波信號時域波形特征(局部極大、極小，過零點和均值線)，將α波、β波和θ波節律分別分解為系列IMF分量，進一步提取出IMF的平均能量和幅值差等特征，基于支持向量機(support vector machine, SVM)分類器實現情感狀態評估，進而為音樂治療提供幫助與支持。
1 數據
本文采用Koelstra等[19]提出的分析人類情感狀態的多通道標準情感數據庫(database for emotion analysis using physiological signals, DEAP)。該實驗數據記錄了32位健康受試者觀看40個時長為1 min的不同音樂視頻時的腦電信號和外周生理信號。受試者的平均年齡26.9歲，男女各半，實驗前被告知實驗的各個細節。該實驗是在照明受控的實驗室中進行，受試者坐在距離屏幕約1 m處，實驗前被問及音量是否合適，并且實驗中所有的視頻約占畫面的2/3，以便盡量減少受試者的不舒適感、眼動等干擾。
實驗中的音樂片段選擇是使用last.fin網站的情感標簽、視頻高亮監測和在線評測工具檢索刺激材料。首先，在last.fin網站利用情感標簽選取每種情感類型中推薦的前10首歌，共得到1084首歌，把這1084首歌按其所屬的情感類型分布到“喚醒度-效價”的二維情感模型中去，根據相關標準每個象限選出15首，另外每個象限再人工選擇15首，共選取120個。然后利用線性回歸方法找到喚醒度和情緒效價評估權值最高的1 min片段，得到包含最大量情感內容的120個片段。最后，使用網上在線測評方法對120個片段進行評估，選出實驗用的40個片段。
每個受試進行40次實驗。實驗開始前有2 min的基線校正(屏幕上出現一個固定光標)，讓受試心情放松，然后開始正式實驗：①顯示2 s音樂片段編號，提示受試者實驗進程；②看屏幕上的固定光標采集5 s基線；③播放1 min音樂片段，記錄所需信號；④受試者根據自我評估系統(self-assessment manikins, SAM)對該音樂片段的喚醒度、效價、喜歡程度以及優勢度進行評估。實驗進行到一半時休息一會兒再繼續。該實驗記錄的腦電信號采樣頻率為512 Hz，腦電信號的采集使用32導AgCl電極，電極的分布根據國際通用的標準10-20系統，同時采集8個電極的外周生理信號，共40個通道。
本文所選的腦電信號實驗數據是已經過預處理(降采樣，去除眼電等噪聲，經過4.0~45 Hz的濾波等)的數據，采樣頻率為128 Hz。每個受試的數據由實驗數據和標簽組成，記錄了時長63 s共計8 064個數據點。這63 s數據中，后60 s為實驗過程中的數據，前3 s為實驗中用于去除基線漂移的數據，因此我們在使用數據時應采用從第4秒開始的后60 s的數據段。每個受試數據形式為40×40×8 064，其中第一個40代表音頻，第二個40代表腦電信號常用通道。本文分析的數據是從某個音頻中的某個腦電通道選取的。
2 方法
2.1 算法流程圖
本文算法總流程如圖 1所示。
圖1
算法流程圖
Figure1.
Algorithm flow chart
2.2 EMD算法
EMD是一種基于自身的時間尺度對非線性、非平穩序列信號進行分解的算法，具有自適應性，不需要設置基函數。EMD的目的就是獲取有限的一系列IMF分量，得到的各階IMF分量隨尺度的增大頻率逐漸降低，有助于突出腦電信號各分量的局部特征。
EMD原理具體步驟如下，假設原始信號為x(t)：
步驟1：找出x(t)信號的所有局部極大值點和極小值點；
步驟2：對步驟1中的所有局部極大值點進行曲線擬合成上包絡線m1(t)，所有極小值點擬合成下包絡線m2(t)；
步驟3：求出上下包絡線的均值曲線u(t)：
| $u\left( t \right) = \frac{1}{2}\left[{{m_1}\left( t \right) + {m_2}\left( t \right)} \right]$ |
步驟4：從原始信號x(t)中分離出包絡均值曲線u(t)，得到一個剩余函數h1(t)：
| ${h_1}\left( t \right) = x\left( t \right) - u\left( t \right)$ |
如果該剩余函數h1(t)滿足IMF分量的兩個條件：Ⅰ.在整個信號段內的過零點數和極值點數最多相差為1，即在IMF中不能出現大于零的極小值，也不能出現小于零的極大值點；Ⅱ.信號的上下包絡線關于時間軸成局部對稱。則該剩余函數h1(t)為第一個IMF分量，否則h1(t)將作為新的原始信號重復步驟1~4，循環K次直到得到的剩余函數滿足IMF分量的兩個條件為止。此時，該剩余函數作為第一個IMF分量，記作：
| ${s_1}\left( t \right) = {h_{1k}}\left( t \right)$ |
本次循環不能為無限循環，Huang最終給出一個類似柯西收斂準則的分量停止條件：
| $SD = \frac{{\sum\limits_{t = 0}^T {{{\left| {{h_{k - 1}}\left( t \right) - {h_k}\left( t \right)} \right|}^2}} }}{{\sum\limits_{t = 0}^T {h_{k - 1}^2\left( t \right)} }}$ |
一般SD取值在0.1~0.3時，迭代停止，篩分過程結束。
步驟5：將IMF分量s1(t)從原始信號x(t)中分離出來，得到一個殘余信號c1(t)：
| ${c_1}\left( t \right) = x\left( t \right) - {s_1}\left( t \right)$ |
步驟6：再將殘余信號c1(t)作為新的原始信號重復步驟1~5，循環n次直到最終得到的殘余函數cn(t)為一個單調函數或常數為止，EMD分解過程結束，由上述公式可易得出：
| $x\left( t \right) = \sum\limits_{i = 1}^n {{s_i}\left( t \right) + {c_n}\left( t \right)} $ |
EMD分解得出各階IMF分量頻率不同，越晚分離出來的IMF分量頻率越低，因此可以根據音樂對腦電節律波頻率的影響，在各階IMF分量中提取出相應的特征。
2.3 特征提取
研究表明，音樂誘發情緒的腦電信號特征主要體現在α、β和θ這三種節律波上[20-23]。當人的情緒興奮時，β節律幅值升高，而悲傷時α節律下降，θ節律上升，即θ節律與抑制性情緒成正比，而β節律與激勵性情緒成正比[24-25]，并且在額葉與顳葉區域表現比較明顯，因此本研究選擇額區F3、F4、顳區T7、T8和中央C3、C4通道。
DEAP數據庫應用“效價-喚醒度”二維情感模型。該情感模型依據效價(valence)將情緒分為正、負兩極，位于正極的稱積極情緒，通常帶來愉悅感受，位于負極的稱消極情緒，通常產生不愉悅感受；同時依據喚醒度(arousal)區分情緒的強弱，喚醒度越大，所產生的情緒就越強烈。消極情緒一般包括悲傷、恐懼、憤怒、焦慮、痛苦、憎恨等[26]，積極情緒包括快樂、滿意、興趣、自豪、感激和愛等[27]。本文選擇四種比較常見的基本情緒，其中包括積極情緒happy、exciting，消極情緒sad、terrible。無情緒選擇安靜基線腦電信號。因60 s數據共7 680個數據點比較大，且MATLAB運行速度較慢，因此截取1 024個數據點進行分析。
首先將32個受試的5種情緒腦電信號利用小波變換，分解并重構出θ波、α波和β波；以happy狀態下的F4數據為例，其通道重構波如圖 2所示。
圖2
重構后的腦電三種節律波
Figure2.
Three electroencephalogram rhythms after reconstruction
將重構后的三種節律波進行EMD分解，得到的IMF分量如圖 3所示，從圖中可以看出happy情緒狀態下的F4通道數據經EMD可以得到7個IMF分量。再利用傅里葉變換轉換到頻域，得到各階IMF的頻譜圖如圖 4所示。
圖3
IMF分量圖
Figure3.
Components of IMFs
圖4
IMF分量頻譜圖
Figure4.
IMFs component spectrum
EMD得到的各階IMF分量頻率不同，由圖 4可以看出，IMF的階數越高，其頻率越低。如果提取所有IMF分量相關特征，那么得到的特征向量維數就會很高。這些特征量中就會包含很多與音樂腦電特征相關性較小的特征值，從而導致情緒識別的準確率降低。本文研究的腦電節律波頻率為4~30 Hz，因此截取前3階IMF分量如圖 5所示，這些IMF分量幾乎包含了信號90%的能量。取前三階IMF分量進行特征提取。
圖5
確定階數的IMF分量頻譜圖
Figure5.
Order number of IMFs spectrum component
由于音樂干預時，腦電的三種節律波的幅度均會發生變化，因此提取出相鄰IMF分量的幅值差作為一個特征值，即：
| ${H_{ij}} = \frac{1}{n}\left( {\sum\limits_{t = 1}^n {\left| {{s_i}\left( t \right) - {s_j}\left( t \right)} \right|} } \right)$ |
其中si(t)是第i個IMF分量，sj(t)是第j個IMF分量，n為IMF分量數據點個數，Hij為第i個IMF與第j個IMF分量的幅值差。各階IMF分量的頻率差異較大，即存在能量差，因此將各階IMF分量的平均能量作為一個特征值，即：
| ${E_l} = \frac{1}{n}\sum\limits_{t = 1}^n {{{\left| {{s_l}\left( l \right)} \right|}^2}} $ |
其中sl是第l個IMF分量，n為IMF分量數據點數。El是第l個IMF分量的平均能量特征。
綜上，本文每種類型音樂的每個通道的一個受試有3種節律波，每種節律波中包含3個平均能量特征值和2個幅度差特征值，因此每個受試有15個特征值。
2.4 LibSVM分類器
SVM分類器的目標是基于訓練數據產出一個模型，用來預測只給出屬性的測試數據的目標值。本文選擇SVM分類模型為C-支持向量分類(C-support vector classification, C-SVC)類型。
C-SVC分類模型算法的基本原理[28-29]：設給定的訓練集T={(Xi, yi)｜Xi∈Rn, yi∈{-1, 1}, i=1, …, k}，尋找函數h(X)，利用決策函數y=f(X)=sgn(h(x))來推斷出輸入向量X對應的輸出Y值。其具體步驟如下：
步驟1：選擇合適的核函數K(Xi, Xj)和懲罰參數C，構造并求解最優化問題
| $\mathop {\min }\limits_a \frac{1}{2}\sum\limits_{i = 1}^k {\sum\limits_{j = 1}^k {{y_i}{y_j}{\alpha _i}{\alpha _j}K\left( {{X_i},{X_j}} \right){\text{ - }}\sum\limits_{j = 1}^k {{\alpha _j}} } } $ |
| $\begin{gathered} s.t.\sum\limits_{i = 1}^k {{y_i}{\alpha _i} = } 0,\hfill \\ 0 \leqslant {\alpha _i} \leqslant C,i = 1,\cdots ,k \hfill \\ \end{gathered} $ |
得到最優解α*=(α1*, …, αk*)T；
步驟2：選擇α*的一個分量0＜αj*＜C，并由此計算b*=；
步驟3：構造決策函數
| $f\left( X \right) = \operatorname{sgn} \left( {\sum\limits_{i = 1}^k {{y_i}\alpha _i^*K\left( {{X_i},X} \right) + {b^*}} } \right)$ |
本文選擇高斯徑向基核函數(Gaussian radial basis function, RBF)作為核函數，其表達式為：，σ是寬度參數。本文選用臺灣林智仁(Chih-Jen Lin)教授2001年開發的LibSVM分類器工具箱實現音樂腦電的分類。利用粒子群優化算法(particle swarm optimization, PSO)進行參數尋優，得到最佳的懲罰因子c和核函數參數g。
3 結果及分析
3.1 無情緒與積極情緒、消極情緒之間
3.1.1 區分無情緒與積極情緒
如表 1所示，當無情緒狀態與積極情緒識別時，單一選用β波的平均能量和幅值差特征值進行分類比單一選用α波分類準確率要高，并且，選用這兩種波的特征值進行分類時，雖然最優準確率為100%，但測試集準確率相對于單一選用β節律降低了，說明腦電信號中的β節律波的特征與音樂誘發產生的積極情緒有關，并且單獨用β波比混合波的分類效果要好得多。
表1
無情緒與積極情緒分類表
Table1.
Classification of positive emotions and no emotions
3.1.2 區分無情緒與消極情緒
如表 2所示，可以看出當區別無情緒與消極情緒時，單獨使用α波特征不如單獨用θ波得到的準確率高，而這兩種波混合特征要比單獨用其中一種波進行分類的準確率高，由此可以說明，音樂誘發產生的消極情緒與α波和θ波都有關，但θ波與消極情緒的相關性要比α波的高。
表2
無情緒與消極情緒分類表
Table2.
Classification of negative emotions and no emotions
3.2 積極情緒與消極情緒之間
如表 3所示，可以看出本文算法在積極情緒之間以及消極情緒之間的特征分類準確率不高，與提取近似熵、小波熵特征得到的準確率相近，多在60%~70%范圍內；但對積極與消極情緒之間的特征分類準確率較高，多在70%~80%，由此說明利用提取腦電信號節律波IMF分量的特征進行分類的效果比較好。
表3
四種情緒之間兩兩分類表
Table3.
Dichotomous classification among four kinds of emotion
1985年，Cole和Ray分別研究β節律的高、低頻段與情緒加工的關系，結果表明積極性情緒與消極性情緒在β節律的高、低頻段所引起的大腦活動不同[30]。該研究只偏重于一種腦電節律。之后Kabuto和袁全等開始研究腦電基本節律的頻域特征與音樂之間的關系，對音樂所誘發的愉快放松情緒狀態下的腦電信號進行分析，發現θ節律波的能量升高，而α波的能量卻有明顯下降[31]。這些研究只是從腦電波的某一節律來分析情緒與腦電之間的聯系，或者只選擇某一情緒狀態下腦電信號作為研究對象，缺少廣泛性。而本文是選擇與情緒相關性較大的三種腦電節律波，分析在積極性情緒與消極性情緒這兩種狀態下節律波的變化，并得到了較好的分類效果，將情緒與腦電之間的研究進一步得到了擴展。
4 結論
本文結合小波分析與經驗模態分解實現了音樂干預腦電的特征提取。分別提取了腦電α、β以及θ波節律，并根據音樂刺激腦電特點，提取了相鄰IMF分量的幅值差，以及前3階IMF分量的平均能量等15個特征量。分別比較了無情緒與積極情緒、無情緒與消極情緒以及積極與消極情緒之間的分類正確率，結果表明，無情緒分別與積極情緒、消極情緒之間正確率可達90%以上；所選的四種情緒之間兩兩分類問題中，利用該腦電特征提取方法得到的分類準確率要比一般的特征提取方法(如近似熵、小波熵)高，可達到70%左右。但是，雖然四種情緒兩兩分類的識別率有所提高，但不能達到完全分離的效果，比如sad情緒與terrible情緒分類效果還是很低。因此該特征提取方法在同性近似情緒之間的識別效力還有待提升改進。
區分聆聽音樂時產生的情緒是積極還是消極，與腦電信號的節律波密切相關。當為積極情緒時β節律變化波動較大，提取β波的特征進行分類效果很好；當為消極情緒時α節律和θ節律波動較大，提取這兩種波的特征進行分類效果較好，并且θ波與消極情緒的相關性要比α波的相關性高。
不同類型的音樂對人情緒的影響不同，治療的結果也存在較大的差異，研究不同音樂下情感狀態的變化是一項有意義的研究，基于不同音樂干預下腦電情感狀態分析，將為音樂治療提供有效幫助。
引言
音樂治療是一門融合音樂、醫學、心理學等多領域的新興、邊緣學科，以無痛苦、無不良反應、易管理等特點在音樂胎教、兒童自閉癥、抑郁癥及腦損傷、阿爾茨海默癥等神經精神性疾病的心理康復治療等方面取得了一定的研究成果[1-3]。
腦電(electroencephalogram, EEG)信號能夠綜合反映人體生理與心理狀態，結合腦電信號的變化分析情感狀態，進而評估音樂治療的效果是一種有效的評估方法。Karthick等[4]通過采集安靜閉眼休息、聆聽Indian Carnatic古典音樂和搖滾音樂這三種狀態下的額葉腦電信號，基于去趨勢波動分析(detrended fluctuation analysis, DFA)算法和多尺度熵(multiscale entropy, MSE)，表明有、無音樂時的腦電信號有顯著差異，大腦處理音樂時腦波的復雜性增大。Mikutta等[5]通過受試者閉眼聆聽貝多芬第五交響曲時記錄腦電信號、不記錄腦電信號而進行主觀喚醒度評估及問卷調查，基于譜功率協方差圖，發現處于高興奮期時，右額葉區域抑制低頻α波活動。Kabuto等[6]提取音樂腦電信號基本節律的頻域特征，發現愉快情緒狀態下，腦電θ波的能量升高，而α波的能量卻明顯下降。王小甜[7]基于時域幅度直方圖、頻域功率譜密度、時頻希爾伯特譜熵，定量、定性分析不同音樂情感下的腦電信號變化，發現額顳區θ節律與抑制性情緒成正比，而β節律與激勵性情緒成正比。賴永秀等[8]發現額區正性情緒的α波功率相對大于負性情緒。陳曾等[9]對情感腦電信號進行小波包分解重構，將β節律用于情感狀態識別。不同的情緒所引起的某種節律波活躍程度不同，因此，利用小波變換將原始信號分解重構獲得基本節律波，提取節律波的腦電特征將能更好地分析音樂情緒腦電的變化。
基于腦電評估音樂治療效果過程中，腦電特征提取是關鍵問題之一。腦電信號特征提取分析方法主要包括：①時頻分析法，如小波分析、短時傅里葉變換(short-time Fourier transform, STFT)、Cabor變換、Wigner-Ville分布(Wigner-Ville distribution, WVD)、Cohen類雙線性時頻分布等；②能量分析法，如幅度能量分析、平均能量分析、小波系數能量特征分析等；③基于非線性分析方法，如分維數、DFA、希爾伯特-黃變換、近似熵、LZ(Lempel-Ziv, LZ)復雜度、經驗模態分解(empirical mode decomposition, EMD)等。EMD是一種特別適用于處理非平穩、非線性信號的分析方法，該方法能根據信號本身的尺度特征對信號進行分解，獲得一系列的固有模態函數(intrinsic mode function, IMF)分量，信號的非線性和非平穩特征能在各階IMF分量中顯示出來[10]。
1998年，Huang等[11]提出EMD算法，應用于腦電信號分析中得到較好的效果[12]。隨后EMD算法迅速得到發展。例如，謝平等[13]將EMD與多尺度熵結合；余煒等[14]利用EMD分解提取IMF分量特征；楊航與郭曉金將EMD與能量特征結合均用于運動想象腦電數據的特征提取中，并且得到了較好的分類效果；Diez等[15]對有關心理任務誘發腦電進行EMD，對每個IMF分量提取出均方根、方差、Shannon熵、LZ復雜度、中心頻率和最大頻率6個特征，得到180組特征向量，然后利用Wilk的λ參數進行特征降維，將得到的最相關的特征向量作為腦電特征。Parvez等[16]將發作期和發作間期的癲癇腦電進行如下處理：①離散余弦變換(discrete cosine transform, DCT)后分離出高頻組分，提取出高頻部分的熵和能量特征；②EMD分解，計算出每個IMF分量的瞬時頻率和幅值，將其兩者結合產生帶寬特征。最后將DCT提取的特征和EMD特征結合用于腦電數據分類，達到了較好的分類結果。Arasteh等[17]把EMD靈敏度的信噪比作為腦電特征。楊鵬圓[18]利用小波包和希爾伯特-黃變換提取情感腦電特征進行了愉悅度識別。EMD在腦電信號的應用主要集中于運動想象腦電、疲勞腦電等，結合小波變換在情感腦電中的應用研究則較少，小波變換結合EMD進行信號去噪有很好的效果，因此，先用小波變換提取所需節律波，再結合EMD方法從情感腦電節律波中提取腦電特征用于情緒分類識別具有一定的研究意義。
本文以不同音樂刺激下的腦電狀態為分析對象，利用小波變換分解重構出與情緒相關程度較大的頻段——α波、β波和θ波，基于EMD的自適應性特色，充分利用所提取的節律波信號時域波形特征(局部極大、極小，過零點和均值線)，將α波、β波和θ波節律分別分解為系列IMF分量，進一步提取出IMF的平均能量和幅值差等特征，基于支持向量機(support vector machine, SVM)分類器實現情感狀態評估，進而為音樂治療提供幫助與支持。
1 數據
本文采用Koelstra等[19]提出的分析人類情感狀態的多通道標準情感數據庫(database for emotion analysis using physiological signals, DEAP)。該實驗數據記錄了32位健康受試者觀看40個時長為1 min的不同音樂視頻時的腦電信號和外周生理信號。受試者的平均年齡26.9歲，男女各半，實驗前被告知實驗的各個細節。該實驗是在照明受控的實驗室中進行，受試者坐在距離屏幕約1 m處，實驗前被問及音量是否合適，并且實驗中所有的視頻約占畫面的2/3，以便盡量減少受試者的不舒適感、眼動等干擾。
實驗中的音樂片段選擇是使用last.fin網站的情感標簽、視頻高亮監測和在線評測工具檢索刺激材料。首先，在last.fin網站利用情感標簽選取每種情感類型中推薦的前10首歌，共得到1084首歌，把這1084首歌按其所屬的情感類型分布到“喚醒度-效價”的二維情感模型中去，根據相關標準每個象限選出15首，另外每個象限再人工選擇15首，共選取120個。然后利用線性回歸方法找到喚醒度和情緒效價評估權值最高的1 min片段，得到包含最大量情感內容的120個片段。最后，使用網上在線測評方法對120個片段進行評估，選出實驗用的40個片段。
每個受試進行40次實驗。實驗開始前有2 min的基線校正(屏幕上出現一個固定光標)，讓受試心情放松，然后開始正式實驗：①顯示2 s音樂片段編號，提示受試者實驗進程；②看屏幕上的固定光標采集5 s基線；③播放1 min音樂片段，記錄所需信號；④受試者根據自我評估系統(self-assessment manikins, SAM)對該音樂片段的喚醒度、效價、喜歡程度以及優勢度進行評估。實驗進行到一半時休息一會兒再繼續。該實驗記錄的腦電信號采樣頻率為512 Hz，腦電信號的采集使用32導AgCl電極，電極的分布根據國際通用的標準10-20系統，同時采集8個電極的外周生理信號，共40個通道。
本文所選的腦電信號實驗數據是已經過預處理(降采樣，去除眼電等噪聲，經過4.0~45 Hz的濾波等)的數據，采樣頻率為128 Hz。每個受試的數據由實驗數據和標簽組成，記錄了時長63 s共計8 064個數據點。這63 s數據中，后60 s為實驗過程中的數據，前3 s為實驗中用于去除基線漂移的數據，因此我們在使用數據時應采用從第4秒開始的后60 s的數據段。每個受試數據形式為40×40×8 064，其中第一個40代表音頻，第二個40代表腦電信號常用通道。本文分析的數據是從某個音頻中的某個腦電通道選取的。
2 方法
2.1 算法流程圖
本文算法總流程如圖 1所示。
圖1
算法流程圖
Figure1.
Algorithm flow chart
2.2 EMD算法
EMD是一種基于自身的時間尺度對非線性、非平穩序列信號進行分解的算法，具有自適應性，不需要設置基函數。EMD的目的就是獲取有限的一系列IMF分量，得到的各階IMF分量隨尺度的增大頻率逐漸降低，有助于突出腦電信號各分量的局部特征。
EMD原理具體步驟如下，假設原始信號為x(t)：
步驟1：找出x(t)信號的所有局部極大值點和極小值點；
步驟2：對步驟1中的所有局部極大值點進行曲線擬合成上包絡線m1(t)，所有極小值點擬合成下包絡線m2(t)；
步驟3：求出上下包絡線的均值曲線u(t)：
| $u\left( t \right) = \frac{1}{2}\left[{{m_1}\left( t \right) + {m_2}\left( t \right)} \right]$ |
步驟4：從原始信號x(t)中分離出包絡均值曲線u(t)，得到一個剩余函數h1(t)：
| ${h_1}\left( t \right) = x\left( t \right) - u\left( t \right)$ |
如果該剩余函數h1(t)滿足IMF分量的兩個條件：Ⅰ.在整個信號段內的過零點數和極值點數最多相差為1，即在IMF中不能出現大于零的極小值，也不能出現小于零的極大值點；Ⅱ.信號的上下包絡線關于時間軸成局部對稱。則該剩余函數h1(t)為第一個IMF分量，否則h1(t)將作為新的原始信號重復步驟1~4，循環K次直到得到的剩余函數滿足IMF分量的兩個條件為止。此時，該剩余函數作為第一個IMF分量，記作：
| ${s_1}\left( t \right) = {h_{1k}}\left( t \right)$ |
本次循環不能為無限循環，Huang最終給出一個類似柯西收斂準則的分量停止條件：
| $SD = \frac{{\sum\limits_{t = 0}^T {{{\left| {{h_{k - 1}}\left( t \right) - {h_k}\left( t \right)} \right|}^2}} }}{{\sum\limits_{t = 0}^T {h_{k - 1}^2\left( t \right)} }}$ |
一般SD取值在0.1~0.3時，迭代停止，篩分過程結束。
步驟5：將IMF分量s1(t)從原始信號x(t)中分離出來，得到一個殘余信號c1(t)：
| ${c_1}\left( t \right) = x\left( t \right) - {s_1}\left( t \right)$ |
步驟6：再將殘余信號c1(t)作為新的原始信號重復步驟1~5，循環n次直到最終得到的殘余函數cn(t)為一個單調函數或常數為止，EMD分解過程結束，由上述公式可易得出：
| $x\left( t \right) = \sum\limits_{i = 1}^n {{s_i}\left( t \right) + {c_n}\left( t \right)} $ |
EMD分解得出各階IMF分量頻率不同，越晚分離出來的IMF分量頻率越低，因此可以根據音樂對腦電節律波頻率的影響，在各階IMF分量中提取出相應的特征。
2.3 特征提取
研究表明，音樂誘發情緒的腦電信號特征主要體現在α、β和θ這三種節律波上[20-23]。當人的情緒興奮時，β節律幅值升高，而悲傷時α節律下降，θ節律上升，即θ節律與抑制性情緒成正比，而β節律與激勵性情緒成正比[24-25]，并且在額葉與顳葉區域表現比較明顯，因此本研究選擇額區F3、F4、顳區T7、T8和中央C3、C4通道。
DEAP數據庫應用“效價-喚醒度”二維情感模型。該情感模型依據效價(valence)將情緒分為正、負兩極，位于正極的稱積極情緒，通常帶來愉悅感受，位于負極的稱消極情緒，通常產生不愉悅感受；同時依據喚醒度(arousal)區分情緒的強弱，喚醒度越大，所產生的情緒就越強烈。消極情緒一般包括悲傷、恐懼、憤怒、焦慮、痛苦、憎恨等[26]，積極情緒包括快樂、滿意、興趣、自豪、感激和愛等[27]。本文選擇四種比較常見的基本情緒，其中包括積極情緒happy、exciting，消極情緒sad、terrible。無情緒選擇安靜基線腦電信號。因60 s數據共7 680個數據點比較大，且MATLAB運行速度較慢，因此截取1 024個數據點進行分析。
首先將32個受試的5種情緒腦電信號利用小波變換，分解并重構出θ波、α波和β波；以happy狀態下的F4數據為例，其通道重構波如圖 2所示。
圖2
重構后的腦電三種節律波
Figure2.
Three electroencephalogram rhythms after reconstruction
將重構后的三種節律波進行EMD分解，得到的IMF分量如圖 3所示，從圖中可以看出happy情緒狀態下的F4通道數據經EMD可以得到7個IMF分量。再利用傅里葉變換轉換到頻域，得到各階IMF的頻譜圖如圖 4所示。
圖3
IMF分量圖
Figure3.
Components of IMFs
圖4
IMF分量頻譜圖
Figure4.
IMFs component spectrum
EMD得到的各階IMF分量頻率不同，由圖 4可以看出，IMF的階數越高，其頻率越低。如果提取所有IMF分量相關特征，那么得到的特征向量維數就會很高。這些特征量中就會包含很多與音樂腦電特征相關性較小的特征值，從而導致情緒識別的準確率降低。本文研究的腦電節律波頻率為4~30 Hz，因此截取前3階IMF分量如圖 5所示，這些IMF分量幾乎包含了信號90%的能量。取前三階IMF分量進行特征提取。
圖5
確定階數的IMF分量頻譜圖
Figure5.
Order number of IMFs spectrum component
由于音樂干預時，腦電的三種節律波的幅度均會發生變化，因此提取出相鄰IMF分量的幅值差作為一個特征值，即：
| ${H_{ij}} = \frac{1}{n}\left( {\sum\limits_{t = 1}^n {\left| {{s_i}\left( t \right) - {s_j}\left( t \right)} \right|} } \right)$ |
其中si(t)是第i個IMF分量，sj(t)是第j個IMF分量，n為IMF分量數據點個數，Hij為第i個IMF與第j個IMF分量的幅值差。各階IMF分量的頻率差異較大，即存在能量差，因此將各階IMF分量的平均能量作為一個特征值，即：
| ${E_l} = \frac{1}{n}\sum\limits_{t = 1}^n {{{\left| {{s_l}\left( l \right)} \right|}^2}} $ |
其中sl是第l個IMF分量，n為IMF分量數據點數。El是第l個IMF分量的平均能量特征。
綜上，本文每種類型音樂的每個通道的一個受試有3種節律波，每種節律波中包含3個平均能量特征值和2個幅度差特征值，因此每個受試有15個特征值。
2.4 LibSVM分類器
SVM分類器的目標是基于訓練數據產出一個模型，用來預測只給出屬性的測試數據的目標值。本文選擇SVM分類模型為C-支持向量分類(C-support vector classification, C-SVC)類型。
C-SVC分類模型算法的基本原理[28-29]：設給定的訓練集T={(Xi, yi)｜Xi∈Rn, yi∈{-1, 1}, i=1, …, k}，尋找函數h(X)，利用決策函數y=f(X)=sgn(h(x))來推斷出輸入向量X對應的輸出Y值。其具體步驟如下：
步驟1：選擇合適的核函數K(Xi, Xj)和懲罰參數C，構造并求解最優化問題
| $\mathop {\min }\limits_a \frac{1}{2}\sum\limits_{i = 1}^k {\sum\limits_{j = 1}^k {{y_i}{y_j}{\alpha _i}{\alpha _j}K\left( {{X_i},{X_j}} \right){\text{ - }}\sum\limits_{j = 1}^k {{\alpha _j}} } } $ |
| $\begin{gathered} s.t.\sum\limits_{i = 1}^k {{y_i}{\alpha _i} = } 0,\hfill \\ 0 \leqslant {\alpha _i} \leqslant C,i = 1,\cdots ,k \hfill \\ \end{gathered} $ |
得到最優解α*=(α1*, …, αk*)T；
步驟2：選擇α*的一個分量0＜αj*＜C，并由此計算b*=；
步驟3：構造決策函數
| $f\left( X \right) = \operatorname{sgn} \left( {\sum\limits_{i = 1}^k {{y_i}\alpha _i^*K\left( {{X_i},X} \right) + {b^*}} } \right)$ |
本文選擇高斯徑向基核函數(Gaussian radial basis function, RBF)作為核函數，其表達式為：，σ是寬度參數。本文選用臺灣林智仁(Chih-Jen Lin)教授2001年開發的LibSVM分類器工具箱實現音樂腦電的分類。利用粒子群優化算法(particle swarm optimization, PSO)進行參數尋優，得到最佳的懲罰因子c和核函數參數g。
3 結果及分析
3.1 無情緒與積極情緒、消極情緒之間
3.1.1 區分無情緒與積極情緒
如表 1所示，當無情緒狀態與積極情緒識別時，單一選用β波的平均能量和幅值差特征值進行分類比單一選用α波分類準確率要高，并且，選用這兩種波的特征值進行分類時，雖然最優準確率為100%，但測試集準確率相對于單一選用β節律降低了，說明腦電信號中的β節律波的特征與音樂誘發產生的積極情緒有關，并且單獨用β波比混合波的分類效果要好得多。
表1
無情緒與積極情緒分類表
Table1.
Classification of positive emotions and no emotions
3.1.2 區分無情緒與消極情緒
如表 2所示，可以看出當區別無情緒與消極情緒時，單獨使用α波特征不如單獨用θ波得到的準確率高，而這兩種波混合特征要比單獨用其中一種波進行分類的準確率高，由此可以說明，音樂誘發產生的消極情緒與α波和θ波都有關，但θ波與消極情緒的相關性要比α波的高。
表2
無情緒與消極情緒分類表
Table2.
Classification of negative emotions and no emotions
3.2 積極情緒與消極情緒之間
如表 3所示，可以看出本文算法在積極情緒之間以及消極情緒之間的特征分類準確率不高，與提取近似熵、小波熵特征得到的準確率相近，多在60%~70%范圍內；但對積極與消極情緒之間的特征分類準確率較高，多在70%~80%，由此說明利用提取腦電信號節律波IMF分量的特征進行分類的效果比較好。
表3
四種情緒之間兩兩分類表
Table3.
Dichotomous classification among four kinds of emotion
1985年，Cole和Ray分別研究β節律的高、低頻段與情緒加工的關系，結果表明積極性情緒與消極性情緒在β節律的高、低頻段所引起的大腦活動不同[30]。該研究只偏重于一種腦電節律。之后Kabuto和袁全等開始研究腦電基本節律的頻域特征與音樂之間的關系，對音樂所誘發的愉快放松情緒狀態下的腦電信號進行分析，發現θ節律波的能量升高，而α波的能量卻有明顯下降[31]。這些研究只是從腦電波的某一節律來分析情緒與腦電之間的聯系，或者只選擇某一情緒狀態下腦電信號作為研究對象，缺少廣泛性。而本文是選擇與情緒相關性較大的三種腦電節律波，分析在積極性情緒與消極性情緒這兩種狀態下節律波的變化，并得到了較好的分類效果，將情緒與腦電之間的研究進一步得到了擴展。
4 結論
本文結合小波分析與經驗模態分解實現了音樂干預腦電的特征提取。分別提取了腦電α、β以及θ波節律，并根據音樂刺激腦電特點，提取了相鄰IMF分量的幅值差，以及前3階IMF分量的平均能量等15個特征量。分別比較了無情緒與積極情緒、無情緒與消極情緒以及積極與消極情緒之間的分類正確率，結果表明，無情緒分別與積極情緒、消極情緒之間正確率可達90%以上；所選的四種情緒之間兩兩分類問題中，利用該腦電特征提取方法得到的分類準確率要比一般的特征提取方法(如近似熵、小波熵)高，可達到70%左右。但是，雖然四種情緒兩兩分類的識別率有所提高，但不能達到完全分離的效果，比如sad情緒與terrible情緒分類效果還是很低。因此該特征提取方法在同性近似情緒之間的識別效力還有待提升改進。
區分聆聽音樂時產生的情緒是積極還是消極，與腦電信號的節律波密切相關。當為積極情緒時β節律變化波動較大，提取β波的特征進行分類效果很好；當為消極情緒時α節律和θ節律波動較大，提取這兩種波的特征進行分類效果較好，并且θ波與消極情緒的相關性要比α波的相關性高。
不同類型的音樂對人情緒的影響不同，治療的結果也存在較大的差異，研究不同音樂下情感狀態的變化是一項有意義的研究，基于不同音樂干預下腦電情感狀態分析，將為音樂治療提供有效幫助。
