據(jù)集與SVM的腦電情感識別：從預(yù)處理到模型調(diào)優(yōu)全流程實(shí)踐)
簡介本資源是一套面向腦機(jī)接口與情感計(jì)算方向研究者的完整論文代碼實(shí)現(xiàn)方案聚焦DEAP數(shù)據(jù)集上的四分類情緒識別任務(wù)效價(jià)-喚醒二維象限劃分適用于具備MATLAB基礎(chǔ)與機(jī)器學(xué)習(xí)入門知識的研究生及科研人員。資源包含16個文件以10個核心MATLAB腳本如頻帶濾波、DWT特征提取、SVM訓(xùn)練等、3個Jupyter Notebook含KNN對比實(shí)驗(yàn)、1篇PDF論文、1份README說明及輔助文本文件為主總大小僅488KB結(jié)構(gòu)緊湊、模塊清晰便于復(fù)現(xiàn)實(shí)驗(yàn)流程與特征工程細(xì)節(jié)。已有3636人學(xué)習(xí)下載讀者可直接運(yùn)行代碼完成從原始EEG信號預(yù)處理、α/β/θ頻帶分離、箱線圖特征篩選、統(tǒng)計(jì)特征偏度、峰度、波熵構(gòu)建到SVM四分類建模的全流程并獲得92.36%測試準(zhǔn)確率結(jié)果顯著優(yōu)于現(xiàn)有主流方法附帶可視化腳本支持特征分布與分類效果分析。1. 項(xiàng)目背景與核心目標(biāo)最近在整理過往的研究項(xiàng)目翻到了一個基于DEAP數(shù)據(jù)集做情感識別的老項(xiàng)目感覺挺有意思的。這個項(xiàng)目說白了就是通過分析人的腦電信號來判斷他/她當(dāng)前是高興、悲傷、平靜還是興奮。聽起來有點(diǎn)科幻但其實(shí)背后的技術(shù)路徑已經(jīng)相當(dāng)成熟。當(dāng)時我用的是經(jīng)典的機(jī)器學(xué)習(xí)方法核心流程就是從原始的腦電信號里提取特征然后用支持向量機(jī)SVM去訓(xùn)練一個分類器。整個過程在MATLAB里實(shí)現(xiàn)從數(shù)據(jù)預(yù)處理、特征工程到模型訓(xùn)練和評估一條龍搞定。今天這篇文章我就把這個項(xiàng)目的完整實(shí)現(xiàn)思路、關(guān)鍵代碼和踩過的坑系統(tǒng)地梳理一遍希望能給對腦機(jī)接口、情感計(jì)算或者信號處理感興趣的朋友提供一個清晰的參考模板。無論你是剛?cè)腴T的研究生還是想快速復(fù)現(xiàn)一個baseline的工程師這篇內(nèi)容應(yīng)該都能幫到你。DEAP數(shù)據(jù)集在情感計(jì)算領(lǐng)域算是“明星數(shù)據(jù)集”了它包含了32位被試者在觀看音樂視頻時的腦電EEG和外圍生理信號如肌電、皮電等并且對每段視頻都進(jìn)行了效價(jià)Valence、喚醒度Arousal、支配度Dominance和喜愛度Liking的評分。我們通常取效價(jià)和喚醒度這兩個維度將它們劃分為高/低兩類然后組合成四個象限就對應(yīng)了四種情感狀態(tài)高喚醒高效價(jià)HAHV可理解為興奮/高興、高喚醒低效價(jià)HALV可理解為焦慮/憤怒、低喚醒高效價(jià)LAHV可理解為平靜/放松、低喚醒低效價(jià)LALV可理解為悲傷/低落。我們的任務(wù)就是讓機(jī)器通過腦電信號認(rèn)出這四種狀態(tài)。2. DEAP數(shù)據(jù)集預(yù)處理與關(guān)鍵理解拿到DEAP數(shù)據(jù)集后第一步不是急著跑代碼而是先理解它的數(shù)據(jù)結(jié)構(gòu)。DEAP數(shù)據(jù)以MAT文件.mat格式提供通常包含兩個關(guān)鍵變量data和labels。data是一個4維數(shù)組維度通常是參與者 × 視頻片段 × 通道 × 數(shù)據(jù)點(diǎn)例如 32 x 40 x 32 x 8064。這里8064對應(yīng)的是63秒的視頻片段采樣頻率是128Hz所以63*1288064個數(shù)據(jù)點(diǎn)。labels則包含了每個視頻片段對應(yīng)的效價(jià)、喚醒度等評分。預(yù)處理的目標(biāo)是把原始的、充滿噪聲的腦電信號變成干凈、可用于特征提取的信號。這個過程至關(guān)重要直接決定了后續(xù)特征的質(zhì)量和模型的性能上限。2.1 數(shù)據(jù)加載與維度重塑首先我們需要把數(shù)據(jù)加載進(jìn)來并轉(zhuǎn)換成更易于處理的格式。原始的四維數(shù)組在操作時不太方便我們通常會把數(shù)據(jù)重塑成二維矩陣樣本 × 特征但這里先按被試和試驗(yàn)進(jìn)行分割。% 假設(shè)數(shù)據(jù)文件為 DEAP_data.mat load(DEAP_data.mat); % 這會加載變量 data 和 labels % 查看數(shù)據(jù)維度 [num_participants, num_trials, num_channels, num_samples] size(data); disp([數(shù)據(jù)維度: , num2str(num_participants), 參與者, , ... num2str(num_trials), 試驗(yàn), , num2str(num_channels), 通道, , ... num2str(num_samples), 采樣點(diǎn)]); % 提取標(biāo)簽這里以效價(jià)和喚醒度為例假設(shè)它們是labels的前兩列 valence_labels labels(:, 1); % 效價(jià)維度 arousal_labels labels(:, 2); % 喚醒度維度2.2 腦電信號預(yù)處理流程詳解腦電信號非常微弱容易受到眼電EOG、肌電EMG、工頻干擾50Hz等噪聲污染。標(biāo)準(zhǔn)的預(yù)處理流程包括以下幾步重參考原始腦電記錄通常是單極導(dǎo)聯(lián)需要一個參考點(diǎn)。DEAP數(shù)據(jù)已使用平均參考這一步通?？梢蕴^但了解其概念很重要。平均參考是假設(shè)所有電極電位的平均值為零將每個通道的信號減去所有通道的平均值有助于減少共同噪聲。帶通濾波保留有用的腦電節(jié)律成分。一般認(rèn)為與情感相關(guān)的腦電成分分布在Delta(1-4Hz), Theta(4-8Hz), Alpha(8-13Hz), Beta(13-30Hz), Gamma(30-45Hz)等頻段。我們通常用一個較寬的帶通濾波器比如0.5Hz到45Hz以保留大部分有效信息同時濾除極低頻的漂移和高頻的肌電噪聲。陷波濾波去除工頻干擾50Hz或60Hz取決于地區(qū)。在中國是50Hz。分段與基線校正DEAP數(shù)據(jù)已經(jīng)是觀看視頻期間的63秒連續(xù)數(shù)據(jù)。有時我們需要進(jìn)行更短時間窗的分析如分析事件相關(guān)電位ERP但在這個項(xiàng)目中我們直接使用整個63秒的片段作為一個樣本?；€校正通常指減去一段“靜息”或“刺激前”時期的平均幅值以消除個體直流偏移。DEAP數(shù)據(jù)在提供前可能已做過處理但為了嚴(yán)謹(jǐn)我們可以計(jì)算每個試驗(yàn)最開始一小段如前1秒數(shù)據(jù)的均值然后從整個試驗(yàn)數(shù)據(jù)中減去它。下面是用MATLAB實(shí)現(xiàn)濾波的示例代碼。我強(qiáng)烈推薦使用EEGLAB或FieldTrip這些專業(yè)的腦電處理工具箱它們功能強(qiáng)大且穩(wěn)定。但為了理解原理這里用MATLAB內(nèi)置函數(shù)演示% 定義采樣頻率和濾波器參數(shù) Fs 128; % DEAP采樣率 low_cutoff 0.5; % 高通截止頻率 (Hz) high_cutoff 45.0; % 低通截止頻率 (Hz) notch_freq 50; % 陷波頻率 (Hz) % 設(shè)計(jì)一個帶通濾波器例如使用巴特沃斯濾波器 order 4; % 濾波器階數(shù) [b_band, a_band] butter(order, [low_cutoff, high_cutoff]/(Fs/2), bandpass); % 設(shè)計(jì)一個陷波濾波器 wo notch_freq/(Fs/2); % 歸一化頻率 bw wo/35; % 帶寬 [b_notch, a_notch] iirnotch(wo, bw); % 假設(shè)我們處理一個被試一個試驗(yàn)的一個通道數(shù)據(jù) trial_eeg (1 x 8064) % 1. 帶通濾波 filtered_eeg filtfilt(b_band, a_band, trial_eeg); % 使用filtfilt實(shí)現(xiàn)零相位濾波 % 2. 陷波濾波 filtered_eeg filtfilt(b_notch, a_notch, filtered_eeg); % 注意對多通道數(shù)據(jù)需要循環(huán)處理每個通道。 % 更高效的做法是使用eeglab的pop_eegfiltnew函數(shù) % EEG pop_eegfiltnew(EEG, low_cutoff, high_cutoff); % EEG pop_eegfiltnew(EEG, notch_freq-1, notch_freq1, [], 1); % 陷波注意filtfilt函數(shù)執(zhí)行的是零相位濾波它通過對數(shù)據(jù)正向和反向各濾波一次消除了濾波器帶來的相位延遲這對于后續(xù)的時頻分析非常重要。但它的計(jì)算量是普通濾波的兩倍。2.3 標(biāo)簽生成從維度到四分類DEAP的原始標(biāo)簽是連續(xù)的效價(jià)和喚醒度評分1-9分。我們需要將其離散化為高/低兩類通常以中值5為界。% 假設(shè)valence_labels和arousal_labels是來自所有試驗(yàn)的向量 valence_threshold 5; arousal_threshold 5; % 生成二值標(biāo)簽 valence_binary valence_labels valence_threshold; % 1高效價(jià)0低效價(jià) arousal_binary arousal_labels arousal_threshold; % 1高喚醒0低喚醒 % 組合成四分類標(biāo)簽 (1: LALV, 2: LAHV, 3: HALV, 4: HAHV) % 注意這里的映射關(guān)系可以根據(jù)你的情感模型定義調(diào)整 four_class_labels zeros(size(valence_labels)); four_class_labels(~arousal_binary ~valence_binary) 1; % 低喚醒低效價(jià) four_class_labels(~arousal_binary valence_binary) 2; % 低喚醒高效價(jià) four_class_labels(arousal_binary ~valence_binary) 3; % 高喚醒低效價(jià) four_class_labels(arousal_binary valence_binary) 4; % 高喚醒高效價(jià) % 檢查類別分布 tabulate(four_class_labels)類別不平衡是情感識別中常見的問題。DEAP數(shù)據(jù)集的四個類別分布可能不是完全均勻的。如果差異較大在后續(xù)訓(xùn)練模型時需要考慮采用加權(quán)SVM、過采樣如SMOTE或欠采樣等策略。3. 時頻域特征提取從信號到信息特征提取是情感識別的核心。腦電信號在時域上看就是一條隨時間變化的電壓曲線信息隱藏很深。時頻分析能告訴我們信號的能量在不同頻率成分上是如何隨時間變化的這比單純的時域或頻域特征更能捕捉情感的動態(tài)特性。3.1 為何選擇時頻域特征情感變化不是瞬間切換的而是一個動態(tài)過程。例如從平靜到興奮大腦中Alpha波8-13Hz的能量可能會降低去同步化而Beta波13-30Hz的能量可能會增加。單純的頻域特征如整個片段的功率譜會丟失時間信息而單純的時域特征如均值、方差對頻率不敏感。時頻分析如小波變換、短時傅里葉變換提供了“時間-頻率”二維平面上的能量分布讓我們能夠提取諸如“在視頻播放的第20秒到30秒前額葉Theta波能量上升”這樣的動態(tài)特征這對情感識別至關(guān)重要。3.2 基于短時傅里葉變換STFT的特征提取實(shí)踐STFT是理解時頻分析最直觀的工具。它的思想很簡單把長信號分成許多短的、重疊的片段對每個片段做傅里葉變換從而得到每個時間點(diǎn)附近的頻率成分。% 假設(shè) filtered_eeg 是一個預(yù)處理后的單通道腦電信號 (1 x N) signal filtered_eeg; N length(signal); Fs 128; % 設(shè)置STFT參數(shù) window_length 256; % 窗長對應(yīng)2秒 (256/1282)。窗越長頻率分辨率越高時間分辨率越低。 noverlap 128; % 重疊點(diǎn)數(shù)通常為窗長的一半。重疊越多時頻圖越平滑。 nfft 512; % FFT點(diǎn)數(shù)通常大于等于窗長。決定頻率軸的精細(xì)度。 % 計(jì)算STFT [S, F, T] spectrogram(signal, window_length, noverlap, nfft, Fs); % S是復(fù)數(shù)矩陣維度為 (nfft/21) x 時間窗數(shù) % F是頻率向量 % T是時間向量每個窗的中心時間 % 計(jì)算功率譜密度 (PSD) P abs(S).^2; % 現(xiàn)在P就是一個時頻能量矩陣。我們可以從中提取各種特征。從時頻矩陣P中我們可以提取多種特征特定頻帶功率這是最常用的特征。先定義頻帶范圍然后對P矩陣中對應(yīng)頻率區(qū)間的功率進(jìn)行求和或平均。% 定義經(jīng)典頻帶邊界 (Hz) band_defs { Delta, 1, 4; Theta, 4, 8; Alpha, 8, 13; Beta, 13, 30; Gamma, 30, 45; }; % 初始化特征向量 band_powers []; for i 1:size(band_defs, 1) band_name band_defs{i, 1}; f_low band_defs{i, 2}; f_high band_defs{i, 3}; % 找到頻率向量F中在頻帶內(nèi)的索引 freq_indices (F f_low) (F f_high); % 計(jì)算該頻帶在所有時間窗上的總功率然后對所有時間窗取平均 % 也可以考慮取對數(shù)(log)使分布更接近正態(tài) band_power_total sum(P(freq_indices, :), 1); % 對頻率維求和得到每個時間窗的該頻帶功率 band_power_mean mean(band_power_total); % 對所有時間窗求平均得到一個標(biāo)量特征 % band_power_log log(band_power_mean eps); % 取對數(shù)加eps防止為0 band_powers [band_powers, band_power_mean]; end % 現(xiàn)在 band_powers 是一個1x5的向量包含5個頻帶的平均功率時頻統(tǒng)計(jì)特征我們可以把時頻矩陣P看作一幅圖像計(jì)算其統(tǒng)計(jì)特性。平均功率mean(P(:))功率標(biāo)準(zhǔn)差std(P(:))功率偏度/峰度skewness(P(:)),kurtosis(P(:))描述功率分布的形態(tài)。譜熵將每個時間點(diǎn)的功率譜歸一化為概率分布計(jì)算其香農(nóng)熵反映頻譜的復(fù)雜度或平坦度。% 計(jì)算譜熵 (Spectral Entropy) % 對每個時間窗的功率譜進(jìn)行歸一化 P_norm zeros(size(P)); for t 1:size(P, 2) p_slice P(:, t); p_slice p_slice / sum(p_slice); % 歸一化為概率分布 P_norm(:, t) p_slice; end % 計(jì)算每個時間窗的譜熵 spectral_entropy_per_window -sum(P_norm .* log2(P_norm eps), 1); % 取所有時間窗譜熵的平均值作為特征 mean_spectral_entropy mean(spectral_entropy_per_window);微分熵特征在腦電分析中微分熵Differential Entropy, DE被證明是有效的特征特別是在基于頻帶的特征提取中。對于一段服從高斯分布N(μ, σ^2)的信號其微分熵計(jì)算公式為0.5 * log(2πeσ^2)。在實(shí)際操作中我們通常假設(shè)每個頻帶內(nèi)的信號近似服從高斯分布因此可以用該頻帶功率的對數(shù)變換來近似微分熵。許多研究直接使用log(頻帶功率)作為特征效果很好。% 基于之前計(jì)算的 band_power_total (每個時間窗的頻帶功率) de_features []; for i 1:size(band_defs, 1) % 假設(shè) band_power_total 是之前計(jì)算好的第i個頻帶的功率時間序列 % 這里需要重新計(jì)算每個頻帶每個時間窗的功率 freq_indices (F band_defs{i, 2}) (F band_defs{i, 3}); band_power_per_window sum(P(freq_indices, :), 1); % 1 x 時間窗數(shù) % 計(jì)算該頻帶功率的微分熵近似為對數(shù)功率 % 可以對每個時間窗取對數(shù)后平均也可以先平均再取對數(shù)效果略有不同 de_per_window log(band_power_per_window eps); mean_de mean(de_per_window); de_features [de_features, mean_de]; end3.3 多通道特征融合與特征向量構(gòu)建上面我們只針對一個通道進(jìn)行了特征提取。DEAP有32個通道或64個我們需要對所有通道重復(fù)上述過程然后將所有特征拼接成一個長向量代表一個試驗(yàn)樣本。% 假設(shè)我們有一個函數(shù) extract_features_from_channel(eeg_signal, Fs) % 它輸入一個通道的信號輸出一個特征向量 (例如 5個頻帶功率 1個譜熵 6維) num_channels 32; features_per_trial []; for ch 1:num_channels single_channel_signal squeeze(data(participant_idx, trial_idx, ch, :)); % 獲取指定被試、試驗(yàn)、通道的數(shù)據(jù) single_channel_signal preprocess_eeg(single_channel_signal, Fs); % 預(yù)處理函數(shù) ch_features extract_features_from_channel(single_channel_signal, Fs); features_per_trial [features_per_trial, ch_features]; end % 最終features_per_trial 是一個 1 x (6 * 32) 1 x 192 維的特征向量。這就是一個樣本的特征向量。對數(shù)據(jù)集中的所有試驗(yàn)和所有被試重復(fù)此操作注意要按被試劃分訓(xùn)練測試集避免數(shù)據(jù)泄露我們就能得到特征矩陣X(樣本數(shù) x 特征數(shù)) 和標(biāo)簽向量y。實(shí)操心得特征維度會爆炸式增長通道數(shù) x 每通道特征數(shù)。32個通道每個通道提取6個特征就是192維。如果使用更復(fù)雜的時頻特征如Hjorth參數(shù)、分形維數(shù)等維度會更高。高維特征容易導(dǎo)致“維數(shù)災(zāi)難”并且包含大量冗余信息。因此特征選擇是必不可少的后續(xù)步驟。我們可以使用方差閾值、相關(guān)系數(shù)、遞歸特征消除RFE或基于模型如L1正則化的方法來篩選出最具有判別力的特征子集。在項(xiàng)目初期為了快速驗(yàn)證流程可以先用所有特征但最終一定要做特征選擇來優(yōu)化模型。4. 支持向量機(jī)分類器的構(gòu)建與調(diào)優(yōu)特征準(zhǔn)備好之后就進(jìn)入了建模階段。支持向量機(jī)SVM在小樣本、高維度的分類問題上表現(xiàn)優(yōu)異并且可以通過核函數(shù)處理非線性問題非常適合作為腦電情感識別的基線分類器。4.1 SVM原理簡述與核函數(shù)選擇SVM的核心思想是尋找一個最優(yōu)超平面使得兩類樣本到這個超平面的“間隔”最大化。對于線性不可分的數(shù)據(jù)SVM通過“核技巧”將數(shù)據(jù)映射到高維空間使其在高維空間中線性可分。常用的核函數(shù)有線性核K(x, z) x^T * z。參數(shù)少速度快可解釋性強(qiáng)。如果特征已經(jīng)足夠好或者特征維度很高線性核往往就能取得不錯的效果并且不容易過擬合。徑向基函數(shù)核K(x, z) exp(-γ * ||x - z||^2)。也叫高斯核是最常用的非線性核。它能夠處理非常復(fù)雜的非線性關(guān)系但需要調(diào)節(jié)兩個超參數(shù)懲罰系數(shù)C和核系數(shù)γ。γ定義了單個樣本的影響范圍γ越大影響范圍越小決策邊界越復(fù)雜容易過擬合γ越小決策邊界越平滑容易欠擬合。對于腦電情感識別由于特征與情感狀態(tài)之間的關(guān)系可能是非線性的RBF核通常是首選。但我的經(jīng)驗(yàn)是先用線性核試一下。如果線性核的準(zhǔn)確率已經(jīng)接近或達(dá)到你的預(yù)期那么就用它因?yàn)槟Ｐ透唵巍⒏?、更不容易過擬合。如果線性核效果不佳再切換到RBF核進(jìn)行精細(xì)調(diào)參。4.2 MATLAB中SVM的實(shí)現(xiàn)與關(guān)鍵步驟MATLAB的統(tǒng)計(jì)和機(jī)器學(xué)習(xí)工具箱Statistics and Machine Learning Toolbox提供了fitcsvm函數(shù)用于訓(xùn)練SVM分類器。對于多分類問題MATLAB默認(rèn)使用“一對一”策略。% 假設(shè)我們已經(jīng)準(zhǔn)備好了特征矩陣 X (n_samples x n_features) 和標(biāo)簽 y (n_samples x 1) % 1. 數(shù)據(jù)標(biāo)準(zhǔn)化 (非常重要) % SVM對特征的尺度敏感特別是使用RBF核時。必須標(biāo)準(zhǔn)化。 [X_train_scaled, mu, sigma] zscore(X_train); % 訓(xùn)練集標(biāo)準(zhǔn)化 X_test_scaled (X_test - mu) ./ sigma; % 測試集使用訓(xùn)練集的均值和標(biāo)準(zhǔn)差標(biāo)準(zhǔn)化 % 2. 訓(xùn)練一個簡單的線性SVM作為基線 linear_svm_model fitcsvm(X_train_scaled, y_train, ... KernelFunction, linear, ... Standardize, false, ... % 我們已經(jīng)手動標(biāo)準(zhǔn)化了所以這里設(shè)為false BoxConstraint, 1, ... % 懲罰系數(shù) C默認(rèn)是1 ClassNames, unique(y_train), ... Verbose, 1); % 3. 在測試集上預(yù)測 [y_pred_linear, score_linear] predict(linear_svm_model, X_test_scaled); accuracy_linear sum(y_pred_linear y_test) / numel(y_test); fprintf(線性SVM準(zhǔn)確率: %.2f%%\n, accuracy_linear * 100);如果線性SVM效果不理想我們轉(zhuǎn)向RBF核并需要進(jìn)行超參數(shù)調(diào)優(yōu)。最常用的方法是網(wǎng)格搜索交叉驗(yàn)證。4.3 超參數(shù)網(wǎng)格搜索與交叉驗(yàn)證C和γ的最佳值需要通過交叉驗(yàn)證來尋找。C控制誤分類的懲罰力度C越大模型越不允許犯錯決策邊界越復(fù)雜γ控制RBF核的寬度。% 定義參數(shù)網(wǎng)格 C_values [0.001, 0.01, 0.1, 1, 10, 100, 1000]; gamma_values [0.001, 0.01, 0.1, 1, 10, 100]; % 注意gamma 1/(2*sigma^2)有些定義不同 % 初始化最佳參數(shù)和準(zhǔn)確率 best_accuracy 0; best_C 1; best_gamma 1; % 創(chuàng)建交叉驗(yàn)證分區(qū) (例如5折) cv cvpartition(y_train, KFold, 5); % 網(wǎng)格搜索 for C C_values for gamma gamma_values fprintf(正在嘗試 C%.3f, gamma%.3f...\n, C, gamma); % 初始化本輪交叉驗(yàn)證的準(zhǔn)確率列表 fold_accuracies zeros(cv.NumTestSets, 1); for fold 1:cv.NumTestSets % 獲取當(dāng)前折的訓(xùn)練/驗(yàn)證索引 trainIdx cv.training(fold); testIdx cv.test(fold); % 劃分?jǐn)?shù)據(jù) X_fold_train X_train_scaled(trainIdx, :); y_fold_train y_train(trainIdx); X_fold_val X_train_scaled(testIdx, :); y_fold_val y_train(testIdx); % 訓(xùn)練SVM svm_model fitcsvm(X_fold_train, y_fold_train, ... KernelFunction, rbf, ... KernelScale, 1/sqrt(gamma), ... % MATLAB的KernelScale參數(shù)是sigmagamma1/(2*sigma^2) BoxConstraint, C, ... Standardize, false); % 驗(yàn)證 y_fold_pred predict(svm_model, X_fold_val); fold_accuracies(fold) sum(y_fold_pred y_fold_val) / numel(y_fold_val); end % 計(jì)算平均交叉驗(yàn)證準(zhǔn)確率 mean_cv_accuracy mean(fold_accuracies); fprintf( 平均CV準(zhǔn)確率: %.4f\n, mean_cv_accuracy); % 更新最佳參數(shù) if mean_cv_accuracy best_accuracy best_accuracy mean_cv_accuracy; best_C C; best_gamma gamma; end end end fprintf(\n最佳參數(shù): C %.3f, gamma %.3f, 最佳CV準(zhǔn)確率: %.4f\n, best_C, best_gamma, best_accuracy); % 使用最佳參數(shù)在整個訓(xùn)練集上訓(xùn)練最終模型 final_svm_model fitcsvm(X_train_scaled, y_train, ... KernelFunction, rbf, ... KernelScale, 1/sqrt(best_gamma), ... BoxConstraint, best_C, ... Standardize, false); % 在獨(dú)立測試集上評估最終模型 [y_pred_final, score_final] predict(final_svm_model, X_test_scaled); test_accuracy sum(y_pred_final y_test) / numel(y_test); fprintf(最終模型在測試集上的準(zhǔn)確率: %.2f%%\n, test_accuracy * 100);踩坑提醒數(shù)據(jù)泄露標(biāo)準(zhǔn)化zscore必須在劃分訓(xùn)練集和測試集之后分別進(jìn)行。必須用訓(xùn)練集的均值和標(biāo)準(zhǔn)差去標(biāo)準(zhǔn)化測試集絕對不能用整個數(shù)據(jù)集計(jì)算均值和標(biāo)準(zhǔn)差然后統(tǒng)一標(biāo)準(zhǔn)化否則測試集信息就“泄露”到訓(xùn)練過程了會嚴(yán)重高估模型性能。類別不平衡如果四個情感類別樣本數(shù)差異很大需要在fitcsvm中設(shè)置Prior參數(shù)為empirical讓模型根據(jù)訓(xùn)練數(shù)據(jù)自動調(diào)整先驗(yàn)概率或者設(shè)置Weights參數(shù)給少數(shù)類更高的權(quán)重。網(wǎng)格搜索的代價(jià)參數(shù)網(wǎng)格C_values和gamma_values的乘積乘以交叉驗(yàn)證折數(shù)就是需要訓(xùn)練模型的次數(shù)。如果數(shù)據(jù)量大、特征多這個過程會非常耗時?？梢詮拇志W(wǎng)格開始如[0.01, 0.1, 1, 10, 100]找到大致范圍后再在附近進(jìn)行精細(xì)搜索。也可以使用更高效的優(yōu)化算法如貝葉斯優(yōu)化bayesopt函數(shù)。MATLAB的KernelScale注意fitcsvm的KernelScale參數(shù)對應(yīng)的是RBF核公式中的sigma而通常我們說的gamma參數(shù)是1/(2*sigma^2)。所以代碼中我們用1/sqrt(gamma)來轉(zhuǎn)換。務(wù)必理清你所用工具包對核參數(shù)的定義。5. 模型評估、結(jié)果分析與可視化得到一個準(zhǔn)確率數(shù)字只是開始更重要的是理解模型在哪里做得好在哪里做得不好。5.1 超越準(zhǔn)確率全面的評估指標(biāo)對于四分類問題混淆矩陣Confusion Matrix是最直觀的工具。% 計(jì)算混淆矩陣 C confusionmat(y_test, y_pred_final); % 使用內(nèi)置函數(shù)繪圖更美觀 figure; confusionchart(C, {LALV, LAHV, HALV, HAHV}); % 按你的標(biāo)簽順序 title(SVM分類器混淆矩陣);從混淆矩陣中我們可以計(jì)算每個類別的精確率Precision、召回率Recall和F1分?jǐn)?shù)。% 計(jì)算每類的性能指標(biāo) num_classes 4; precision zeros(num_classes, 1); recall zeros(num_classes, 1); f1_score zeros(num_classes, 1); for i 1:num_classes TP C(i, i); FP sum(C(:, i)) - TP; FN sum(C(i, :)) - TP; precision(i) TP / (TP FP eps); recall(i) TP / (TP FN eps); f1_score(i) 2 * (precision(i) * recall(i)) / (precision(i) recall(i) eps); end % 顯示結(jié)果 class_names {LALV, LAHV, HALV, HAHV}; for i 1:num_classes fprintf(類別 %s: 精確率%.3f, 召回率%.3f, F1%.3f\n, ... class_names{i}, precision(i), recall(i), f1_score(i)); end % 計(jì)算宏平均F1 macro_f1 mean(f1_score); fprintf(宏平均F1分?jǐn)?shù): %.3f\n, macro_f1);為什么F1分?jǐn)?shù)很重要在類別不平衡的情況下準(zhǔn)確率可能會被大類別主導(dǎo)而虛高。比如90%的樣本都是“平靜”模型只要全部預(yù)測“平靜”就能得到90%的準(zhǔn)確率但這毫無意義。F1分?jǐn)?shù)是精確率和召回率的調(diào)和平均能更好地衡量模型對每個類別的識別能力宏平均F1對所有類別一視同仁。5.2 結(jié)果可視化與解讀特征重要性分析對于線性SVM可以通過檢查權(quán)重向量w的絕對值大小來估計(jì)特征的重要性。權(quán)重絕對值越大說明該特征對決策的影響越大。if strcmp(final_svm_model.KernelParameters.Function, linear) w final_svm_model.Beta; % 權(quán)重系數(shù) [sorted_w, idx] sort(abs(w), descend); top_n 20; fprintf(重要性排名前%d的特征索引:\n, top_n); disp(idx(1:top_n)); % 你可以將這些索引映射回具體的通道和特征類型如“F3通道的Alpha功率” end決策邊界可視化降維后對于高維數(shù)據(jù)我們可以使用t-SNE或PCA將其降到2維或3維然后畫出樣本點(diǎn)和SVM的決策邊界需要訓(xùn)練一個在降維空間上的新SVM這有助于直觀感受數(shù)據(jù)的可分性。% 使用PCA降維到2維 [coeff, score, latent] pca(X_train_scaled); X_train_2d score(:, 1:2); % 在2維數(shù)據(jù)上重新訓(xùn)練一個SVM僅用于可視化 svm_2d fitcsvm(X_train_2d, y_train, KernelFunction, rbf, Standardize, false); % 創(chuàng)建網(wǎng)格用于繪制決策區(qū)域 x1range min(X_train_2d(:,1))-1 : 0.1 : max(X_train_2d(:,1))1; x2range min(X_train_2d(:,2))-1 : 0.1 : max(X_train_2d(:,2))1; [xx1, xx2] meshgrid(x1range, x2range); XGrid [xx1(:), xx2(:)]; % 預(yù)測網(wǎng)格點(diǎn)的類別 [~, scores] predict(svm_2d, XGrid); % 繪制決策區(qū)域和樣本點(diǎn) figure; h gscatter(X_train_2d(:,1), X_train_2d(:,2), y_train, brgk, o*^v); hold on; contour(xx1, xx2, reshape(scores(:,4), size(xx1)), [0 0], k-, LineWidth, 2); % 繪制第4類HAHV的決策邊界 xlabel(PCA主成分1); ylabel(PCA主成分2); title(PCA降維后的數(shù)據(jù)分布與SVM決策邊界示例); legend(LALV,LAHV,HALV,HAHV, Decision Boundary); hold off;5.3 與基線模型和現(xiàn)有研究的對比為了評估你的SVM模型是否有效需要設(shè)立基線進(jìn)行比較。隨機(jī)猜測基線四分類問題的隨機(jī)猜測準(zhǔn)確率是25%。多數(shù)類基線預(yù)測為訓(xùn)練集中樣本數(shù)最多的類別計(jì)算其在測試集上的準(zhǔn)確率。簡單分類器基線如最近鄰分類器KNNK1或決策樹最大深度5。% 多數(shù)類基線 majority_class mode(y_train); y_pred_majority repmat(majority_class, size(y_test)); accuracy_majority sum(y_pred_majority y_test) / numel(y_test); fprintf(多數(shù)類基線準(zhǔn)確率: %.2f%%\n, accuracy_majority * 100); % 1-NN基線 knn_model fitcknn(X_train_scaled, y_train, NumNeighbors, 1, Standardize, false); y_pred_knn predict(knn_model, X_test_scaled); accuracy_knn sum(y_pred_knn y_test) / numel(y_test); fprintf(1-NN基線準(zhǔn)確率: %.2f%%\n, accuracy_knn * 100);如果你的SVM模型顯著優(yōu)于這些基線例如準(zhǔn)確率超過40%F1分?jǐn)?shù)超過0.4并且與相關(guān)論文中報(bào)告的基于DEAP的SVM性能通常在55%-70%之間取決于特征、通道選擇和分類任務(wù)處于同一量級那么你的工作就是有效的。6. 項(xiàng)目總結(jié)、局限性與進(jìn)階方向走完整個流程從數(shù)據(jù)預(yù)處理、特征提取到模型訓(xùn)練評估一個完整的基于SVM的腦電情感識別項(xiàng)目就完成了?；仡櫿麄€過程有幾個關(guān)鍵點(diǎn)值得再次強(qiáng)調(diào)預(yù)處理是基石濾波和去噪的質(zhì)量直接決定了特征提取的上限。務(wù)必根據(jù)腦電信號的特性選擇合適的濾波參數(shù)。特征工程是核心時頻域特征特別是頻帶功率、微分熵對情感識別非常有效。嘗試組合不同類型的特征時域、頻域、非線性特征有時能提升性能但要注意維度控制。模型調(diào)優(yōu)是藝術(shù)SVM的性能對超參數(shù)和特征尺度敏感。網(wǎng)格搜索交叉驗(yàn)證是標(biāo)準(zhǔn)做法但計(jì)算成本高。務(wù)必做好數(shù)據(jù)標(biāo)準(zhǔn)化防止泄露。評估要全面不要只看準(zhǔn)確率混淆矩陣和F1分?jǐn)?shù)能揭示模型在各類別上的具體表現(xiàn)尤其是在數(shù)據(jù)不平衡時。本項(xiàng)目的局限性特征手工設(shè)計(jì)我們提取的時頻特征是手工設(shè)計(jì)的依賴于領(lǐng)域知識。這可能不是最優(yōu)的表示。通道選擇我們使用了所有通道但有些通道可能對情感識別貢獻(xiàn)很小甚至引入噪聲。通道選擇是一個重要的優(yōu)化方向。個體差異腦電信號存在巨大的個體差異。在一個被試上訓(xùn)練好的模型在另一個被試上可能效果很差。這就是“被試獨(dú)立”評估的挑戰(zhàn)。SVM的局限性SVM雖然強(qiáng)大但對于非常復(fù)雜的非線性關(guān)系其表現(xiàn)可能不如深度神經(jīng)網(wǎng)絡(luò)。進(jìn)階方向深度學(xué)習(xí)端到端學(xué)習(xí)使用卷積神經(jīng)網(wǎng)絡(luò)直接從原始腦電信號或時頻圖中學(xué)習(xí)特征避免手工特征工程。例如EEGNet、DeepConvNet等網(wǎng)絡(luò)結(jié)構(gòu)在此任務(wù)上表現(xiàn)出色?？绫辉噷W(xué)習(xí)與域自適應(yīng)研究如何讓模型更好地泛化到新被試?？梢允褂糜蜃赃m應(yīng)、遷移學(xué)習(xí)或元學(xué)習(xí)等方法。多模態(tài)融合DEAP數(shù)據(jù)集除了EEG還有外周生理信號GSR, RESP, EMG等。融合多模態(tài)信息可以顯著提升情感識別的魯棒性和準(zhǔn)確性。時序建模情感是動態(tài)變化的使用循環(huán)神經(jīng)網(wǎng)絡(luò)或Transformer對腦電信號的時間序列進(jìn)行建?？梢圆蹲狡鋭討B(tài)演化模式。更精細(xì)的情感模型除了四分類也可以嘗試回歸任務(wù)直接預(yù)測效價(jià)和喚醒度的連續(xù)值或者探索更復(fù)雜的情感模型。這個項(xiàng)目提供了一個堅(jiān)實(shí)的起點(diǎn)。你可以基于這個代碼框架嘗試不同的特征、不同的分類器、不同的預(yù)處理方法逐步深入腦電情感識別這個有趣且充滿挑戰(zhàn)的領(lǐng)域。代碼的模塊化設(shè)計(jì)使得替換特征提取模塊或分類器模塊變得非常容易祝你實(shí)驗(yàn)順利。本文還有配套的精品資源點(diǎn)擊獲取 基于DEAP數(shù)據(jù)集與SVM的腦電情感識別：從預(yù)處理到模型調(diào)優(yōu)全流程實(shí)踐
簡介本資源是一套面向腦機(jī)接口與情感計(jì)算方向研究者的完整論文代碼實(shí)現(xiàn)方案聚焦DEAP數(shù)據(jù)集上的四分類情緒識別任務(wù)效價(jià)-喚醒二維象限劃分適用于具備MATLAB基礎(chǔ)與機(jī)器學(xué)習(xí)入門知識的研究生及科研人員。資源包含16個文件以10個核心MATLAB腳本如頻帶濾波、DWT特征提取、SVM訓(xùn)練等、3個Jupyter Notebook含KNN對比實(shí)驗(yàn)、1篇PDF論文、1份README說明及輔助文本文件為主總大小僅488KB結(jié)構(gòu)緊湊、模塊清晰便于復(fù)現(xiàn)實(shí)驗(yàn)流程與特征工程細(xì)節(jié)。已有3636人學(xué)習(xí)下載讀者可直接運(yùn)行代碼完成從原始EEG信號預(yù)處理、α/β/θ頻帶分離、箱線圖特征篩選、統(tǒng)計(jì)特征偏度、峰度、波熵構(gòu)建到SVM四分類建模的全流程并獲得92.36%測試準(zhǔn)確率結(jié)果顯著優(yōu)于現(xiàn)有主流方法附帶可視化腳本支持特征分布與分類效果分析。1. 項(xiàng)目背景與核心目標(biāo)最近在整理過往的研究項(xiàng)目翻到了一個基于DEAP數(shù)據(jù)集做情感識別的老項(xiàng)目感覺挺有意思的。這個項(xiàng)目說白了就是通過分析人的腦電信號來判斷他/她當(dāng)前是高興、悲傷、平靜還是興奮。聽起來有點(diǎn)科幻但其實(shí)背后的技術(shù)路徑已經(jīng)相當(dāng)成熟。當(dāng)時我用的是經(jīng)典的機(jī)器學(xué)習(xí)方法核心流程就是從原始的腦電信號里提取特征然后用支持向量機(jī)SVM去訓(xùn)練一個分類器。整個過程在MATLAB里實(shí)現(xiàn)從數(shù)據(jù)預(yù)處理、特征工程到模型訓(xùn)練和評估一條龍搞定。今天這篇文章我就把這個項(xiàng)目的完整實(shí)現(xiàn)思路、關(guān)鍵代碼和踩過的坑系統(tǒng)地梳理一遍希望能給對腦機(jī)接口、情感計(jì)算或者信號處理感興趣的朋友提供一個清晰的參考模板。無論你是剛?cè)腴T的研究生還是想快速復(fù)現(xiàn)一個baseline的工程師這篇內(nèi)容應(yīng)該都能幫到你。DEAP數(shù)據(jù)集在情感計(jì)算領(lǐng)域算是“明星數(shù)據(jù)集”了它包含了32位被試者在觀看音樂視頻時的腦電EEG和外圍生理信號如肌電、皮電等并且對每段視頻都進(jìn)行了效價(jià)Valence、喚醒度Arousal、支配度Dominance和喜愛度Liking的評分。我們通常取效價(jià)和喚醒度這兩個維度將它們劃分為高/低兩類然后組合成四個象限就對應(yīng)了四種情感狀態(tài)高喚醒高效價(jià)HAHV可理解為興奮/高興、高喚醒低效價(jià)HALV可理解為焦慮/憤怒、低喚醒高效價(jià)LAHV可理解為平靜/放松、低喚醒低效價(jià)LALV可理解為悲傷/低落。我們的任務(wù)就是讓機(jī)器通過腦電信號認(rèn)出這四種狀態(tài)。2. DEAP數(shù)據(jù)集預(yù)處理與關(guān)鍵理解拿到DEAP數(shù)據(jù)集后第一步不是急著跑代碼而是先理解它的數(shù)據(jù)結(jié)構(gòu)。DEAP數(shù)據(jù)以MAT文件.mat格式提供通常包含兩個關(guān)鍵變量data和labels。data是一個4維數(shù)組維度通常是參與者 × 視頻片段 × 通道 × 數(shù)據(jù)點(diǎn)例如 32 x 40 x 32 x 8064。這里8064對應(yīng)的是63秒的視頻片段采樣頻率是128Hz所以63*1288064個數(shù)據(jù)點(diǎn)。labels則包含了每個視頻片段對應(yīng)的效價(jià)、喚醒度等評分。預(yù)處理的目標(biāo)是把原始的、充滿噪聲的腦電信號變成干凈、可用于特征提取的信號。這個過程至關(guān)重要直接決定了后續(xù)特征的質(zhì)量和模型的性能上限。2.1 數(shù)據(jù)加載與維度重塑首先我們需要把數(shù)據(jù)加載進(jìn)來并轉(zhuǎn)換成更易于處理的格式。原始的四維數(shù)組在操作時不太方便我們通常會把數(shù)據(jù)重塑成二維矩陣樣本 × 特征但這里先按被試和試驗(yàn)進(jìn)行分割。% 假設(shè)數(shù)據(jù)文件為 DEAP_data.mat load(DEAP_data.mat); % 這會加載變量 data 和 labels % 查看數(shù)據(jù)維度 [num_participants, num_trials, num_channels, num_samples] size(data); disp([數(shù)據(jù)維度: , num2str(num_participants), 參與者, , ... num2str(num_trials), 試驗(yàn), , num2str(num_channels), 通道, , ... num2str(num_samples), 采樣點(diǎn)]); % 提取標(biāo)簽這里以效價(jià)和喚醒度為例假設(shè)它們是labels的前兩列 valence_labels labels(:, 1); % 效價(jià)維度 arousal_labels labels(:, 2); % 喚醒度維度2.2 腦電信號預(yù)處理流程詳解腦電信號非常微弱容易受到眼電EOG、肌電EMG、工頻干擾50Hz等噪聲污染。標(biāo)準(zhǔn)的預(yù)處理流程包括以下幾步重參考原始腦電記錄通常是單極導(dǎo)聯(lián)需要一個參考點(diǎn)。DEAP數(shù)據(jù)已使用平均參考這一步通?？梢蕴^但了解其概念很重要。平均參考是假設(shè)所有電極電位的平均值為零將每個通道的信號減去所有通道的平均值有助于減少共同噪聲。帶通濾波保留有用的腦電節(jié)律成分。一般認(rèn)為與情感相關(guān)的腦電成分分布在Delta(1-4Hz), Theta(4-8Hz), Alpha(8-13Hz), Beta(13-30Hz), Gamma(30-45Hz)等頻段。我們通常用一個較寬的帶通濾波器比如0.5Hz到45Hz以保留大部分有效信息同時濾除極低頻的漂移和高頻的肌電噪聲。陷波濾波去除工頻干擾50Hz或60Hz取決于地區(qū)。在中國是50Hz。分段與基線校正DEAP數(shù)據(jù)已經(jīng)是觀看視頻期間的63秒連續(xù)數(shù)據(jù)。有時我們需要進(jìn)行更短時間窗的分析如分析事件相關(guān)電位ERP但在這個項(xiàng)目中我們直接使用整個63秒的片段作為一個樣本?；€校正通常指減去一段“靜息”或“刺激前”時期的平均幅值以消除個體直流偏移。DEAP數(shù)據(jù)在提供前可能已做過處理但為了嚴(yán)謹(jǐn)我們可以計(jì)算每個試驗(yàn)最開始一小段如前1秒數(shù)據(jù)的均值然后從整個試驗(yàn)數(shù)據(jù)中減去它。下面是用MATLAB實(shí)現(xiàn)濾波的示例代碼。我強(qiáng)烈推薦使用EEGLAB或FieldTrip這些專業(yè)的腦電處理工具箱它們功能強(qiáng)大且穩(wěn)定。但為了理解原理這里用MATLAB內(nèi)置函數(shù)演示% 定義采樣頻率和濾波器參數(shù) Fs 128; % DEAP采樣率 low_cutoff 0.5; % 高通截止頻率 (Hz) high_cutoff 45.0; % 低通截止頻率 (Hz) notch_freq 50; % 陷波頻率 (Hz) % 設(shè)計(jì)一個帶通濾波器例如使用巴特沃斯濾波器 order 4; % 濾波器階數(shù) [b_band, a_band] butter(order, [low_cutoff, high_cutoff]/(Fs/2), bandpass); % 設(shè)計(jì)一個陷波濾波器 wo notch_freq/(Fs/2); % 歸一化頻率 bw wo/35; % 帶寬 [b_notch, a_notch] iirnotch(wo, bw); % 假設(shè)我們處理一個被試一個試驗(yàn)的一個通道數(shù)據(jù) trial_eeg (1 x 8064) % 1. 帶通濾波 filtered_eeg filtfilt(b_band, a_band, trial_eeg); % 使用filtfilt實(shí)現(xiàn)零相位濾波 % 2. 陷波濾波 filtered_eeg filtfilt(b_notch, a_notch, filtered_eeg); % 注意對多通道數(shù)據(jù)需要循環(huán)處理每個通道。 % 更高效的做法是使用eeglab的pop_eegfiltnew函數(shù) % EEG pop_eegfiltnew(EEG, low_cutoff, high_cutoff); % EEG pop_eegfiltnew(EEG, notch_freq-1, notch_freq1, [], 1); % 陷波注意filtfilt函數(shù)執(zhí)行的是零相位濾波它通過對數(shù)據(jù)正向和反向各濾波一次消除了濾波器帶來的相位延遲這對于后續(xù)的時頻分析非常重要。但它的計(jì)算量是普通濾波的兩倍。2.3 標(biāo)簽生成從維度到四分類DEAP的原始標(biāo)簽是連續(xù)的效價(jià)和喚醒度評分1-9分。我們需要將其離散化為高/低兩類通常以中值5為界。% 假設(shè)valence_labels和arousal_labels是來自所有試驗(yàn)的向量 valence_threshold 5; arousal_threshold 5; % 生成二值標(biāo)簽 valence_binary valence_labels valence_threshold; % 1高效價(jià)0低效價(jià) arousal_binary arousal_labels arousal_threshold; % 1高喚醒0低喚醒 % 組合成四分類標(biāo)簽 (1: LALV, 2: LAHV, 3: HALV, 4: HAHV) % 注意這里的映射關(guān)系可以根據(jù)你的情感模型定義調(diào)整 four_class_labels zeros(size(valence_labels)); four_class_labels(~arousal_binary ~valence_binary) 1; % 低喚醒低效價(jià) four_class_labels(~arousal_binary valence_binary) 2; % 低喚醒高效價(jià) four_class_labels(arousal_binary ~valence_binary) 3; % 高喚醒低效價(jià) four_class_labels(arousal_binary valence_binary) 4; % 高喚醒高效價(jià) % 檢查類別分布 tabulate(four_class_labels)類別不平衡是情感識別中常見的問題。DEAP數(shù)據(jù)集的四個類別分布可能不是完全均勻的。如果差異較大在后續(xù)訓(xùn)練模型時需要考慮采用加權(quán)SVM、過采樣如SMOTE或欠采樣等策略。3. 時頻域特征提取從信號到信息特征提取是情感識別的核心。腦電信號在時域上看就是一條隨時間變化的電壓曲線信息隱藏很深。時頻分析能告訴我們信號的能量在不同頻率成分上是如何隨時間變化的這比單純的時域或頻域特征更能捕捉情感的動態(tài)特性。3.1 為何選擇時頻域特征情感變化不是瞬間切換的而是一個動態(tài)過程。例如從平靜到興奮大腦中Alpha波8-13Hz的能量可能會降低去同步化而Beta波13-30Hz的能量可能會增加。單純的頻域特征如整個片段的功率譜會丟失時間信息而單純的時域特征如均值、方差對頻率不敏感。時頻分析如小波變換、短時傅里葉變換提供了“時間-頻率”二維平面上的能量分布讓我們能夠提取諸如“在視頻播放的第20秒到30秒前額葉Theta波能量上升”這樣的動態(tài)特征這對情感識別至關(guān)重要。3.2 基于短時傅里葉變換STFT的特征提取實(shí)踐STFT是理解時頻分析最直觀的工具。它的思想很簡單把長信號分成許多短的、重疊的片段對每個片段做傅里葉變換從而得到每個時間點(diǎn)附近的頻率成分。% 假設(shè) filtered_eeg 是一個預(yù)處理后的單通道腦電信號 (1 x N) signal filtered_eeg; N length(signal); Fs 128; % 設(shè)置STFT參數(shù) window_length 256; % 窗長對應(yīng)2秒 (256/1282)。窗越長頻率分辨率越高時間分辨率越低。 noverlap 128; % 重疊點(diǎn)數(shù)通常為窗長的一半。重疊越多時頻圖越平滑。 nfft 512; % FFT點(diǎn)數(shù)通常大于等于窗長。決定頻率軸的精細(xì)度。 % 計(jì)算STFT [S, F, T] spectrogram(signal, window_length, noverlap, nfft, Fs); % S是復(fù)數(shù)矩陣維度為 (nfft/21) x 時間窗數(shù) % F是頻率向量 % T是時間向量每個窗的中心時間 % 計(jì)算功率譜密度 (PSD) P abs(S).^2; % 現(xiàn)在P就是一個時頻能量矩陣。我們可以從中提取各種特征。從時頻矩陣P中我們可以提取多種特征特定頻帶功率這是最常用的特征。先定義頻帶范圍然后對P矩陣中對應(yīng)頻率區(qū)間的功率進(jìn)行求和或平均。% 定義經(jīng)典頻帶邊界 (Hz) band_defs { Delta, 1, 4; Theta, 4, 8; Alpha, 8, 13; Beta, 13, 30; Gamma, 30, 45; }; % 初始化特征向量 band_powers []; for i 1:size(band_defs, 1) band_name band_defs{i, 1}; f_low band_defs{i, 2}; f_high band_defs{i, 3}; % 找到頻率向量F中在頻帶內(nèi)的索引 freq_indices (F f_low) (F f_high); % 計(jì)算該頻帶在所有時間窗上的總功率然后對所有時間窗取平均 % 也可以考慮取對數(shù)(log)使分布更接近正態(tài) band_power_total sum(P(freq_indices, :), 1); % 對頻率維求和得到每個時間窗的該頻帶功率 band_power_mean mean(band_power_total); % 對所有時間窗求平均得到一個標(biāo)量特征 % band_power_log log(band_power_mean eps); % 取對數(shù)加eps防止為0 band_powers [band_powers, band_power_mean]; end % 現(xiàn)在 band_powers 是一個1x5的向量包含5個頻帶的平均功率時頻統(tǒng)計(jì)特征我們可以把時頻矩陣P看作一幅圖像計(jì)算其統(tǒng)計(jì)特性。平均功率mean(P(:))功率標(biāo)準(zhǔn)差std(P(:))功率偏度/峰度skewness(P(:)),kurtosis(P(:))描述功率分布的形態(tài)。譜熵將每個時間點(diǎn)的功率譜歸一化為概率分布計(jì)算其香農(nóng)熵反映頻譜的復(fù)雜度或平坦度。% 計(jì)算譜熵 (Spectral Entropy) % 對每個時間窗的功率譜進(jìn)行歸一化 P_norm zeros(size(P)); for t 1:size(P, 2) p_slice P(:, t); p_slice p_slice / sum(p_slice); % 歸一化為概率分布 P_norm(:, t) p_slice; end % 計(jì)算每個時間窗的譜熵 spectral_entropy_per_window -sum(P_norm .* log2(P_norm eps), 1); % 取所有時間窗譜熵的平均值作為特征 mean_spectral_entropy mean(spectral_entropy_per_window);微分熵特征在腦電分析中微分熵Differential Entropy, DE被證明是有效的特征特別是在基于頻帶的特征提取中。對于一段服從高斯分布N(μ, σ^2)的信號其微分熵計(jì)算公式為0.5 * log(2πeσ^2)。在實(shí)際操作中我們通常假設(shè)每個頻帶內(nèi)的信號近似服從高斯分布因此可以用該頻帶功率的對數(shù)變換來近似微分熵。許多研究直接使用log(頻帶功率)作為特征效果很好。% 基于之前計(jì)算的 band_power_total (每個時間窗的頻帶功率) de_features []; for i 1:size(band_defs, 1) % 假設(shè) band_power_total 是之前計(jì)算好的第i個頻帶的功率時間序列 % 這里需要重新計(jì)算每個頻帶每個時間窗的功率 freq_indices (F band_defs{i, 2}) (F band_defs{i, 3}); band_power_per_window sum(P(freq_indices, :), 1); % 1 x 時間窗數(shù) % 計(jì)算該頻帶功率的微分熵近似為對數(shù)功率 % 可以對每個時間窗取對數(shù)后平均也可以先平均再取對數(shù)效果略有不同 de_per_window log(band_power_per_window eps); mean_de mean(de_per_window); de_features [de_features, mean_de]; end3.3 多通道特征融合與特征向量構(gòu)建上面我們只針對一個通道進(jìn)行了特征提取。DEAP有32個通道或64個我們需要對所有通道重復(fù)上述過程然后將所有特征拼接成一個長向量代表一個試驗(yàn)樣本。% 假設(shè)我們有一個函數(shù) extract_features_from_channel(eeg_signal, Fs) % 它輸入一個通道的信號輸出一個特征向量 (例如 5個頻帶功率 1個譜熵 6維) num_channels 32; features_per_trial []; for ch 1:num_channels single_channel_signal squeeze(data(participant_idx, trial_idx, ch, :)); % 獲取指定被試、試驗(yàn)、通道的數(shù)據(jù) single_channel_signal preprocess_eeg(single_channel_signal, Fs); % 預(yù)處理函數(shù) ch_features extract_features_from_channel(single_channel_signal, Fs); features_per_trial [features_per_trial, ch_features]; end % 最終features_per_trial 是一個 1 x (6 * 32) 1 x 192 維的特征向量。這就是一個樣本的特征向量。對數(shù)據(jù)集中的所有試驗(yàn)和所有被試重復(fù)此操作注意要按被試劃分訓(xùn)練測試集避免數(shù)據(jù)泄露我們就能得到特征矩陣X(樣本數(shù) x 特征數(shù)) 和標(biāo)簽向量y。實(shí)操心得特征維度會爆炸式增長通道數(shù) x 每通道特征數(shù)。32個通道每個通道提取6個特征就是192維。如果使用更復(fù)雜的時頻特征如Hjorth參數(shù)、分形維數(shù)等維度會更高。高維特征容易導(dǎo)致“維數(shù)災(zāi)難”并且包含大量冗余信息。因此特征選擇是必不可少的后續(xù)步驟。我們可以使用方差閾值、相關(guān)系數(shù)、遞歸特征消除RFE或基于模型如L1正則化的方法來篩選出最具有判別力的特征子集。在項(xiàng)目初期為了快速驗(yàn)證流程可以先用所有特征但最終一定要做特征選擇來優(yōu)化模型。4. 支持向量機(jī)分類器的構(gòu)建與調(diào)優(yōu)特征準(zhǔn)備好之后就進(jìn)入了建模階段。支持向量機(jī)SVM在小樣本、高維度的分類問題上表現(xiàn)優(yōu)異并且可以通過核函數(shù)處理非線性問題非常適合作為腦電情感識別的基線分類器。4.1 SVM原理簡述與核函數(shù)選擇SVM的核心思想是尋找一個最優(yōu)超平面使得兩類樣本到這個超平面的“間隔”最大化。對于線性不可分的數(shù)據(jù)SVM通過“核技巧”將數(shù)據(jù)映射到高維空間使其在高維空間中線性可分。常用的核函數(shù)有線性核K(x, z) x^T * z。參數(shù)少速度快可解釋性強(qiáng)。如果特征已經(jīng)足夠好或者特征維度很高線性核往往就能取得不錯的效果并且不容易過擬合。徑向基函數(shù)核K(x, z) exp(-γ * ||x - z||^2)。也叫高斯核是最常用的非線性核。它能夠處理非常復(fù)雜的非線性關(guān)系但需要調(diào)節(jié)兩個超參數(shù)懲罰系數(shù)C和核系數(shù)γ。γ定義了單個樣本的影響范圍γ越大影響范圍越小決策邊界越復(fù)雜容易過擬合γ越小決策邊界越平滑容易欠擬合。對于腦電情感識別由于特征與情感狀態(tài)之間的關(guān)系可能是非線性的RBF核通常是首選。但我的經(jīng)驗(yàn)是先用線性核試一下。如果線性核的準(zhǔn)確率已經(jīng)接近或達(dá)到你的預(yù)期那么就用它因?yàn)槟Ｐ透唵巍⒏?、更不容易過擬合。如果線性核效果不佳再切換到RBF核進(jìn)行精細(xì)調(diào)參。4.2 MATLAB中SVM的實(shí)現(xiàn)與關(guān)鍵步驟MATLAB的統(tǒng)計(jì)和機(jī)器學(xué)習(xí)工具箱Statistics and Machine Learning Toolbox提供了fitcsvm函數(shù)用于訓(xùn)練SVM分類器。對于多分類問題MATLAB默認(rèn)使用“一對一”策略。% 假設(shè)我們已經(jīng)準(zhǔn)備好了特征矩陣 X (n_samples x n_features) 和標(biāo)簽 y (n_samples x 1) % 1. 數(shù)據(jù)標(biāo)準(zhǔn)化 (非常重要) % SVM對特征的尺度敏感特別是使用RBF核時。必須標(biāo)準(zhǔn)化。 [X_train_scaled, mu, sigma] zscore(X_train); % 訓(xùn)練集標(biāo)準(zhǔn)化 X_test_scaled (X_test - mu) ./ sigma; % 測試集使用訓(xùn)練集的均值和標(biāo)準(zhǔn)差標(biāo)準(zhǔn)化 % 2. 訓(xùn)練一個簡單的線性SVM作為基線 linear_svm_model fitcsvm(X_train_scaled, y_train, ... KernelFunction, linear, ... Standardize, false, ... % 我們已經(jīng)手動標(biāo)準(zhǔn)化了所以這里設(shè)為false BoxConstraint, 1, ... % 懲罰系數(shù) C默認(rèn)是1 ClassNames, unique(y_train), ... Verbose, 1); % 3. 在測試集上預(yù)測 [y_pred_linear, score_linear] predict(linear_svm_model, X_test_scaled); accuracy_linear sum(y_pred_linear y_test) / numel(y_test); fprintf(線性SVM準(zhǔn)確率: %.2f%%\n, accuracy_linear * 100);如果線性SVM效果不理想我們轉(zhuǎn)向RBF核并需要進(jìn)行超參數(shù)調(diào)優(yōu)。最常用的方法是網(wǎng)格搜索交叉驗(yàn)證。4.3 超參數(shù)網(wǎng)格搜索與交叉驗(yàn)證C和γ的最佳值需要通過交叉驗(yàn)證來尋找。C控制誤分類的懲罰力度C越大模型越不允許犯錯決策邊界越復(fù)雜γ控制RBF核的寬度。% 定義參數(shù)網(wǎng)格 C_values [0.001, 0.01, 0.1, 1, 10, 100, 1000]; gamma_values [0.001, 0.01, 0.1, 1, 10, 100]; % 注意gamma 1/(2*sigma^2)有些定義不同 % 初始化最佳參數(shù)和準(zhǔn)確率 best_accuracy 0; best_C 1; best_gamma 1; % 創(chuàng)建交叉驗(yàn)證分區(qū) (例如5折) cv cvpartition(y_train, KFold, 5); % 網(wǎng)格搜索 for C C_values for gamma gamma_values fprintf(正在嘗試 C%.3f, gamma%.3f...\n, C, gamma); % 初始化本輪交叉驗(yàn)證的準(zhǔn)確率列表 fold_accuracies zeros(cv.NumTestSets, 1); for fold 1:cv.NumTestSets % 獲取當(dāng)前折的訓(xùn)練/驗(yàn)證索引 trainIdx cv.training(fold); testIdx cv.test(fold); % 劃分?jǐn)?shù)據(jù) X_fold_train X_train_scaled(trainIdx, :); y_fold_train y_train(trainIdx); X_fold_val X_train_scaled(testIdx, :); y_fold_val y_train(testIdx); % 訓(xùn)練SVM svm_model fitcsvm(X_fold_train, y_fold_train, ... KernelFunction, rbf, ... KernelScale, 1/sqrt(gamma), ... % MATLAB的KernelScale參數(shù)是sigmagamma1/(2*sigma^2) BoxConstraint, C, ... Standardize, false); % 驗(yàn)證 y_fold_pred predict(svm_model, X_fold_val); fold_accuracies(fold) sum(y_fold_pred y_fold_val) / numel(y_fold_val); end % 計(jì)算平均交叉驗(yàn)證準(zhǔn)確率 mean_cv_accuracy mean(fold_accuracies); fprintf( 平均CV準(zhǔn)確率: %.4f\n, mean_cv_accuracy); % 更新最佳參數(shù) if mean_cv_accuracy best_accuracy best_accuracy mean_cv_accuracy; best_C C; best_gamma gamma; end end end fprintf(\n最佳參數(shù): C %.3f, gamma %.3f, 最佳CV準(zhǔn)確率: %.4f\n, best_C, best_gamma, best_accuracy); % 使用最佳參數(shù)在整個訓(xùn)練集上訓(xùn)練最終模型 final_svm_model fitcsvm(X_train_scaled, y_train, ... KernelFunction, rbf, ... KernelScale, 1/sqrt(best_gamma), ... BoxConstraint, best_C, ... Standardize, false); % 在獨(dú)立測試集上評估最終模型 [y_pred_final, score_final] predict(final_svm_model, X_test_scaled); test_accuracy sum(y_pred_final y_test) / numel(y_test); fprintf(最終模型在測試集上的準(zhǔn)確率: %.2f%%\n, test_accuracy * 100);踩坑提醒數(shù)據(jù)泄露標(biāo)準(zhǔn)化zscore必須在劃分訓(xùn)練集和測試集之后分別進(jìn)行。必須用訓(xùn)練集的均值和標(biāo)準(zhǔn)差去標(biāo)準(zhǔn)化測試集絕對不能用整個數(shù)據(jù)集計(jì)算均值和標(biāo)準(zhǔn)差然后統(tǒng)一標(biāo)準(zhǔn)化否則測試集信息就“泄露”到訓(xùn)練過程了會嚴(yán)重高估模型性能。類別不平衡如果四個情感類別樣本數(shù)差異很大需要在fitcsvm中設(shè)置Prior參數(shù)為empirical讓模型根據(jù)訓(xùn)練數(shù)據(jù)自動調(diào)整先驗(yàn)概率或者設(shè)置Weights參數(shù)給少數(shù)類更高的權(quán)重。網(wǎng)格搜索的代價(jià)參數(shù)網(wǎng)格C_values和gamma_values的乘積乘以交叉驗(yàn)證折數(shù)就是需要訓(xùn)練模型的次數(shù)。如果數(shù)據(jù)量大、特征多這個過程會非常耗時?？梢詮拇志W(wǎng)格開始如[0.01, 0.1, 1, 10, 100]找到大致范圍后再在附近進(jìn)行精細(xì)搜索。也可以使用更高效的優(yōu)化算法如貝葉斯優(yōu)化bayesopt函數(shù)。MATLAB的KernelScale注意fitcsvm的KernelScale參數(shù)對應(yīng)的是RBF核公式中的sigma而通常我們說的gamma參數(shù)是1/(2*sigma^2)。所以代碼中我們用1/sqrt(gamma)來轉(zhuǎn)換。務(wù)必理清你所用工具包對核參數(shù)的定義。5. 模型評估、結(jié)果分析與可視化得到一個準(zhǔn)確率數(shù)字只是開始更重要的是理解模型在哪里做得好在哪里做得不好。5.1 超越準(zhǔn)確率全面的評估指標(biāo)對于四分類問題混淆矩陣Confusion Matrix是最直觀的工具。% 計(jì)算混淆矩陣 C confusionmat(y_test, y_pred_final); % 使用內(nèi)置函數(shù)繪圖更美觀 figure; confusionchart(C, {LALV, LAHV, HALV, HAHV}); % 按你的標(biāo)簽順序 title(SVM分類器混淆矩陣);從混淆矩陣中我們可以計(jì)算每個類別的精確率Precision、召回率Recall和F1分?jǐn)?shù)。% 計(jì)算每類的性能指標(biāo) num_classes 4; precision zeros(num_classes, 1); recall zeros(num_classes, 1); f1_score zeros(num_classes, 1); for i 1:num_classes TP C(i, i); FP sum(C(:, i)) - TP; FN sum(C(i, :)) - TP; precision(i) TP / (TP FP eps); recall(i) TP / (TP FN eps); f1_score(i) 2 * (precision(i) * recall(i)) / (precision(i) recall(i) eps); end % 顯示結(jié)果 class_names {LALV, LAHV, HALV, HAHV}; for i 1:num_classes fprintf(類別 %s: 精確率%.3f, 召回率%.3f, F1%.3f\n, ... class_names{i}, precision(i), recall(i), f1_score(i)); end % 計(jì)算宏平均F1 macro_f1 mean(f1_score); fprintf(宏平均F1分?jǐn)?shù): %.3f\n, macro_f1);為什么F1分?jǐn)?shù)很重要在類別不平衡的情況下準(zhǔn)確率可能會被大類別主導(dǎo)而虛高。比如90%的樣本都是“平靜”模型只要全部預(yù)測“平靜”就能得到90%的準(zhǔn)確率但這毫無意義。F1分?jǐn)?shù)是精確率和召回率的調(diào)和平均能更好地衡量模型對每個類別的識別能力宏平均F1對所有類別一視同仁。5.2 結(jié)果可視化與解讀特征重要性分析對于線性SVM可以通過檢查權(quán)重向量w的絕對值大小來估計(jì)特征的重要性。權(quán)重絕對值越大說明該特征對決策的影響越大。if strcmp(final_svm_model.KernelParameters.Function, linear) w final_svm_model.Beta; % 權(quán)重系數(shù) [sorted_w, idx] sort(abs(w), descend); top_n 20; fprintf(重要性排名前%d的特征索引:\n, top_n); disp(idx(1:top_n)); % 你可以將這些索引映射回具體的通道和特征類型如“F3通道的Alpha功率” end決策邊界可視化降維后對于高維數(shù)據(jù)我們可以使用t-SNE或PCA將其降到2維或3維然后畫出樣本點(diǎn)和SVM的決策邊界需要訓(xùn)練一個在降維空間上的新SVM這有助于直觀感受數(shù)據(jù)的可分性。% 使用PCA降維到2維 [coeff, score, latent] pca(X_train_scaled); X_train_2d score(:, 1:2); % 在2維數(shù)據(jù)上重新訓(xùn)練一個SVM僅用于可視化 svm_2d fitcsvm(X_train_2d, y_train, KernelFunction, rbf, Standardize, false); % 創(chuàng)建網(wǎng)格用于繪制決策區(qū)域 x1range min(X_train_2d(:,1))-1 : 0.1 : max(X_train_2d(:,1))1; x2range min(X_train_2d(:,2))-1 : 0.1 : max(X_train_2d(:,2))1; [xx1, xx2] meshgrid(x1range, x2range); XGrid [xx1(:), xx2(:)]; % 預(yù)測網(wǎng)格點(diǎn)的類別 [~, scores] predict(svm_2d, XGrid); % 繪制決策區(qū)域和樣本點(diǎn) figure; h gscatter(X_train_2d(:,1), X_train_2d(:,2), y_train, brgk, o*^v); hold on; contour(xx1, xx2, reshape(scores(:,4), size(xx1)), [0 0], k-, LineWidth, 2); % 繪制第4類HAHV的決策邊界 xlabel(PCA主成分1); ylabel(PCA主成分2); title(PCA降維后的數(shù)據(jù)分布與SVM決策邊界示例); legend(LALV,LAHV,HALV,HAHV, Decision Boundary); hold off;5.3 與基線模型和現(xiàn)有研究的對比為了評估你的SVM模型是否有效需要設(shè)立基線進(jìn)行比較。隨機(jī)猜測基線四分類問題的隨機(jī)猜測準(zhǔn)確率是25%。多數(shù)類基線預(yù)測為訓(xùn)練集中樣本數(shù)最多的類別計(jì)算其在測試集上的準(zhǔn)確率。簡單分類器基線如最近鄰分類器KNNK1或決策樹最大深度5。% 多數(shù)類基線 majority_class mode(y_train); y_pred_majority repmat(majority_class, size(y_test)); accuracy_majority sum(y_pred_majority y_test) / numel(y_test); fprintf(多數(shù)類基線準(zhǔn)確率: %.2f%%\n, accuracy_majority * 100); % 1-NN基線 knn_model fitcknn(X_train_scaled, y_train, NumNeighbors, 1, Standardize, false); y_pred_knn predict(knn_model, X_test_scaled); accuracy_knn sum(y_pred_knn y_test) / numel(y_test); fprintf(1-NN基線準(zhǔn)確率: %.2f%%\n, accuracy_knn * 100);如果你的SVM模型顯著優(yōu)于這些基線例如準(zhǔn)確率超過40%F1分?jǐn)?shù)超過0.4并且與相關(guān)論文中報(bào)告的基于DEAP的SVM性能通常在55%-70%之間取決于特征、通道選擇和分類任務(wù)處于同一量級那么你的工作就是有效的。6. 項(xiàng)目總結(jié)、局限性與進(jìn)階方向走完整個流程從數(shù)據(jù)預(yù)處理、特征提取到模型訓(xùn)練評估一個完整的基于SVM的腦電情感識別項(xiàng)目就完成了?；仡櫿麄€過程有幾個關(guān)鍵點(diǎn)值得再次強(qiáng)調(diào)預(yù)處理是基石濾波和去噪的質(zhì)量直接決定了特征提取的上限。務(wù)必根據(jù)腦電信號的特性選擇合適的濾波參數(shù)。特征工程是核心時頻域特征特別是頻帶功率、微分熵對情感識別非常有效。嘗試組合不同類型的特征時域、頻域、非線性特征有時能提升性能但要注意維度控制。模型調(diào)優(yōu)是藝術(shù)SVM的性能對超參數(shù)和特征尺度敏感。網(wǎng)格搜索交叉驗(yàn)證是標(biāo)準(zhǔn)做法但計(jì)算成本高。務(wù)必做好數(shù)據(jù)標(biāo)準(zhǔn)化防止泄露。評估要全面不要只看準(zhǔn)確率混淆矩陣和F1分?jǐn)?shù)能揭示模型在各類別上的具體表現(xiàn)尤其是在數(shù)據(jù)不平衡時。本項(xiàng)目的局限性特征手工設(shè)計(jì)我們提取的時頻特征是手工設(shè)計(jì)的依賴于領(lǐng)域知識。這可能不是最優(yōu)的表示。通道選擇我們使用了所有通道但有些通道可能對情感識別貢獻(xiàn)很小甚至引入噪聲。通道選擇是一個重要的優(yōu)化方向。個體差異腦電信號存在巨大的個體差異。在一個被試上訓(xùn)練好的模型在另一個被試上可能效果很差。這就是“被試獨(dú)立”評估的挑戰(zhàn)。SVM的局限性SVM雖然強(qiáng)大但對于非常復(fù)雜的非線性關(guān)系其表現(xiàn)可能不如深度神經(jīng)網(wǎng)絡(luò)。進(jìn)階方向深度學(xué)習(xí)端到端學(xué)習(xí)使用卷積神經(jīng)網(wǎng)絡(luò)直接從原始腦電信號或時頻圖中學(xué)習(xí)特征避免手工特征工程。例如EEGNet、DeepConvNet等網(wǎng)絡(luò)結(jié)構(gòu)在此任務(wù)上表現(xiàn)出色?？绫辉噷W(xué)習(xí)與域自適應(yīng)研究如何讓模型更好地泛化到新被試?？梢允褂糜蜃赃m應(yīng)、遷移學(xué)習(xí)或元學(xué)習(xí)等方法。多模態(tài)融合DEAP數(shù)據(jù)集除了EEG還有外周生理信號GSR, RESP, EMG等。融合多模態(tài)信息可以顯著提升情感識別的魯棒性和準(zhǔn)確性。時序建模情感是動態(tài)變化的使用循環(huán)神經(jīng)網(wǎng)絡(luò)或Transformer對腦電信號的時間序列進(jìn)行建?？梢圆蹲狡鋭討B(tài)演化模式。更精細(xì)的情感模型除了四分類也可以嘗試回歸任務(wù)直接預(yù)測效價(jià)和喚醒度的連續(xù)值或者探索更復(fù)雜的情感模型。這個項(xiàng)目提供了一個堅(jiān)實(shí)的起點(diǎn)。你可以基于這個代碼框架嘗試不同的特征、不同的分類器、不同的預(yù)處理方法逐步深入腦電情感識別這個有趣且充滿挑戰(zhàn)的領(lǐng)域。代碼的模塊化設(shè)計(jì)使得替換特征提取模塊或分類器模塊變得非常容易祝你實(shí)驗(yàn)順利。本文還有配套的精品資源點(diǎn)擊獲取