戰(zhàn)：從信號(hào)預(yù)處理到特征提取的完整流程)
看到“心音分類器”這個(gè)題目可能有朋友第一反應(yīng)是這不是典型的深度學(xué)習(xí)活兒?jiǎn)崮脗€(gè)預(yù)訓(xùn)練網(wǎng)絡(luò)懟上去精度刷到99%完事。但真當(dāng)你拿到一份真實(shí)的生理信號(hào)數(shù)據(jù)時(shí)會(huì)發(fā)現(xiàn)事情遠(yuǎn)沒有這么簡(jiǎn)單。信號(hào)里有環(huán)境噪聲、有呼吸音、有傳感器摩擦聲甚至還有醫(yī)生操作時(shí)不小心碰觸麥克風(fēng)的聲音直接丟給任何模型都是災(zāi)難。我去年在社團(tuán)周會(huì)上帶新人做機(jī)器學(xué)習(xí)入門項(xiàng)目選了一圈之后最后定在了“MATLAB 心音分類”這條線上折騰了大概三周最大的感受是用MATLAB做這個(gè)項(xiàng)目真正的價(jià)值點(diǎn)根本不在“分類器”那一環(huán)而是前面的信號(hào)預(yù)處理和特征提取——那才是決定模型上限的地方。這篇就把我們的完整流程、核心代碼和踩過的坑一次說清楚適合剛接觸機(jī)器學(xué)習(xí)、想在真實(shí)數(shù)據(jù)上完整跑通“信號(hào) → 特征 → 模型 → 評(píng)估”這條鏈路的朋友參考。1. 為什么是MATLAB為什么是心音信號(hào)1.1 這個(gè)項(xiàng)目在做什么把生物信號(hào)變成機(jī)器學(xué)習(xí)問題心音分類在最簡(jiǎn)單版本里是一個(gè)二分類問題給定一段心音錄音判斷它是“正?！边€是“異常”。異常心音包括各種雜音、瓣膜問題導(dǎo)致的額外音、心律不齊等等。醫(yī)學(xué)上醫(yī)生用聽診器聽診靠經(jīng)驗(yàn)判斷而我們要做的就是用機(jī)器學(xué)習(xí)模型去模擬這個(gè)判斷過程。選這個(gè)題作為機(jī)器學(xué)習(xí)初探項(xiàng)目有幾個(gè)天然優(yōu)勢(shì)。第一信號(hào)非常直觀——你把它畫出來正常人能看出明顯的“嘣噠”兩下心跳節(jié)律異常信號(hào)往往波形亂、幅值不均勻這個(gè)可視化過程對(duì)建立直覺特別有幫助。第二公開數(shù)據(jù)集很成熟PhysioNet 2016心音分類挑戰(zhàn)賽的數(shù)據(jù)集有三千多條錄音wav格式采樣率2000Hz文件組織和標(biāo)簽都齊。第三問題邊界清晰特征提取有明確的信號(hào)處理理論支撐分類器選擇可以從最基礎(chǔ)的決策樹一直玩到集成學(xué)習(xí)梯度非常平滑不會(huì)出現(xiàn)“要么完全不會(huì)、要么直接上深度學(xué)習(xí)黑盒”的斷層感。這個(gè)項(xiàng)目我們最終的目標(biāo)不是刷高精度數(shù)值而是讓社團(tuán)里的同學(xué)完整經(jīng)歷一遍機(jī)器學(xué)習(xí)的標(biāo)準(zhǔn)流程拿到原始數(shù)據(jù)、做清洗和預(yù)處理、人工設(shè)計(jì)特征、訓(xùn)練模型、用交叉驗(yàn)證評(píng)估泛化能力。走完這一步后面換任何領(lǐng)域的數(shù)據(jù)——比如電機(jī)振動(dòng)信號(hào)、語音命令、心電圖——套路都是通的。1.2 MATLAB在這條鏈路里的獨(dú)特優(yōu)勢(shì)我知道現(xiàn)在一說機(jī)器學(xué)習(xí)大家第一反應(yīng)是Python。但在這個(gè)項(xiàng)目里MATLAB其實(shí)是更有性價(jià)比的選擇尤其是在做信號(hào)處理和特征提取的階段。Python做信號(hào)處理不是不行但需要numpy、scipy、librosa、pandas、scikit-learn等等好幾個(gè)包配合著裝版本之間存在兼容性問題光是給新人配環(huán)境就能耗掉半天。而MATLAB一個(gè)軟件就把這些全干了信號(hào)處理工具箱里有濾波、傅里葉變換、小波分析、倒譜分析的完整函數(shù)統(tǒng)計(jì)與機(jī)器學(xué)習(xí)工具箱里有fitctree決策樹、fitcecoc多分類SVM、fitcknnKNN、fitcensemble集成學(xué)習(xí)等現(xiàn)成模型畫圖方面MATLAB的原生繪圖函數(shù)比matplotlib默認(rèn)樣式好看一個(gè)量級(jí)而且操作不需要寫一堆樣式調(diào)整代碼。另外還有一點(diǎn)在初學(xué)階段特別關(guān)鍵MATLAB的變量工作區(qū)是可視化的。你可以隨時(shí)雙擊變量看數(shù)據(jù)長(zhǎng)什么樣、尺寸是多少、是單精度還是雙精度這在調(diào)試特征矩陣維度的時(shí)侯簡(jiǎn)直救命。Python IDE里你還要print shape或者用debuggerMATLAB直接一目了然。MATLAB當(dāng)然也有它的缺點(diǎn)比如閉源、貴、深度學(xué)習(xí)的生態(tài)確實(shí)不如Python豐富。但對(duì)于“機(jī)器學(xué)習(xí)初探”這個(gè)定位來說這些缺點(diǎn)完全不影響。學(xué)生用校園授權(quán)版或者開源社區(qū)里跟Octave的兼容性方案都?jí)蛴?。我個(gè)人的態(tài)度是工具是手段不是信仰哪個(gè)能讓你最快把核心邏輯跑通就用哪個(gè)。這個(gè)項(xiàng)目的核心邏輯是信號(hào)處理和特征工程用MATLAB能讓你的精力百分之百集中在核心邏輯而不是環(huán)境配置上這就是我推薦它的理由。2. 心音數(shù)據(jù)的獲取與預(yù)處理先把“雜音”里的信息洗出來2.1 數(shù)據(jù)要從哪里來PhysioNet 2016數(shù)據(jù)集我們用的數(shù)據(jù)集是PhysioNet/Computing in Cardiology Challenge 2016提供的心音分類公開數(shù)據(jù)集這是一套專門用于心音自動(dòng)分析研究的基準(zhǔn)數(shù)據(jù)。數(shù)據(jù)分為training-a到training-e五個(gè)子文件夾一共包含三千多條心音錄音wav文件每個(gè)文件夾下面還有一個(gè)HEART_SOUND.csv文件里面是每條錄音對(duì)應(yīng)的標(biāo)簽1表示正常-1表示異常。先說幾個(gè)基本參數(shù)錄音統(tǒng)一是2000Hz采樣率也就是說每秒鐘有2000個(gè)采樣點(diǎn)人耳能聽到的頻率范圍是20Hz到20000Hz而心音的絕大部分能量集中在20Hz到400Hz這個(gè)低頻段2000Hz采樣率對(duì)心音來說信息量完全足夠。讀取文件用MATLAB的audioread函數(shù)即可。% 讀取單條心音錄音 [y, fs] audioread(training-a/a0001.wav); % y是采樣點(diǎn)數(shù)組fs是采樣率(2000Hz) disp(size(y)); disp(fs);這里有一個(gè)新手常犯的錯(cuò)誤錄音文件不是等長(zhǎng)的。有些錄音3秒有些8秒文件總時(shí)長(zhǎng)差異很大。如果直接把整段錄音作為輸入不同樣本的特征維度就對(duì)不上。所以預(yù)處理階段要統(tǒng)一處理到固定長(zhǎng)度或者做分段切割后面我會(huì)細(xì)說。2.2 預(yù)處理三步走降采樣、濾波、歸一化拿到原始錄音之后不要急著提特征。原始信號(hào)里包含了大量與心音無關(guān)的“臟東西”我把它稱為噪聲三件套環(huán)境噪聲空調(diào)聲、房間混響、生理干擾呼吸聲、肌肉信號(hào)、傳感器噪聲麥克風(fēng)接觸摩擦、電噪聲。此外還有50Hz/60Hz的工頻干擾。我用的預(yù)處理流程是三步降采樣可選、帶通濾波、幅值歸一化。降采樣既然心音主要能量在400Hz以下根據(jù)奈奎斯特采樣定理理論上800Hz采樣率就足夠重建信號(hào)。把2000Hz降采樣到1000Hz甚至800Hz好處是數(shù)據(jù)量減半后續(xù)計(jì)算更快。不過對(duì)入門項(xiàng)目來說2000Hz的數(shù)據(jù)量本身也不大降不降采樣對(duì)最終結(jié)果影響很小所以這一步我們后來基本跳過。帶通濾波這是我們預(yù)處理里的核心步驟。心音的診斷信息主要位于25Hz到400Hz之間因此設(shè)計(jì)了一個(gè)四階巴特沃斯帶通濾波器保留這個(gè)頻段濾掉低頻漂移和高頻噪聲。低頻漂移特別值得注意因?yàn)閭鞲衅鹘佑|或呼吸動(dòng)作會(huì)產(chǎn)生一種非常慢的基線偏移如果不濾掉后面提取的時(shí)域特征比如均值和方差會(huì)被這個(gè)漂移帶偏。% 設(shè)計(jì)四階巴特沃斯帶通濾波器通帶25Hz~400Hz fs 2000; fL 25; fH 400; [b, a] butter(4, [fL fH]/(fs/2), bandpass); % 使用filtfilt零相位濾波避免相位失真 y_filtered filtfilt(b, a, y);為什么用filtfilt而不用filter這是很多初學(xué)會(huì)忽略的細(xì)節(jié)。filter是常規(guī)濾波會(huì)帶來相位偏移filtfilt先把信號(hào)正向?yàn)V一遍、再反向?yàn)V一遍兩次的相位偏移互相抵消得到零相位失真的結(jié)果。對(duì)心音這種看重波形形態(tài)的分析來說絕對(duì)有必要。幅值歸一化因?yàn)椴煌浺舻匿浿圃O(shè)備、深淺閾值、病人體型都不同音量幅值差異非常大。歸一化就是把每個(gè)樣本的最大幅值縮放為1。這樣后續(xù)特征比如峰值、能量才有跨樣本可比性。% 將幅值歸一化到[-1, 1] y_normalized y_filtered / max(abs(y_filtered));我見過一個(gè)很隱蔽的坑如果你先歸一化、再濾波濾波后的信號(hào)幅值可能又超過1了因?yàn)闉V波本身會(huì)改變信號(hào)的能量分布。所以順序應(yīng)該是先濾波、再歸一化。這個(gè)小細(xì)節(jié)不影響大的方向但會(huì)讓你的特征數(shù)值更干凈。2.3 信號(hào)分段與樣本構(gòu)建從連續(xù)波形到一個(gè)個(gè)“樣本”現(xiàn)在手里是清洗過的完整錄音。但問題來了一條錄音少說也有好幾秒心音是典型的周期信號(hào)——正常心跳每分鐘60到100次大約每0.6秒到1秒一個(gè)周期。整段錄音當(dāng)成一個(gè)樣本信息太稠、特征糊成一片分類器根本找不到判別性信息。更合理的做法是分段。我們把每條錄音切成了固定長(zhǎng)度為3秒的片段當(dāng)作獨(dú)立的樣本參與后續(xù)訓(xùn)練。為什么是3秒經(jīng)驗(yàn)值是3秒至少包含3個(gè)心動(dòng)周期既能捕捉到周期性和一致性又不至于因?yàn)榘芷谔喽旬惓＜?xì)節(jié)平均掉。% 將信號(hào)切分為3秒一段 segLen 3 * fs; % 3秒 numSeg floor(length(y_normalized) / segLen); for i 1:numSeg seg y_normalized((i-1)*segLen 1 : i*segLen); % 將seg作為一條樣本存入后續(xù)的特征提取環(huán)節(jié) end這里有個(gè)取舍需要說清楚。理想情況下我們應(yīng)該是先做心音分割識(shí)別出S1、收縮期、S2、舒張期這四段再?gòu)拿總€(gè)心動(dòng)周期提特征這樣特征和生理機(jī)制對(duì)應(yīng)得更準(zhǔn)。但心音分割本身就是一個(gè)不小的課題需要做包絡(luò)檢測(cè)、峰值定位、邏輯判斷相當(dāng)于在分類問題之前又疊了一個(gè)完整項(xiàng)目對(duì)初學(xué)者來說步子邁得太大。我們最終采取的方案是“簡(jiǎn)單粗暴按時(shí)間窗切”訓(xùn)練出來的模型效果其實(shí)已經(jīng)夠用了。初探階段先跑通流程比什么都重要。3. 特征提取把心跳聲翻譯成機(jī)器認(rèn)識(shí)的數(shù)字3.1 時(shí)域特征如何描述一段波形的“性格”信號(hào)本身是幾千個(gè)采樣點(diǎn)機(jī)器學(xué)習(xí)模型不能直接吃這種原始波形。原始波形就算長(zhǎng)度統(tǒng)一維度也太大3秒×2000采樣率6000維而且每個(gè)采樣點(diǎn)的具體數(shù)值受噪聲影響很大不夠魯棒。特征提取的目標(biāo)就是把一段波形壓縮成一組維度不高、但能描述其本質(zhì)屬性的數(shù)字。我管時(shí)域特征叫“波形性格參數(shù)”它們是直接在時(shí)間軸上計(jì)算出的統(tǒng)計(jì)量均方根值RMS描述信號(hào)整體能量大小。心音異常時(shí)往往伴隨額外雜音信號(hào)能量分布會(huì)變化。方差描述信號(hào)波動(dòng)程度。異常心音經(jīng)常更“亂”方差會(huì)變大。峰值/谷值描述單次心跳的最大沖擊力。S1第一心音和S2第二心音的峰值比例在不同疾病下會(huì)有變化。過零率信號(hào)在單位時(shí)間內(nèi)穿越零點(diǎn)的次數(shù)。過零率越高說明信號(hào)越“高頻”異常雜音通常是高頻的會(huì)讓過零率顯著上升。峰值因子峰值除以RMS描述波形“尖銳”程度。心音異常時(shí)心臟射血狀態(tài)改變波形形態(tài)隨之改變。function feat extractTimeFeatures(seg) feat zeros(1, 5); feat(1) rms(seg); % 均方根值 feat(2) var(seg); % 方差 feat(3) max(abs(seg)); % 峰值 feat(4) sum(abs(diff(sign(seg)))) / (length(seg)-1); % 過零率 feat(5) max(abs(seg)) / (rms(seg) eps); % 峰值因子 end這幾個(gè)人人都能算、看起來“沒什么技術(shù)含量”的特征在實(shí)際分類中的貢獻(xiàn)往往出乎意料地大。尤其是過零率和方差訓(xùn)練出來的隨機(jī)森林模型里重要性排得都很靠前。原因也簡(jiǎn)單異常心音中高頻雜音成分多信號(hào)振蕩更頻繁過零率自然高。你不需要理解它的生理機(jī)理只需要信任它能區(qū)分正常和異常。3.2 頻域特征與倒譜MFCC為什么能跨界單獨(dú)靠時(shí)域特征還不夠——不同病人心音的時(shí)域波形即使同是“正?！毙螒B(tài)差異也很大直接用原始幅值統(tǒng)計(jì)量做特征很難做到跨個(gè)體泛化。頻域特征的加入相當(dāng)于從另一個(gè)角度描述信號(hào)不再是“這個(gè)波形長(zhǎng)什么樣”而是“這個(gè)波形由哪些頻率成分構(gòu)成”。首先做短時(shí)傅里葉變換得到語譜圖?；谡Z譜圖可以計(jì)算頻域特征頻譜質(zhì)心信號(hào)能量在頻率上的加權(quán)平均描述“聲音聽起來亮不亮”。頻譜帶寬能量在質(zhì)心周圍的散布程度異常雜音會(huì)讓帶寬變大。頻譜滾降點(diǎn)信號(hào)能量有多少集中在低頻端的閾值點(diǎn)心音是低頻為主滾降點(diǎn)通常在低頻。然后是MFCC梅爾頻率倒譜系數(shù)。這個(gè)詞大家可能在語音識(shí)別里聽過但心音也是聲學(xué)信號(hào)MFCC同樣適用。MFCC把頻譜按人耳感知的梅爾刻度映射到對(duì)數(shù)刻度然后做離散余弦變換得到的倒譜系數(shù)能非常緊湊地表示頻譜包絡(luò)的形狀。心音的頻譜包絡(luò)跟心臟結(jié)構(gòu)、瓣膜狀態(tài)有關(guān)異常狀態(tài)下包絡(luò)形狀變化明顯MFCC能更好捕捉這種變化。% 計(jì)算13維MFCC特征 [coeff, ~] mfcc(seg, fs, NumCoeffs, 13, ... WindowLength, round(0.03*fs), OverlapLength, round(0.02*fs)); % 沿時(shí)間軸取平均值和標(biāo)準(zhǔn)差得到26維特征 mfccMean mean(coeff, 1); mfccStd std(coeff, 0, 1); featMFCC [mfccMean, mfccStd];MFCC的重要性在于它經(jīng)過語音識(shí)別領(lǐng)域幾十年驗(yàn)證是一種對(duì)“聲音”來說極其高效的壓縮表示。心音本質(zhì)上也是聲音共享這套表示完全合理。用在我們的模型里加入MFCC后交叉驗(yàn)證準(zhǔn)確率大約提升了4到5個(gè)百分點(diǎn)是貢獻(xiàn)最大的一類特征。3.3 用MATLAB組裝特征矩陣別在維度上翻車把所有片段的特征拼成一個(gè)矩陣這是機(jī)器學(xué)習(xí)訓(xùn)練的前置工作也是很多新手第一次接觸“特征矩陣”這個(gè)概念的地方。這里有一個(gè)非常重要的規(guī)范每個(gè)樣本提取的特征必須是相同長(zhǎng)度的向量。比如時(shí)域特征選了5個(gè)頻域特征選了5個(gè)MFCC選了26個(gè)那這個(gè)樣本的特征向量就是1×36。所有樣本都提這36個(gè)特征最后拼出一個(gè)N×36的特征矩陣。featureMatrix []; labels []; for i 1:length(fileList) [y, fs] audioread(fileList{i}); % 預(yù)處理濾波、歸一化 y filtfilt(b, a, y); y y / max(abs(y)); % 分段 segLen 3 * fs; numSeg floor(length(y) / segLen); for j 1:numSeg seg y((j-1)*segLen1 : j*segLen); featT extractTimeFeatures(seg); % 5維 featF extractFreqFeatures(seg, fs); % 5維 featM extractMFCCFeatures(seg, fs); % 26維 featAll [featT, featF, featM]; % 36維 featureMatrix [featureMatrix; featAll]; labels [labels; labelVal]; % 0或1 end end % 保存到mat文件供后續(xù)訓(xùn)練復(fù)用 save(heart_feature.mat, featureMatrix, labels, -v7.3);這段代碼里我特意標(biāo)了一個(gè)重要的壞習(xí)慣示例特征矩陣和標(biāo)簽用[matrix; newRow]這種橫向拼接方式。在小數(shù)據(jù)量下它能工作但樣本上萬的時(shí)候這種寫法會(huì)反復(fù)復(fù)制內(nèi)存慢得讓人抓狂。正確做法是預(yù)分配矩陣或者把每個(gè)樣本的結(jié)果存到cell里最后一次性合并。被教育過一次之后我現(xiàn)在寫的代碼都是先zeros(N, 36)初始化好再填充。這算是個(gè)性能層面的教訓(xùn)。組裝好之后建議立刻做一步特征標(biāo)準(zhǔn)化。什么是標(biāo)準(zhǔn)化就是把每個(gè)特征的均值變?yōu)?、標(biāo)準(zhǔn)差變?yōu)?。featureMatrixNorm normalize(featureMatrix);為什么標(biāo)準(zhǔn)化是必須的因?yàn)椴煌卣髁烤V差異巨大——MFCC的數(shù)值可能是幾十而RMS只有零點(diǎn)幾。很多分類器內(nèi)部計(jì)算距離或正則化項(xiàng)會(huì)讓量綱大的特征主導(dǎo)結(jié)果。標(biāo)準(zhǔn)化就是消除量綱影響。需要提醒的是標(biāo)準(zhǔn)化參數(shù)只能在訓(xùn)練集上計(jì)算測(cè)試集要用同樣的均值和標(biāo)準(zhǔn)差去轉(zhuǎn)換這個(gè)細(xì)節(jié)放到第5節(jié)交叉驗(yàn)證的時(shí)候再具體展開。4. 分類器選型與訓(xùn)練從基線模型開始的進(jìn)階路線4.1 為什么不要一上來就上深度學(xué)習(xí)把特征矩陣準(zhǔn)備好之后接下來就是模型訓(xùn)練。但這里我要反復(fù)強(qiáng)調(diào)一個(gè)初學(xué)階段的觀念不要一上來就上深度學(xué)習(xí)。神經(jīng)網(wǎng)絡(luò)是一個(gè)高度非線性的“萬能函數(shù)”但它需要海量數(shù)據(jù)、大量調(diào)參、以及在驗(yàn)證集上的迭代監(jiān)控。我們的心音數(shù)據(jù)集經(jīng)過分段后樣本量大約幾千條特征維度只有36維這種“小樣本中等維度”的設(shè)定本身就是傳統(tǒng)機(jī)器學(xué)習(xí)的主場(chǎng)。用傳統(tǒng)機(jī)器學(xué)習(xí)有幾個(gè)看得見摸得著的好處第一訓(xùn)練快幾十秒出結(jié)果可以反復(fù)試第二可解釋性強(qiáng)比如決策樹你可以把樹畫出來看它到底學(xué)了什么規(guī)則第三對(duì)特征工程的好壞能立刻得到反饋。你加了一組特征、準(zhǔn)確率提升了你知道是這組特征起了作用你換了個(gè)亂七八糟的特征、準(zhǔn)確率下降了你也知道是這組特征在拖后腿。這種反饋對(duì)建立直覺極其重要是深度學(xué)習(xí)的黑盒給不了的。我在表里列了四種適合初學(xué)的分類器以及它們?cè)诒卷?xiàng)目上的表現(xiàn)對(duì)比分類器原理一句話本項(xiàng)目初測(cè)準(zhǔn)確率訓(xùn)練速度調(diào)參難度決策樹用一系列“特征閾值”的判斷做分類約78%極快低K近鄰KNN看距離最近的K個(gè)樣本屬于哪類約82%訓(xùn)練快、預(yù)測(cè)慢低支持向量機(jī)SVM在高維空間找一個(gè)最優(yōu)分割超平面約84%快中隨機(jī)森林訓(xùn)練多棵決策樹投票決定結(jié)果約88%快中低AdaBoost集成用多個(gè)弱分類器加權(quán)組合成強(qiáng)分類器約90%中等中從這個(gè)表能明顯看出規(guī)律集成方法隨機(jī)森林、AdaBoost比單模型表現(xiàn)好。原因也很直觀心音特征內(nèi)部往往存在復(fù)雜的交互關(guān)系比如某個(gè)時(shí)域特征在結(jié)合某個(gè)頻域特征的時(shí)候才有判別意義單棵決策樹對(duì)這個(gè)交互的建模能力是有限的但多棵樹從不同角度切分?jǐn)?shù)據(jù)再投票能把這種交互關(guān)系綜合起來。4.2 三種主流分類器的MATLAB實(shí)現(xiàn)與調(diào)用MATLAB的fitc系列函數(shù)把傳統(tǒng)機(jī)器學(xué)習(xí)封裝得極其友好。下面給出我們項(xiàng)目里用到的決策樹、SVM和集成學(xué)習(xí)三種實(shí)現(xiàn)的對(duì)照代碼讀者拿去改數(shù)據(jù)路徑就能跑通。% 假設(shè)featureNorm是N×36特征矩陣label是N×1標(biāo)簽0或1 % 1. 決策樹 mdlTree fitctree(featureNorm, label, CrossVal, off); accuracyTree 1 - resubLoss(mdlTree); % 注意這只是訓(xùn)練集回代精度 % 2. SVMRBF核 mdlSVM fitcecoc(featureNorm, label, Learners, svm); % fitcecoc是錯(cuò)誤糾錯(cuò)輸出碼用于多分類二分類時(shí)等價(jià)于一個(gè)SVM % 3. 集成學(xué)習(xí)AdaBoostM1 決策樹弱學(xué)習(xí)器 t templateTree(MaxNumSplits, 10); mdlEnsemble fitcensemble(featureNorm, label, ... Method, AdaBoostM2, ... NumLearningCycles, 50, ... Learners, t);fitcensemble是這里面的明星函數(shù)。它內(nèi)部自動(dòng)處理了自適應(yīng)增強(qiáng)的各個(gè)細(xì)節(jié)你只需要指定弱學(xué)習(xí)器和方法。我們指定MaxNumSplits10是為了限制每棵弱樹不要太深太深的樹容易過擬合保留一些“弱能力”的模型反而讓Adaboost更有效。訓(xùn)練完之后預(yù)測(cè)新數(shù)據(jù)的邏輯更簡(jiǎn)單% 對(duì)測(cè)試集預(yù)測(cè) predLabel predict(mdlEnsemble, testFeatureNorm); % 計(jì)算準(zhǔn)確率 acc sum(predLabel testLabel) / length(testLabel);4.3 集成模型的參數(shù)調(diào)節(jié)經(jīng)驗(yàn)機(jī)器學(xué)習(xí)里有個(gè)概念叫“偏差-方差權(quán)衡”模型太簡(jiǎn)單會(huì)欠擬合偏差高模型太復(fù)雜會(huì)過擬合方差高。集成模型的設(shè)計(jì)初衷就是通過組合多個(gè)弱模型在保持低偏差的同時(shí)降低方差。用fitcensemble的時(shí)候有兩個(gè)參數(shù)最值得調(diào)NumLearningCycles弱學(xué)習(xí)器數(shù)量和Learners里的深度限制。我自己的經(jīng)驗(yàn)是弱學(xué)習(xí)器數(shù)量從10加到50時(shí)效果提升最明顯到100以上提升就趨于平緩了。收益邊際遞減這個(gè)現(xiàn)象在集成學(xué)習(xí)里非常典型。弱樹的深度控制在5到10層比較合適。太深容易單棵樹就過擬合了太淺每棵樹都太笨集成出來的模型不夠聰明。在深度學(xué)習(xí)特征的標(biāo)準(zhǔn)化問題上對(duì)樹模型來說標(biāo)準(zhǔn)化影響不大因?yàn)闃渲魂P(guān)心閾值切分不關(guān)心距離但對(duì)SVM和KNN這種基于距離的模型影響巨大。我們的做法是找出全部特征直接標(biāo)準(zhǔn)化一視同仁。調(diào)參的時(shí)候不要靠猜可以用MATLAB的fitcensemble配合OptimizeHyperparameters做自動(dòng)超參數(shù)搜索。不過要注意超參數(shù)搜索會(huì)跑很久數(shù)據(jù)量大的時(shí)候建議先抽取一部分?jǐn)?shù)據(jù)做快速實(shí)驗(yàn)確定參數(shù)范圍后再全量訓(xùn)練。% 自動(dòng)超參數(shù)調(diào)優(yōu)示例小心運(yùn)行時(shí)間 mdlOpt fitcensemble(featureNorm, label, ... Method, AdaBoostM2, ... OptimizeHyperparameters, {NumLearningCycles, LearnRate}, ... HyperparameterOptimizationOptions, struct(MaxObjectiveEvaluations, 30));5. 評(píng)估與調(diào)參準(zhǔn)確率之外還要看敏感度和特異性5.1 混淆矩陣與三類指標(biāo)別被準(zhǔn)確率蒙蔽訓(xùn)練完模型很多人只會(huì)看一個(gè)準(zhǔn)確率數(shù)字。但在這個(gè)項(xiàng)目里我要強(qiáng)調(diào)準(zhǔn)確率是一個(gè)不夠用的指標(biāo)。原因在于正負(fù)樣本不均衡心音數(shù)據(jù)集中正常樣本比異常樣本略多即使你無腦把全部樣本輸出成“正常”準(zhǔn)確率也能到差不多52%。在這種背景下準(zhǔn)確率不能描述模型在“少數(shù)類”上的表現(xiàn)。我們需要看混淆矩陣它把預(yù)測(cè)結(jié)果和真實(shí)標(biāo)簽的對(duì)應(yīng)關(guān)系展開成一個(gè)2×2表真實(shí)\預(yù)測(cè)預(yù)測(cè)正常預(yù)測(cè)異常正常TN真陰性FP假陽性異常FN假陰性TP真陽性從混淆矩陣可以推導(dǎo)出一組指標(biāo)敏感度/召回率 TP / (TP FN)異常樣本被正確找出來的比例。漏掉一個(gè)異常樣本意味著漏診一個(gè)可能有心臟問題的病人這是最不能接受的錯(cuò)誤。特異性 TN / (TN FP)正常樣本被正確識(shí)別為正常的比例。誤報(bào)一個(gè)正常樣本會(huì)讓人白白焦慮但代價(jià)相對(duì)可控。F1分?jǐn)?shù) 2 × 精確率 × 召回率 / (精確率 召回率)精確率和召回率的調(diào)和平均兼顧兩方面的綜合評(píng)價(jià)。在我給社團(tuán)同學(xué)演示的過程中出現(xiàn)過這樣一幕某位同學(xué)的隨機(jī)森林模型準(zhǔn)確率89%看起來很亮眼但打開混淆矩陣發(fā)現(xiàn)異常類少數(shù)類的召回率只有62%也就是說十個(gè)真正異常的心音里它漏掉了將近四個(gè)。如果不是看了混淆矩陣這個(gè)嚴(yán)重缺陷根本不會(huì)被發(fā)現(xiàn)。后來我們調(diào)整了分類閾值默認(rèn)0.5改為0.4偏向于判斷為異常召回率提上去了整體F1分?jǐn)?shù)也確實(shí)變好了。5.2 交叉驗(yàn)證的正確姿勢(shì)防止“自欺欺人”的評(píng)估在最初跑通流程時(shí)我們犯過一個(gè)錯(cuò)誤把預(yù)處理后的全部特征矩陣隨機(jī)劃分成70%訓(xùn)練、30%測(cè)試一次訓(xùn)練一次測(cè)試拿到90%的準(zhǔn)確率就發(fā)朋友圈慶祝了。后來復(fù)盤發(fā)現(xiàn)一個(gè)問題同一個(gè)病人的好幾種錄音片段可能同時(shí)出現(xiàn)在訓(xùn)練集和測(cè)試集中模型相當(dāng)于是“見過這個(gè)人”之后去答題成績(jī)虛高不奇怪。正確的做法是按“錄音文件”分組劃分讓同一個(gè)病人的所有片段要么全部在訓(xùn)練集、要么全部在測(cè)試集。% 按錄音ID分組進(jìn)行分層交叉驗(yàn)證 patientIDs unique(fileIDs); cv cvpartition(patientIDs, KFold, 5); % 按病人做5折交叉驗(yàn)證還有一個(gè)很容易被忽視的操作細(xì)節(jié)特征標(biāo)準(zhǔn)化必須在每一折訓(xùn)練集上單獨(dú)進(jìn)行再用訓(xùn)練集的均值/標(biāo)準(zhǔn)差去轉(zhuǎn)換驗(yàn)證集和測(cè)試集。我在前文提到過這里的核心問題是數(shù)據(jù)泄漏。如果先把全量數(shù)據(jù)標(biāo)準(zhǔn)化再劃分那么驗(yàn)證集的信息就已經(jīng)“泄露”到訓(xùn)練集里了。更糟糕的是如果你用全量數(shù)據(jù)算均值和標(biāo)準(zhǔn)差在測(cè)試集上使用每次劃分均值標(biāo)準(zhǔn)差都不一樣這種不一致會(huì)讓你對(duì)結(jié)果的解讀越來越混亂。accFold zeros(5,1); for k 1:5 trainIdx training(cv, k); testIdx test(cv, k); % 注意只對(duì)訓(xùn)練集擬合標(biāo)準(zhǔn)化參數(shù) mu mean(featureMatrix(trainIdx,:), 1); sig std(featureMatrix(trainIdx,:), 0, 1); Xtr (featureMatrix(trainIdx,:) - mu) ./ sig; Xte (featureMatrix(testIdx,:) - mu) ./ sig; mdl fitcensemble(Xtr, labels(trainIdx), Method, AdaBoostM2); pred predict(mdl, Xte); accFold(k) sum(pred labels(testIdx)) / length(testIdx); end disp([5折交叉驗(yàn)證準(zhǔn)確率, num2str(mean(accFold))]);5.3 我踩過的幾個(gè)坑寫出來希望你別再踩第一個(gè)坑是原始波形直接喂模型。有同學(xué)想展示自己的“深度學(xué)習(xí)能力”把3秒時(shí)間窗的6000維特征直接扔給神經(jīng)網(wǎng)絡(luò)而不是先做特征提取。結(jié)果是訓(xùn)練集準(zhǔn)確率奇高、測(cè)試集準(zhǔn)確率慘不忍睹。原因很簡(jiǎn)單6000維中包含的絕大多數(shù)是噪聲和個(gè)體差異信息有效的結(jié)構(gòu)性信息占比極低模型學(xué)到的是“記樣本”而不是“學(xué)規(guī)律”。特征提取的價(jià)值就是把這6000維壓縮到36維、去除冗余、讓分類器的學(xué)習(xí)變得可能。第二個(gè)坑是帶通濾波頻率設(shè)置太保守。我們一開始參考了一些文獻(xiàn)把濾波范圍設(shè)成50Hz到300Hz結(jié)果S1心音里一個(gè)重要的低頻“thump”成分被濾掉了所有樣本的低頻形態(tài)被抹平特征區(qū)分度直接下降。后來對(duì)照原始波形的頻譜圖示波發(fā)現(xiàn)心音在25Hz到50Hz之間還有能量分布就把下界放寬到25Hz。這個(gè)改動(dòng)讓隨機(jī)森林的測(cè)試準(zhǔn)確率提升了3個(gè)百分點(diǎn)。第三個(gè)坑是混淆正常與異常標(biāo)簽的含義。PhysioNet 2016的原始標(biāo)簽里除了正常和異常之外其實(shí)還有第三類“不確定”。我們?cè)谡順?biāo)簽的時(shí)候一開始直接把“不確定”當(dāng)成正常處理導(dǎo)致一部分測(cè)試集上噪聲極大。后來仔細(xì)讀了數(shù)據(jù)集文檔把“不確定”類樣本直接丟棄用干凈的二分類數(shù)據(jù)訓(xùn)練指標(biāo)才穩(wěn)定下來。讀原始數(shù)據(jù)說明文檔這件事雖然枯燥但能省下后面無數(shù)的排查時(shí)間。第四個(gè)坑是版本兼容性。MATLAB從R2020a開始信號(hào)處理工具箱的某些函數(shù)簽名有變動(dòng)。我們的代碼里有同學(xué)用舊版本跑不了mfcc函數(shù)舊版本用的是melcepst或者去File Exchange找額外工具包。統(tǒng)一到R2022b之后才沒出幺蛾子。做項(xiàng)目之前先確認(rèn)團(tuán)隊(duì)軟件版本一致這句話再?gòu)?qiáng)調(diào)一遍都不過分。要說我做這個(gè)項(xiàng)目最大的收獲倒不是最終模型跑到了多少準(zhǔn)確率——畢竟心音分類本身是一個(gè)已經(jīng)有大量成熟研究的經(jīng)典課題我們做的這些遠(yuǎn)遠(yuǎn)談不上前沿。更值得說的是這套“信號(hào)清洗 → 人工特征設(shè)計(jì) → 基線模型 → 集成提升 → 嚴(yán)謹(jǐn)評(píng)估”的流程它像一個(gè)手藝活你不親手把數(shù)據(jù)讀進(jìn)來、畫出來、濾一遍、算特征、喂模型、調(diào)閾值你是很難理解什么叫“信號(hào)處理”和“機(jī)器學(xué)習(xí)”之間的銜接的。手里有了一套可以復(fù)用的代碼基底之后以后再遇到其他任何“信號(hào)分類”任務(wù)——軸承故障診斷、語音識(shí)別、腦電情緒分析——把這個(gè)流程換湯不換藥地搬過去就行模型結(jié)構(gòu)不用改、評(píng)估框架不用改最多改改濾波參數(shù)和特征組合。這種可遷移的能力才是這個(gè)入門項(xiàng)目真正想給你的東西。 MATLAB心音分類實(shí)戰(zhàn)：從信號(hào)預(yù)處理到特征提取的完整流程
看到“心音分類器”這個(gè)題目可能有朋友第一反應(yīng)是這不是典型的深度學(xué)習(xí)活兒?jiǎn)崮脗€(gè)預(yù)訓(xùn)練網(wǎng)絡(luò)懟上去精度刷到99%完事。但真當(dāng)你拿到一份真實(shí)的生理信號(hào)數(shù)據(jù)時(shí)會(huì)發(fā)現(xiàn)事情遠(yuǎn)沒有這么簡(jiǎn)單。信號(hào)里有環(huán)境噪聲、有呼吸音、有傳感器摩擦聲甚至還有醫(yī)生操作時(shí)不小心碰觸麥克風(fēng)的聲音直接丟給任何模型都是災(zāi)難。我去年在社團(tuán)周會(huì)上帶新人做機(jī)器學(xué)習(xí)入門項(xiàng)目選了一圈之后最后定在了“MATLAB 心音分類”這條線上折騰了大概三周最大的感受是用MATLAB做這個(gè)項(xiàng)目真正的價(jià)值點(diǎn)根本不在“分類器”那一環(huán)而是前面的信號(hào)預(yù)處理和特征提取——那才是決定模型上限的地方。這篇就把我們的完整流程、核心代碼和踩過的坑一次說清楚適合剛接觸機(jī)器學(xué)習(xí)、想在真實(shí)數(shù)據(jù)上完整跑通“信號(hào) → 特征 → 模型 → 評(píng)估”這條鏈路的朋友參考。1. 為什么是MATLAB為什么是心音信號(hào)1.1 這個(gè)項(xiàng)目在做什么把生物信號(hào)變成機(jī)器學(xué)習(xí)問題心音分類在最簡(jiǎn)單版本里是一個(gè)二分類問題給定一段心音錄音判斷它是“正?！边€是“異常”。異常心音包括各種雜音、瓣膜問題導(dǎo)致的額外音、心律不齊等等。醫(yī)學(xué)上醫(yī)生用聽診器聽診靠經(jīng)驗(yàn)判斷而我們要做的就是用機(jī)器學(xué)習(xí)模型去模擬這個(gè)判斷過程。選這個(gè)題作為機(jī)器學(xué)習(xí)初探項(xiàng)目有幾個(gè)天然優(yōu)勢(shì)。第一信號(hào)非常直觀——你把它畫出來正常人能看出明顯的“嘣噠”兩下心跳節(jié)律異常信號(hào)往往波形亂、幅值不均勻這個(gè)可視化過程對(duì)建立直覺特別有幫助。第二公開數(shù)據(jù)集很成熟PhysioNet 2016心音分類挑戰(zhàn)賽的數(shù)據(jù)集有三千多條錄音wav格式采樣率2000Hz文件組織和標(biāo)簽都齊。第三問題邊界清晰特征提取有明確的信號(hào)處理理論支撐分類器選擇可以從最基礎(chǔ)的決策樹一直玩到集成學(xué)習(xí)梯度非常平滑不會(huì)出現(xiàn)“要么完全不會(huì)、要么直接上深度學(xué)習(xí)黑盒”的斷層感。這個(gè)項(xiàng)目我們最終的目標(biāo)不是刷高精度數(shù)值而是讓社團(tuán)里的同學(xué)完整經(jīng)歷一遍機(jī)器學(xué)習(xí)的標(biāo)準(zhǔn)流程拿到原始數(shù)據(jù)、做清洗和預(yù)處理、人工設(shè)計(jì)特征、訓(xùn)練模型、用交叉驗(yàn)證評(píng)估泛化能力。走完這一步后面換任何領(lǐng)域的數(shù)據(jù)——比如電機(jī)振動(dòng)信號(hào)、語音命令、心電圖——套路都是通的。1.2 MATLAB在這條鏈路里的獨(dú)特優(yōu)勢(shì)我知道現(xiàn)在一說機(jī)器學(xué)習(xí)大家第一反應(yīng)是Python。但在這個(gè)項(xiàng)目里MATLAB其實(shí)是更有性價(jià)比的選擇尤其是在做信號(hào)處理和特征提取的階段。Python做信號(hào)處理不是不行但需要numpy、scipy、librosa、pandas、scikit-learn等等好幾個(gè)包配合著裝版本之間存在兼容性問題光是給新人配環(huán)境就能耗掉半天。而MATLAB一個(gè)軟件就把這些全干了信號(hào)處理工具箱里有濾波、傅里葉變換、小波分析、倒譜分析的完整函數(shù)統(tǒng)計(jì)與機(jī)器學(xué)習(xí)工具箱里有fitctree決策樹、fitcecoc多分類SVM、fitcknnKNN、fitcensemble集成學(xué)習(xí)等現(xiàn)成模型畫圖方面MATLAB的原生繪圖函數(shù)比matplotlib默認(rèn)樣式好看一個(gè)量級(jí)而且操作不需要寫一堆樣式調(diào)整代碼。另外還有一點(diǎn)在初學(xué)階段特別關(guān)鍵MATLAB的變量工作區(qū)是可視化的。你可以隨時(shí)雙擊變量看數(shù)據(jù)長(zhǎng)什么樣、尺寸是多少、是單精度還是雙精度這在調(diào)試特征矩陣維度的時(shí)侯簡(jiǎn)直救命。Python IDE里你還要print shape或者用debuggerMATLAB直接一目了然。MATLAB當(dāng)然也有它的缺點(diǎn)比如閉源、貴、深度學(xué)習(xí)的生態(tài)確實(shí)不如Python豐富。但對(duì)于“機(jī)器學(xué)習(xí)初探”這個(gè)定位來說這些缺點(diǎn)完全不影響。學(xué)生用校園授權(quán)版或者開源社區(qū)里跟Octave的兼容性方案都?jí)蛴?。我個(gè)人的態(tài)度是工具是手段不是信仰哪個(gè)能讓你最快把核心邏輯跑通就用哪個(gè)。這個(gè)項(xiàng)目的核心邏輯是信號(hào)處理和特征工程用MATLAB能讓你的精力百分之百集中在核心邏輯而不是環(huán)境配置上這就是我推薦它的理由。2. 心音數(shù)據(jù)的獲取與預(yù)處理先把“雜音”里的信息洗出來2.1 數(shù)據(jù)要從哪里來PhysioNet 2016數(shù)據(jù)集我們用的數(shù)據(jù)集是PhysioNet/Computing in Cardiology Challenge 2016提供的心音分類公開數(shù)據(jù)集這是一套專門用于心音自動(dòng)分析研究的基準(zhǔn)數(shù)據(jù)。數(shù)據(jù)分為training-a到training-e五個(gè)子文件夾一共包含三千多條心音錄音wav文件每個(gè)文件夾下面還有一個(gè)HEART_SOUND.csv文件里面是每條錄音對(duì)應(yīng)的標(biāo)簽1表示正常-1表示異常。先說幾個(gè)基本參數(shù)錄音統(tǒng)一是2000Hz采樣率也就是說每秒鐘有2000個(gè)采樣點(diǎn)人耳能聽到的頻率范圍是20Hz到20000Hz而心音的絕大部分能量集中在20Hz到400Hz這個(gè)低頻段2000Hz采樣率對(duì)心音來說信息量完全足夠。讀取文件用MATLAB的audioread函數(shù)即可。% 讀取單條心音錄音 [y, fs] audioread(training-a/a0001.wav); % y是采樣點(diǎn)數(shù)組fs是采樣率(2000Hz) disp(size(y)); disp(fs);這里有一個(gè)新手常犯的錯(cuò)誤錄音文件不是等長(zhǎng)的。有些錄音3秒有些8秒文件總時(shí)長(zhǎng)差異很大。如果直接把整段錄音作為輸入不同樣本的特征維度就對(duì)不上。所以預(yù)處理階段要統(tǒng)一處理到固定長(zhǎng)度或者做分段切割后面我會(huì)細(xì)說。2.2 預(yù)處理三步走降采樣、濾波、歸一化拿到原始錄音之后不要急著提特征。原始信號(hào)里包含了大量與心音無關(guān)的“臟東西”我把它稱為噪聲三件套環(huán)境噪聲空調(diào)聲、房間混響、生理干擾呼吸聲、肌肉信號(hào)、傳感器噪聲麥克風(fēng)接觸摩擦、電噪聲。此外還有50Hz/60Hz的工頻干擾。我用的預(yù)處理流程是三步降采樣可選、帶通濾波、幅值歸一化。降采樣既然心音主要能量在400Hz以下根據(jù)奈奎斯特采樣定理理論上800Hz采樣率就足夠重建信號(hào)。把2000Hz降采樣到1000Hz甚至800Hz好處是數(shù)據(jù)量減半后續(xù)計(jì)算更快。不過對(duì)入門項(xiàng)目來說2000Hz的數(shù)據(jù)量本身也不大降不降采樣對(duì)最終結(jié)果影響很小所以這一步我們后來基本跳過。帶通濾波這是我們預(yù)處理里的核心步驟。心音的診斷信息主要位于25Hz到400Hz之間因此設(shè)計(jì)了一個(gè)四階巴特沃斯帶通濾波器保留這個(gè)頻段濾掉低頻漂移和高頻噪聲。低頻漂移特別值得注意因?yàn)閭鞲衅鹘佑|或呼吸動(dòng)作會(huì)產(chǎn)生一種非常慢的基線偏移如果不濾掉后面提取的時(shí)域特征比如均值和方差會(huì)被這個(gè)漂移帶偏。% 設(shè)計(jì)四階巴特沃斯帶通濾波器通帶25Hz~400Hz fs 2000; fL 25; fH 400; [b, a] butter(4, [fL fH]/(fs/2), bandpass); % 使用filtfilt零相位濾波避免相位失真 y_filtered filtfilt(b, a, y);為什么用filtfilt而不用filter這是很多初學(xué)會(huì)忽略的細(xì)節(jié)。filter是常規(guī)濾波會(huì)帶來相位偏移filtfilt先把信號(hào)正向?yàn)V一遍、再反向?yàn)V一遍兩次的相位偏移互相抵消得到零相位失真的結(jié)果。對(duì)心音這種看重波形形態(tài)的分析來說絕對(duì)有必要。幅值歸一化因?yàn)椴煌浺舻匿浿圃O(shè)備、深淺閾值、病人體型都不同音量幅值差異非常大。歸一化就是把每個(gè)樣本的最大幅值縮放為1。這樣后續(xù)特征比如峰值、能量才有跨樣本可比性。% 將幅值歸一化到[-1, 1] y_normalized y_filtered / max(abs(y_filtered));我見過一個(gè)很隱蔽的坑如果你先歸一化、再濾波濾波后的信號(hào)幅值可能又超過1了因?yàn)闉V波本身會(huì)改變信號(hào)的能量分布。所以順序應(yīng)該是先濾波、再歸一化。這個(gè)小細(xì)節(jié)不影響大的方向但會(huì)讓你的特征數(shù)值更干凈。2.3 信號(hào)分段與樣本構(gòu)建從連續(xù)波形到一個(gè)個(gè)“樣本”現(xiàn)在手里是清洗過的完整錄音。但問題來了一條錄音少說也有好幾秒心音是典型的周期信號(hào)——正常心跳每分鐘60到100次大約每0.6秒到1秒一個(gè)周期。整段錄音當(dāng)成一個(gè)樣本信息太稠、特征糊成一片分類器根本找不到判別性信息。更合理的做法是分段。我們把每條錄音切成了固定長(zhǎng)度為3秒的片段當(dāng)作獨(dú)立的樣本參與后續(xù)訓(xùn)練。為什么是3秒經(jīng)驗(yàn)值是3秒至少包含3個(gè)心動(dòng)周期既能捕捉到周期性和一致性又不至于因?yàn)榘芷谔喽旬惓＜?xì)節(jié)平均掉。% 將信號(hào)切分為3秒一段 segLen 3 * fs; % 3秒 numSeg floor(length(y_normalized) / segLen); for i 1:numSeg seg y_normalized((i-1)*segLen 1 : i*segLen); % 將seg作為一條樣本存入后續(xù)的特征提取環(huán)節(jié) end這里有個(gè)取舍需要說清楚。理想情況下我們應(yīng)該是先做心音分割識(shí)別出S1、收縮期、S2、舒張期這四段再?gòu)拿總€(gè)心動(dòng)周期提特征這樣特征和生理機(jī)制對(duì)應(yīng)得更準(zhǔn)。但心音分割本身就是一個(gè)不小的課題需要做包絡(luò)檢測(cè)、峰值定位、邏輯判斷相當(dāng)于在分類問題之前又疊了一個(gè)完整項(xiàng)目對(duì)初學(xué)者來說步子邁得太大。我們最終采取的方案是“簡(jiǎn)單粗暴按時(shí)間窗切”訓(xùn)練出來的模型效果其實(shí)已經(jīng)夠用了。初探階段先跑通流程比什么都重要。3. 特征提取把心跳聲翻譯成機(jī)器認(rèn)識(shí)的數(shù)字3.1 時(shí)域特征如何描述一段波形的“性格”信號(hào)本身是幾千個(gè)采樣點(diǎn)機(jī)器學(xué)習(xí)模型不能直接吃這種原始波形。原始波形就算長(zhǎng)度統(tǒng)一維度也太大3秒×2000采樣率6000維而且每個(gè)采樣點(diǎn)的具體數(shù)值受噪聲影響很大不夠魯棒。特征提取的目標(biāo)就是把一段波形壓縮成一組維度不高、但能描述其本質(zhì)屬性的數(shù)字。我管時(shí)域特征叫“波形性格參數(shù)”它們是直接在時(shí)間軸上計(jì)算出的統(tǒng)計(jì)量均方根值RMS描述信號(hào)整體能量大小。心音異常時(shí)往往伴隨額外雜音信號(hào)能量分布會(huì)變化。方差描述信號(hào)波動(dòng)程度。異常心音經(jīng)常更“亂”方差會(huì)變大。峰值/谷值描述單次心跳的最大沖擊力。S1第一心音和S2第二心音的峰值比例在不同疾病下會(huì)有變化。過零率信號(hào)在單位時(shí)間內(nèi)穿越零點(diǎn)的次數(shù)。過零率越高說明信號(hào)越“高頻”異常雜音通常是高頻的會(huì)讓過零率顯著上升。峰值因子峰值除以RMS描述波形“尖銳”程度。心音異常時(shí)心臟射血狀態(tài)改變波形形態(tài)隨之改變。function feat extractTimeFeatures(seg) feat zeros(1, 5); feat(1) rms(seg); % 均方根值 feat(2) var(seg); % 方差 feat(3) max(abs(seg)); % 峰值 feat(4) sum(abs(diff(sign(seg)))) / (length(seg)-1); % 過零率 feat(5) max(abs(seg)) / (rms(seg) eps); % 峰值因子 end這幾個(gè)人人都能算、看起來“沒什么技術(shù)含量”的特征在實(shí)際分類中的貢獻(xiàn)往往出乎意料地大。尤其是過零率和方差訓(xùn)練出來的隨機(jī)森林模型里重要性排得都很靠前。原因也簡(jiǎn)單異常心音中高頻雜音成分多信號(hào)振蕩更頻繁過零率自然高。你不需要理解它的生理機(jī)理只需要信任它能區(qū)分正常和異常。3.2 頻域特征與倒譜MFCC為什么能跨界單獨(dú)靠時(shí)域特征還不夠——不同病人心音的時(shí)域波形即使同是“正?！毙螒B(tài)差異也很大直接用原始幅值統(tǒng)計(jì)量做特征很難做到跨個(gè)體泛化。頻域特征的加入相當(dāng)于從另一個(gè)角度描述信號(hào)不再是“這個(gè)波形長(zhǎng)什么樣”而是“這個(gè)波形由哪些頻率成分構(gòu)成”。首先做短時(shí)傅里葉變換得到語譜圖?；谡Z譜圖可以計(jì)算頻域特征頻譜質(zhì)心信號(hào)能量在頻率上的加權(quán)平均描述“聲音聽起來亮不亮”。頻譜帶寬能量在質(zhì)心周圍的散布程度異常雜音會(huì)讓帶寬變大。頻譜滾降點(diǎn)信號(hào)能量有多少集中在低頻端的閾值點(diǎn)心音是低頻為主滾降點(diǎn)通常在低頻。然后是MFCC梅爾頻率倒譜系數(shù)。這個(gè)詞大家可能在語音識(shí)別里聽過但心音也是聲學(xué)信號(hào)MFCC同樣適用。MFCC把頻譜按人耳感知的梅爾刻度映射到對(duì)數(shù)刻度然后做離散余弦變換得到的倒譜系數(shù)能非常緊湊地表示頻譜包絡(luò)的形狀。心音的頻譜包絡(luò)跟心臟結(jié)構(gòu)、瓣膜狀態(tài)有關(guān)異常狀態(tài)下包絡(luò)形狀變化明顯MFCC能更好捕捉這種變化。% 計(jì)算13維MFCC特征 [coeff, ~] mfcc(seg, fs, NumCoeffs, 13, ... WindowLength, round(0.03*fs), OverlapLength, round(0.02*fs)); % 沿時(shí)間軸取平均值和標(biāo)準(zhǔn)差得到26維特征 mfccMean mean(coeff, 1); mfccStd std(coeff, 0, 1); featMFCC [mfccMean, mfccStd];MFCC的重要性在于它經(jīng)過語音識(shí)別領(lǐng)域幾十年驗(yàn)證是一種對(duì)“聲音”來說極其高效的壓縮表示。心音本質(zhì)上也是聲音共享這套表示完全合理。用在我們的模型里加入MFCC后交叉驗(yàn)證準(zhǔn)確率大約提升了4到5個(gè)百分點(diǎn)是貢獻(xiàn)最大的一類特征。3.3 用MATLAB組裝特征矩陣別在維度上翻車把所有片段的特征拼成一個(gè)矩陣這是機(jī)器學(xué)習(xí)訓(xùn)練的前置工作也是很多新手第一次接觸“特征矩陣”這個(gè)概念的地方。這里有一個(gè)非常重要的規(guī)范每個(gè)樣本提取的特征必須是相同長(zhǎng)度的向量。比如時(shí)域特征選了5個(gè)頻域特征選了5個(gè)MFCC選了26個(gè)那這個(gè)樣本的特征向量就是1×36。所有樣本都提這36個(gè)特征最后拼出一個(gè)N×36的特征矩陣。featureMatrix []; labels []; for i 1:length(fileList) [y, fs] audioread(fileList{i}); % 預(yù)處理濾波、歸一化 y filtfilt(b, a, y); y y / max(abs(y)); % 分段 segLen 3 * fs; numSeg floor(length(y) / segLen); for j 1:numSeg seg y((j-1)*segLen1 : j*segLen); featT extractTimeFeatures(seg); % 5維 featF extractFreqFeatures(seg, fs); % 5維 featM extractMFCCFeatures(seg, fs); % 26維 featAll [featT, featF, featM]; % 36維 featureMatrix [featureMatrix; featAll]; labels [labels; labelVal]; % 0或1 end end % 保存到mat文件供后續(xù)訓(xùn)練復(fù)用 save(heart_feature.mat, featureMatrix, labels, -v7.3);這段代碼里我特意標(biāo)了一個(gè)重要的壞習(xí)慣示例特征矩陣和標(biāo)簽用[matrix; newRow]這種橫向拼接方式。在小數(shù)據(jù)量下它能工作但樣本上萬的時(shí)候這種寫法會(huì)反復(fù)復(fù)制內(nèi)存慢得讓人抓狂。正確做法是預(yù)分配矩陣或者把每個(gè)樣本的結(jié)果存到cell里最后一次性合并。被教育過一次之后我現(xiàn)在寫的代碼都是先zeros(N, 36)初始化好再填充。這算是個(gè)性能層面的教訓(xùn)。組裝好之后建議立刻做一步特征標(biāo)準(zhǔn)化。什么是標(biāo)準(zhǔn)化就是把每個(gè)特征的均值變?yōu)?、標(biāo)準(zhǔn)差變?yōu)?。featureMatrixNorm normalize(featureMatrix);為什么標(biāo)準(zhǔn)化是必須的因?yàn)椴煌卣髁烤V差異巨大——MFCC的數(shù)值可能是幾十而RMS只有零點(diǎn)幾。很多分類器內(nèi)部計(jì)算距離或正則化項(xiàng)會(huì)讓量綱大的特征主導(dǎo)結(jié)果。標(biāo)準(zhǔn)化就是消除量綱影響。需要提醒的是標(biāo)準(zhǔn)化參數(shù)只能在訓(xùn)練集上計(jì)算測(cè)試集要用同樣的均值和標(biāo)準(zhǔn)差去轉(zhuǎn)換這個(gè)細(xì)節(jié)放到第5節(jié)交叉驗(yàn)證的時(shí)候再具體展開。4. 分類器選型與訓(xùn)練從基線模型開始的進(jìn)階路線4.1 為什么不要一上來就上深度學(xué)習(xí)把特征矩陣準(zhǔn)備好之后接下來就是模型訓(xùn)練。但這里我要反復(fù)強(qiáng)調(diào)一個(gè)初學(xué)階段的觀念不要一上來就上深度學(xué)習(xí)。神經(jīng)網(wǎng)絡(luò)是一個(gè)高度非線性的“萬能函數(shù)”但它需要海量數(shù)據(jù)、大量調(diào)參、以及在驗(yàn)證集上的迭代監(jiān)控。我們的心音數(shù)據(jù)集經(jīng)過分段后樣本量大約幾千條特征維度只有36維這種“小樣本中等維度”的設(shè)定本身就是傳統(tǒng)機(jī)器學(xué)習(xí)的主場(chǎng)。用傳統(tǒng)機(jī)器學(xué)習(xí)有幾個(gè)看得見摸得著的好處第一訓(xùn)練快幾十秒出結(jié)果可以反復(fù)試第二可解釋性強(qiáng)比如決策樹你可以把樹畫出來看它到底學(xué)了什么規(guī)則第三對(duì)特征工程的好壞能立刻得到反饋。你加了一組特征、準(zhǔn)確率提升了你知道是這組特征起了作用你換了個(gè)亂七八糟的特征、準(zhǔn)確率下降了你也知道是這組特征在拖后腿。這種反饋對(duì)建立直覺極其重要是深度學(xué)習(xí)的黑盒給不了的。我在表里列了四種適合初學(xué)的分類器以及它們?cè)诒卷?xiàng)目上的表現(xiàn)對(duì)比分類器原理一句話本項(xiàng)目初測(cè)準(zhǔn)確率訓(xùn)練速度調(diào)參難度決策樹用一系列“特征閾值”的判斷做分類約78%極快低K近鄰KNN看距離最近的K個(gè)樣本屬于哪類約82%訓(xùn)練快、預(yù)測(cè)慢低支持向量機(jī)SVM在高維空間找一個(gè)最優(yōu)分割超平面約84%快中隨機(jī)森林訓(xùn)練多棵決策樹投票決定結(jié)果約88%快中低AdaBoost集成用多個(gè)弱分類器加權(quán)組合成強(qiáng)分類器約90%中等中從這個(gè)表能明顯看出規(guī)律集成方法隨機(jī)森林、AdaBoost比單模型表現(xiàn)好。原因也很直觀心音特征內(nèi)部往往存在復(fù)雜的交互關(guān)系比如某個(gè)時(shí)域特征在結(jié)合某個(gè)頻域特征的時(shí)候才有判別意義單棵決策樹對(duì)這個(gè)交互的建模能力是有限的但多棵樹從不同角度切分?jǐn)?shù)據(jù)再投票能把這種交互關(guān)系綜合起來。4.2 三種主流分類器的MATLAB實(shí)現(xiàn)與調(diào)用MATLAB的fitc系列函數(shù)把傳統(tǒng)機(jī)器學(xué)習(xí)封裝得極其友好。下面給出我們項(xiàng)目里用到的決策樹、SVM和集成學(xué)習(xí)三種實(shí)現(xiàn)的對(duì)照代碼讀者拿去改數(shù)據(jù)路徑就能跑通。% 假設(shè)featureNorm是N×36特征矩陣label是N×1標(biāo)簽0或1 % 1. 決策樹 mdlTree fitctree(featureNorm, label, CrossVal, off); accuracyTree 1 - resubLoss(mdlTree); % 注意這只是訓(xùn)練集回代精度 % 2. SVMRBF核 mdlSVM fitcecoc(featureNorm, label, Learners, svm); % fitcecoc是錯(cuò)誤糾錯(cuò)輸出碼用于多分類二分類時(shí)等價(jià)于一個(gè)SVM % 3. 集成學(xué)習(xí)AdaBoostM1 決策樹弱學(xué)習(xí)器 t templateTree(MaxNumSplits, 10); mdlEnsemble fitcensemble(featureNorm, label, ... Method, AdaBoostM2, ... NumLearningCycles, 50, ... Learners, t);fitcensemble是這里面的明星函數(shù)。它內(nèi)部自動(dòng)處理了自適應(yīng)增強(qiáng)的各個(gè)細(xì)節(jié)你只需要指定弱學(xué)習(xí)器和方法。我們指定MaxNumSplits10是為了限制每棵弱樹不要太深太深的樹容易過擬合保留一些“弱能力”的模型反而讓Adaboost更有效。訓(xùn)練完之后預(yù)測(cè)新數(shù)據(jù)的邏輯更簡(jiǎn)單% 對(duì)測(cè)試集預(yù)測(cè) predLabel predict(mdlEnsemble, testFeatureNorm); % 計(jì)算準(zhǔn)確率 acc sum(predLabel testLabel) / length(testLabel);4.3 集成模型的參數(shù)調(diào)節(jié)經(jīng)驗(yàn)機(jī)器學(xué)習(xí)里有個(gè)概念叫“偏差-方差權(quán)衡”模型太簡(jiǎn)單會(huì)欠擬合偏差高模型太復(fù)雜會(huì)過擬合方差高。集成模型的設(shè)計(jì)初衷就是通過組合多個(gè)弱模型在保持低偏差的同時(shí)降低方差。用fitcensemble的時(shí)候有兩個(gè)參數(shù)最值得調(diào)NumLearningCycles弱學(xué)習(xí)器數(shù)量和Learners里的深度限制。我自己的經(jīng)驗(yàn)是弱學(xué)習(xí)器數(shù)量從10加到50時(shí)效果提升最明顯到100以上提升就趨于平緩了。收益邊際遞減這個(gè)現(xiàn)象在集成學(xué)習(xí)里非常典型。弱樹的深度控制在5到10層比較合適。太深容易單棵樹就過擬合了太淺每棵樹都太笨集成出來的模型不夠聰明。在深度學(xué)習(xí)特征的標(biāo)準(zhǔn)化問題上對(duì)樹模型來說標(biāo)準(zhǔn)化影響不大因?yàn)闃渲魂P(guān)心閾值切分不關(guān)心距離但對(duì)SVM和KNN這種基于距離的模型影響巨大。我們的做法是找出全部特征直接標(biāo)準(zhǔn)化一視同仁。調(diào)參的時(shí)候不要靠猜可以用MATLAB的fitcensemble配合OptimizeHyperparameters做自動(dòng)超參數(shù)搜索。不過要注意超參數(shù)搜索會(huì)跑很久數(shù)據(jù)量大的時(shí)候建議先抽取一部分?jǐn)?shù)據(jù)做快速實(shí)驗(yàn)確定參數(shù)范圍后再全量訓(xùn)練。% 自動(dòng)超參數(shù)調(diào)優(yōu)示例小心運(yùn)行時(shí)間 mdlOpt fitcensemble(featureNorm, label, ... Method, AdaBoostM2, ... OptimizeHyperparameters, {NumLearningCycles, LearnRate}, ... HyperparameterOptimizationOptions, struct(MaxObjectiveEvaluations, 30));5. 評(píng)估與調(diào)參準(zhǔn)確率之外還要看敏感度和特異性5.1 混淆矩陣與三類指標(biāo)別被準(zhǔn)確率蒙蔽訓(xùn)練完模型很多人只會(huì)看一個(gè)準(zhǔn)確率數(shù)字。但在這個(gè)項(xiàng)目里我要強(qiáng)調(diào)準(zhǔn)確率是一個(gè)不夠用的指標(biāo)。原因在于正負(fù)樣本不均衡心音數(shù)據(jù)集中正常樣本比異常樣本略多即使你無腦把全部樣本輸出成“正常”準(zhǔn)確率也能到差不多52%。在這種背景下準(zhǔn)確率不能描述模型在“少數(shù)類”上的表現(xiàn)。我們需要看混淆矩陣它把預(yù)測(cè)結(jié)果和真實(shí)標(biāo)簽的對(duì)應(yīng)關(guān)系展開成一個(gè)2×2表真實(shí)\預(yù)測(cè)預(yù)測(cè)正常預(yù)測(cè)異常正常TN真陰性FP假陽性異常FN假陰性TP真陽性從混淆矩陣可以推導(dǎo)出一組指標(biāo)敏感度/召回率 TP / (TP FN)異常樣本被正確找出來的比例。漏掉一個(gè)異常樣本意味著漏診一個(gè)可能有心臟問題的病人這是最不能接受的錯(cuò)誤。特異性 TN / (TN FP)正常樣本被正確識(shí)別為正常的比例。誤報(bào)一個(gè)正常樣本會(huì)讓人白白焦慮但代價(jià)相對(duì)可控。F1分?jǐn)?shù) 2 × 精確率 × 召回率 / (精確率 召回率)精確率和召回率的調(diào)和平均兼顧兩方面的綜合評(píng)價(jià)。在我給社團(tuán)同學(xué)演示的過程中出現(xiàn)過這樣一幕某位同學(xué)的隨機(jī)森林模型準(zhǔn)確率89%看起來很亮眼但打開混淆矩陣發(fā)現(xiàn)異常類少數(shù)類的召回率只有62%也就是說十個(gè)真正異常的心音里它漏掉了將近四個(gè)。如果不是看了混淆矩陣這個(gè)嚴(yán)重缺陷根本不會(huì)被發(fā)現(xiàn)。后來我們調(diào)整了分類閾值默認(rèn)0.5改為0.4偏向于判斷為異常召回率提上去了整體F1分?jǐn)?shù)也確實(shí)變好了。5.2 交叉驗(yàn)證的正確姿勢(shì)防止“自欺欺人”的評(píng)估在最初跑通流程時(shí)我們犯過一個(gè)錯(cuò)誤把預(yù)處理后的全部特征矩陣隨機(jī)劃分成70%訓(xùn)練、30%測(cè)試一次訓(xùn)練一次測(cè)試拿到90%的準(zhǔn)確率就發(fā)朋友圈慶祝了。后來復(fù)盤發(fā)現(xiàn)一個(gè)問題同一個(gè)病人的好幾種錄音片段可能同時(shí)出現(xiàn)在訓(xùn)練集和測(cè)試集中模型相當(dāng)于是“見過這個(gè)人”之后去答題成績(jī)虛高不奇怪。正確的做法是按“錄音文件”分組劃分讓同一個(gè)病人的所有片段要么全部在訓(xùn)練集、要么全部在測(cè)試集。% 按錄音ID分組進(jìn)行分層交叉驗(yàn)證 patientIDs unique(fileIDs); cv cvpartition(patientIDs, KFold, 5); % 按病人做5折交叉驗(yàn)證還有一個(gè)很容易被忽視的操作細(xì)節(jié)特征標(biāo)準(zhǔn)化必須在每一折訓(xùn)練集上單獨(dú)進(jìn)行再用訓(xùn)練集的均值/標(biāo)準(zhǔn)差去轉(zhuǎn)換驗(yàn)證集和測(cè)試集。我在前文提到過這里的核心問題是數(shù)據(jù)泄漏。如果先把全量數(shù)據(jù)標(biāo)準(zhǔn)化再劃分那么驗(yàn)證集的信息就已經(jīng)“泄露”到訓(xùn)練集里了。更糟糕的是如果你用全量數(shù)據(jù)算均值和標(biāo)準(zhǔn)差在測(cè)試集上使用每次劃分均值標(biāo)準(zhǔn)差都不一樣這種不一致會(huì)讓你對(duì)結(jié)果的解讀越來越混亂。accFold zeros(5,1); for k 1:5 trainIdx training(cv, k); testIdx test(cv, k); % 注意只對(duì)訓(xùn)練集擬合標(biāo)準(zhǔn)化參數(shù) mu mean(featureMatrix(trainIdx,:), 1); sig std(featureMatrix(trainIdx,:), 0, 1); Xtr (featureMatrix(trainIdx,:) - mu) ./ sig; Xte (featureMatrix(testIdx,:) - mu) ./ sig; mdl fitcensemble(Xtr, labels(trainIdx), Method, AdaBoostM2); pred predict(mdl, Xte); accFold(k) sum(pred labels(testIdx)) / length(testIdx); end disp([5折交叉驗(yàn)證準(zhǔn)確率, num2str(mean(accFold))]);5.3 我踩過的幾個(gè)坑寫出來希望你別再踩第一個(gè)坑是原始波形直接喂模型。有同學(xué)想展示自己的“深度學(xué)習(xí)能力”把3秒時(shí)間窗的6000維特征直接扔給神經(jīng)網(wǎng)絡(luò)而不是先做特征提取。結(jié)果是訓(xùn)練集準(zhǔn)確率奇高、測(cè)試集準(zhǔn)確率慘不忍睹。原因很簡(jiǎn)單6000維中包含的絕大多數(shù)是噪聲和個(gè)體差異信息有效的結(jié)構(gòu)性信息占比極低模型學(xué)到的是“記樣本”而不是“學(xué)規(guī)律”。特征提取的價(jià)值就是把這6000維壓縮到36維、去除冗余、讓分類器的學(xué)習(xí)變得可能。第二個(gè)坑是帶通濾波頻率設(shè)置太保守。我們一開始參考了一些文獻(xiàn)把濾波范圍設(shè)成50Hz到300Hz結(jié)果S1心音里一個(gè)重要的低頻“thump”成分被濾掉了所有樣本的低頻形態(tài)被抹平特征區(qū)分度直接下降。后來對(duì)照原始波形的頻譜圖示波發(fā)現(xiàn)心音在25Hz到50Hz之間還有能量分布就把下界放寬到25Hz。這個(gè)改動(dòng)讓隨機(jī)森林的測(cè)試準(zhǔn)確率提升了3個(gè)百分點(diǎn)。第三個(gè)坑是混淆正常與異常標(biāo)簽的含義。PhysioNet 2016的原始標(biāo)簽里除了正常和異常之外其實(shí)還有第三類“不確定”。我們?cè)谡順?biāo)簽的時(shí)候一開始直接把“不確定”當(dāng)成正常處理導(dǎo)致一部分測(cè)試集上噪聲極大。后來仔細(xì)讀了數(shù)據(jù)集文檔把“不確定”類樣本直接丟棄用干凈的二分類數(shù)據(jù)訓(xùn)練指標(biāo)才穩(wěn)定下來。讀原始數(shù)據(jù)說明文檔這件事雖然枯燥但能省下后面無數(shù)的排查時(shí)間。第四個(gè)坑是版本兼容性。MATLAB從R2020a開始信號(hào)處理工具箱的某些函數(shù)簽名有變動(dòng)。我們的代碼里有同學(xué)用舊版本跑不了mfcc函數(shù)舊版本用的是melcepst或者去File Exchange找額外工具包。統(tǒng)一到R2022b之后才沒出幺蛾子。做項(xiàng)目之前先確認(rèn)團(tuán)隊(duì)軟件版本一致這句話再?gòu)?qiáng)調(diào)一遍都不過分。要說我做這個(gè)項(xiàng)目最大的收獲倒不是最終模型跑到了多少準(zhǔn)確率——畢竟心音分類本身是一個(gè)已經(jīng)有大量成熟研究的經(jīng)典課題我們做的這些遠(yuǎn)遠(yuǎn)談不上前沿。更值得說的是這套“信號(hào)清洗 → 人工特征設(shè)計(jì) → 基線模型 → 集成提升 → 嚴(yán)謹(jǐn)評(píng)估”的流程它像一個(gè)手藝活你不親手把數(shù)據(jù)讀進(jìn)來、畫出來、濾一遍、算特征、喂模型、調(diào)閾值你是很難理解什么叫“信號(hào)處理”和“機(jī)器學(xué)習(xí)”之間的銜接的。手里有了一套可以復(fù)用的代碼基底之后以后再遇到其他任何“信號(hào)分類”任務(wù)——軸承故障診斷、語音識(shí)別、腦電情緒分析——把這個(gè)流程換湯不換藥地搬過去就行模型結(jié)構(gòu)不用改、評(píng)估框架不用改最多改改濾波參數(shù)和特征組合。這種可遷移的能力才是這個(gè)入門項(xiàng)目真正想給你的東西。