簡(jiǎn)介面向從事信號(hào)處理、數(shù)據(jù)分析的科研人員與工程師這份資源系統(tǒng)講解小波分析在MATLAB中的落地方法重點(diǎn)解決非平穩(wěn)信號(hào)的時(shí)頻分析、去噪與特征提取問(wèn)題。內(nèi)容從傅里葉變換與小波的對(duì)比入手結(jié)合db4基函數(shù)、多尺度分解、閾值去噪、信號(hào)重構(gòu)與信噪比評(píng)估并討論邊界效應(yīng)和分解層數(shù)選擇等常見難點(diǎn)附有大量可直接練習(xí)的MATLAB代碼示例。資源包內(nèi)含1個(gè)PDF文檔大小約109KB目前已有105人瀏覽學(xué)習(xí)。閱讀后可快速掌握心電圖、腦電圖、機(jī)械振動(dòng)等信號(hào)的處理思路也能根據(jù)場(chǎng)景合理選擇小波基與分解層數(shù)提升實(shí)際項(xiàng)目中的分析效率。1. 小波分析是信號(hào)處理的動(dòng)態(tài) CT從心電圖和振動(dòng)信號(hào)說(shuō)起做信號(hào)處理久了你會(huì)發(fā)現(xiàn)一件事傅里葉變換在平穩(wěn)信號(hào)上有多香在非平穩(wěn)信號(hào)上就有多容易翻車。一段心電信號(hào)呼吸引起的基線漂移疊著 50Hz 工頻干擾QRS 波群的頻率特征還隨時(shí)間變化——這時(shí)候直接上 FFT頻域里什么都看不清時(shí)域里噪聲又蓋住了關(guān)鍵波形。小波分解不一樣它相當(dāng)于給信號(hào)做了一次動(dòng)態(tài) CT 掃描同時(shí)保留時(shí)間和頻率兩個(gè)維度的信息。這份資源講的就是怎么在 MATLAB 里用小波分析做去噪、特征提取和故障識(shí)別選什么基函數(shù)、分解幾層、閾值怎么定還有我實(shí)際調(diào)試時(shí)踩過(guò)的坑。適合手里有非平穩(wěn)信號(hào)、想盡快落地的工程師和研究人員。2. 基函數(shù)與多尺度分解db4、symlets 怎么選wavedec 返回了什么2.1 為什么是 db4小波基的形態(tài)匹配邏輯與三種信號(hào)選型小波分解的第一步不是寫代碼是選基函數(shù)。這步選錯(cuò)了后面所有系數(shù)解讀都得推翻重來(lái)。常見做法是按信號(hào)形態(tài)去匹配Daubechies 系列在生物醫(yī)學(xué)信號(hào)里用得最多尤其是 db4它的尺度函數(shù)形狀和心電信號(hào)的 QRS 波群有一定相似度分解出來(lái)的細(xì)節(jié)系數(shù)能較好地保留波形尖銳特征。機(jī)械振動(dòng)信號(hào)里 symlets 系列表現(xiàn)更好sym5 比 db4 的相位更線性重構(gòu)波形平滑度更高適合找軸承固有頻率這類周期性成分。腦電信號(hào)我遇到過(guò) db4 把 alpha 波切得支離破碎的情況換成 bior3.5 才解決——雙正交小波在非線性相位場(chǎng)景下有優(yōu)勢(shì)但計(jì)算復(fù)雜度會(huì)高一些。選型時(shí)可以寫個(gè)簡(jiǎn)單腳本對(duì)比% 加載一段含噪聲的實(shí)測(cè)振動(dòng)信號(hào)采樣率1000Hz時(shí)長(zhǎng)5秒 load noissin; signal noissin(1:5000); t linspace(0, 5, 5000); [c1, l1] wavedec(signal, 5, db4); [c2, l2] wavedec(signal, 5, sym5); subplot(211); plot(wrcoef(a, c1, l1, db4, 5)); title(db4 低頻重構(gòu)); subplot(212); plot(wrcoef(a, c2, l2, sym5, 5)); title(sym5 低頻重構(gòu));這段代碼的邏輯是控制變量——分解層數(shù)、信號(hào)、重構(gòu)層級(jí)都一致只換基函數(shù)對(duì)比低頻逼近的形態(tài)差異。我一般會(huì)同時(shí)看兩個(gè)指標(biāo)重構(gòu)波形的平滑度以及是否保留原始信號(hào)中的瞬態(tài)突變。平滑度高的適合做趨勢(shì)提取保留細(xì)節(jié)多的適合做特征定位。參數(shù)說(shuō)明wavedec的第一個(gè)參數(shù)是信號(hào)向量第二個(gè)是分解層數(shù)第三個(gè)是基函數(shù)名。實(shí)際項(xiàng)目里我會(huì)把信號(hào)換成自己的 CSV 數(shù)據(jù)用importdata或readmatrix讀進(jìn)來(lái)波形形態(tài)不同小波基的選擇結(jié)論也會(huì)變。2.2 wavedec 的返回結(jié)構(gòu)c 和 l 到底怎么讀新手最容易搞混的是wavedec返回的c和l。c不是按層數(shù)排的它是按照從最后一層低頻系數(shù)到各層高頻系數(shù)的順序拼接的一維數(shù)組。l則是位置向量記錄每一段在c里的起始位置和長(zhǎng)度。比如三層分解c的結(jié)構(gòu)是第三層近似系數(shù)段、第三層細(xì)節(jié)系數(shù)段、第二層細(xì)節(jié)系數(shù)段、第一層細(xì)節(jié)系數(shù)段。l有 5 個(gè)元素分別對(duì)應(yīng)原始信號(hào)長(zhǎng)度、cA3 長(zhǎng)度、cD3 長(zhǎng)度、cD2 長(zhǎng)度、cD1 長(zhǎng)度。% 三層小波分解基函數(shù)用 db4 [c, l] wavedec(ecg, 3, db4); % 提取各層系數(shù) approx appcoef(c, l, db4); % 第三層近似系數(shù) [cd1, cd2, cd3] detcoef(c, l, [1 2 3]); % 一、二、三層細(xì)節(jié)系數(shù) % 查看各層長(zhǎng)度 disp(l);這段代碼里appcoef默認(rèn)返回最后一層的近似系數(shù)如果只寫appcoef(c, l, db4)等價(jià)于appcoef(c, l, db4, 3)。detcoef的第三個(gè)參數(shù)傳入向量[1 2 3]一次性取出三層細(xì)節(jié)系數(shù)比逐個(gè)調(diào)用高效。參數(shù)說(shuō)明l在后續(xù)wrcoef和waverec里都必須原樣傳入它是重構(gòu)時(shí)定位系數(shù)段的依據(jù)不能隨意修改。實(shí)踐里我遇到過(guò)有人把l弄丟然后自己猜長(zhǎng)度去拼接系數(shù)的情況重構(gòu)出來(lái)波形完全不對(duì)查了半天才發(fā)現(xiàn)是這個(gè)原因。2.3 wrcoef 單支重構(gòu)把 QRS 波群和軸承頻率從噪聲里拎出來(lái)拿到各層系數(shù)后最直觀的分析方式是單支重構(gòu)——把某一層系數(shù)單獨(dú)映射回時(shí)域其他層置零。這樣可以精確觀察特定頻段的信號(hào)形態(tài)。比如心電信號(hào)里第三層細(xì)節(jié)系數(shù) d3 對(duì)應(yīng)的頻段如果落在 0.5-1Hz很可能藏著 QRS 波群的包絡(luò)特征。單支重構(gòu)用wrcoef它比手動(dòng)拼接再waverec更簡(jiǎn)潔因?yàn)閮?nèi)部已經(jīng)幫你做了系數(shù)段的定位和補(bǔ)零。% 對(duì)三層分解結(jié)果做單支重構(gòu) a3 wrcoef(a, c, l, db4, 3); % 第三層近似分量 d3 wrcoef(d, c, l, db4, 3); % 第三層細(xì)節(jié)分量 d1 wrcoef(d, c, l, db4, 1); % 第一層細(xì)節(jié)分量 subplot(311); plot(a3); title(三階近似信號(hào) A3); subplot(312); plot(d3); title(三階細(xì)節(jié)信號(hào) D3); subplot(313); plot(d1); title(一階細(xì)節(jié)信號(hào) D1);這段代碼的邏輯是同時(shí)觀察低頻趨勢(shì)和高頻細(xì)節(jié)。wrcoef的第一個(gè)參數(shù)傳a表示重構(gòu)近似分量傳d表示重構(gòu)細(xì)節(jié)分量最后一個(gè)數(shù)字指定層數(shù)。注意重構(gòu)出來(lái)的信號(hào)長(zhǎng)度和原始信號(hào)一致wrcoef內(nèi)部做了系數(shù)延拓和濾波。實(shí)際使用中我一般先把各層重構(gòu)信號(hào)畫出來(lái)和原始波形對(duì)齊看確認(rèn)哪一層的形態(tài)和故障特征吻合再?zèng)Q定后續(xù)是去噪還是做特征提取。這比直接看系數(shù)數(shù)組直觀得多也更容易跟非技術(shù)背景的同事解釋結(jié)果。3. 閾值去噪從自動(dòng)閾值到手動(dòng)調(diào)參把信噪比提上去3.1 去噪鏈路ddencmp/wdencmp 與 wthrmngr/wthresh 兩條路線小波去噪的核心邏輯是噪聲主要集中在高頻細(xì)節(jié)系數(shù)里把絕對(duì)值小于閾值的系數(shù)壓掉或收縮再重構(gòu)就能在保留信號(hào)突變特征的同時(shí)濾掉噪聲。MATLAB 里有兩套路線一套是ddencmp配合wdencmp的封裝式調(diào)用適合快速驗(yàn)證另一套是wthrmngr配合wthresh的手動(dòng)控制適合精細(xì)調(diào)參。我實(shí)際項(xiàng)目里先用第一條摸清閾值量級(jí)再切到第二條做細(xì)調(diào)。% 路線一自動(dòng)計(jì)算閾值并去噪 [thr, sorh, keepapp] ddencmp(den, wv, signal); clean wdencmp(gbl, signal, db4, 3, thr, sorh, keepapp); % 路線二手動(dòng)管理每層閾值 cd1 detcoef(c, l, 1); thr1 wthrmngr(dw2ddenoLVL, penalhi, cd1); cd1_denoised wthresh(cd1, s, thr1);路線一適合信號(hào)質(zhì)量尚可、噪聲類型不復(fù)雜的情況ddencmp會(huì)返回三個(gè)關(guān)鍵參數(shù)thr是全局閾值sorh是軟硬閾值標(biāo)記s或hkeepapp控制是否保留近似系數(shù)。keepapp 設(shè)置為 1 意味著只處理細(xì)節(jié)系數(shù)這樣可以避免低頻有用分量被誤傷。路線二的優(yōu)勢(shì)在于可以逐層設(shè)置不同閾值——比如第一層細(xì)節(jié)系數(shù)噪聲占比高閾值調(diào)大第三層細(xì)節(jié)可能藏著有用特征閾值調(diào)小甚至不做處理。參數(shù)說(shuō)明wthrmngr的penalhi是一種基于懲罰函數(shù)的閾值規(guī)則適合脈沖型噪聲如果你處理的是高斯白噪聲換成sqtwolog或rigrsure效果可能更好。3.2 軟閾值 vs 硬閾值脈沖噪聲選哪個(gè)光滑信號(hào)選哪個(gè)閾值處理方式直接決定重構(gòu)信號(hào)的形態(tài)。硬閾值是把低于閾值的系數(shù)直接置零超過(guò)閾值的保持原值波形會(huì)保留尖銳邊緣但可能在跳變點(diǎn)附近產(chǎn)生振鈴。軟閾值是把所有系數(shù)都往零方向收縮一個(gè)閾值量波形更光滑但會(huì)略微壓縮信號(hào)的幅度。具體選哪個(gè)取決于下游任務(wù)——如果是要檢測(cè) QRS 波群的精確位置硬閾值的定位更準(zhǔn)如果要提取信號(hào)包絡(luò)或做疲勞分析軟閾值的平滑性更合適。% 比較軟硬閾值對(duì)同一段含噪信號(hào)的影響 thr 0.2; % 先固定閾值簡(jiǎn)化對(duì)比 d_hard wthresh(cd1, h, thr); d_soft wthresh(cd1, s, thr); % 重構(gòu)后對(duì)比噪聲殘余量 rec_hard waverec([approx; d_hard; cd2; cd3], l, db4); rec_soft waverec([approx; d_soft; cd2; cd3], l, db4); noise_hard signal - rec_hard; noise_soft signal - rec_soft; fprintf(硬閾值殘余噪聲能量%.4f\n, sum(noise_hard.^2)); fprintf(軟閾值殘余噪聲能量%.4f\n, sum(noise_soft.^2));這段代碼的邏輯是固定其他條件只換閾值方式用殘余噪聲能量做量化對(duì)比。wthresh的第二個(gè)參數(shù)h表示硬閾值s表示軟閾值第三個(gè)是閾值大小。waverec接收的是處理后的系數(shù)拼接數(shù)組注意[approx; d_hard; cd2; cd3]的順序必須和原來(lái)wavedec輸出的順序一致否則重構(gòu)就是錯(cuò)的。我實(shí)際測(cè)試時(shí)軟閾值在脈沖噪聲場(chǎng)景下殘余噪聲能量通常比硬閾值低 10% 左右但如果是處理含瞬態(tài)沖擊的信號(hào)硬閾值的峰值保持能力更強(qiáng)需要自己在任務(wù)指標(biāo)上做平衡。3.3 去噪效果驗(yàn)證snr 計(jì)算與重構(gòu)波形對(duì)比去噪不是看波形順眼就完了得量化。常見做法是計(jì)算信噪比提升但前提是你有干凈的參考信號(hào)。實(shí)驗(yàn)室里可以用已知信號(hào)疊加合成噪聲來(lái)驗(yàn)證工程現(xiàn)場(chǎng)的實(shí)測(cè)數(shù)據(jù)沒有真值就得靠多個(gè)頻段的重構(gòu)對(duì)比來(lái)判斷。% 計(jì)算去噪前后的信噪比 snr_before snr(signal); snr_after snr(rec_soft); fprintf(去噪前 SNR%.2f dB去噪后 SNR%.2f dB\n, snr_before, snr_after); % 如果手頭有干凈的參考信號(hào)用均方誤差曲線更準(zhǔn)確 % 假設(shè) clean_signal 是干凈的原始波形 mse_before mean((signal - clean_signal).^2); mse_after mean((rec_soft - clean_signal).^2); fprintf(MSE 從 %.4f 降到 %.4f\n, mse_before, mse_after);snr函數(shù)在 MATLAB 里可以直接對(duì)信號(hào)計(jì)算以 dB 為單位的信噪比默認(rèn)把信號(hào)均值當(dāng)直流分量處理。注意事項(xiàng)如果信號(hào)本身包含強(qiáng)直流成分snr會(huì)把直流也算進(jìn)信號(hào)功率里導(dǎo)致信噪比虛高。我在處理振動(dòng)信號(hào)時(shí)習(xí)慣先做去趨勢(shì)或高通濾波再算。信噪比提升 5dB 左右對(duì)心電信號(hào)來(lái)說(shuō)是能明顯感知的變化但對(duì)故障診斷來(lái)說(shuō)更重要的是特征頻段的能量是否被保留——有時(shí)候全局信噪比提升了關(guān)鍵頻率成分卻衰減了得不償失。所以驗(yàn)證環(huán)節(jié)一定得把去噪前后的頻譜圖疊在一起對(duì)比確認(rèn)無(wú)損。4. 特征提取與可視化用小波系數(shù)定位故障頻段4.1 分解系數(shù)鋪開看subplot 四聯(lián)圖里的頻段信息做完分解和去噪接下來(lái)是特征提取。這一步的核心是判斷信號(hào)里哪些頻段攜帶了有用信息哪些是純?cè)肼?。我的?xí)慣是把各層重構(gòu)信號(hào)按從低頻到高頻的順序鋪開畫在同一個(gè) figure 里先做一次人工頻段掃描。高頻細(xì)節(jié) D1 如果出現(xiàn)穩(wěn)定的周期性波動(dòng)這往往不是噪聲而是某些高頻部件的故障特征低頻近似 A3 則反映了信號(hào)的基線趨勢(shì)或緩慢漂移。figure; subplot(4, 1, 1); plot(t, recA3); title(近似信號(hào) A3低頻趨勢(shì)); subplot(4, 1, 2); plot(t, recD3); title(細(xì)節(jié)信號(hào) D3); subplot(4, 1, 3); plot(t, recD2); title(細(xì)節(jié)信號(hào) D2); subplot(4, 1, 4); plot(t, recD1); title(細(xì)節(jié)信號(hào) D1); xlabel(時(shí)間 (s));這段代碼是標(biāo)準(zhǔn)的時(shí)間對(duì)齊檢查方法。t是和時(shí)間序列對(duì)應(yīng)的橫坐標(biāo)從linspace(0, 5, 5000)生成。注意subplot的排列順序——從上到下頻率遞增這樣高頻信息總是排在最后方便對(duì)比各層之間的波形關(guān)系。實(shí)際項(xiàng)目中我遇到過(guò) D2 層出現(xiàn)周期性波動(dòng)的情況對(duì)應(yīng)軸承外圈的固有頻率單純看原始波形根本發(fā)現(xiàn)不了因?yàn)楸桓哳l的噪聲蓋住了。通過(guò)小波分層故障頻率就暴露出來(lái)了。另外注意信號(hào)長(zhǎng)度與層數(shù)的關(guān)系如果你分解后某一段長(zhǎng)度不足檢查是不是層數(shù)超過(guò)了log2(N)的極限。4.2 cwt 連續(xù)小波時(shí)頻圖尺度-時(shí)間熱圖怎么看多層分解看的是概貌連續(xù)小波變換cwt則能把時(shí)間-尺度頻率關(guān)系變成一張二維熱圖看的是細(xì)節(jié)。它特別適合瞬態(tài)特征的定位——比如旋轉(zhuǎn)機(jī)械啟動(dòng)瞬間的沖擊、電機(jī)堵轉(zhuǎn)時(shí)的頻率跳變。和wavedec的多尺度分解不同cwt會(huì)遍歷設(shè)定范圍內(nèi)所有尺度輸出一個(gè)二維系數(shù)矩陣對(duì)應(yīng)橫軸為時(shí)間、縱軸為尺度或頻率。% 連續(xù)小波變換觀察信號(hào)時(shí)頻分布 scales 1:64; coefs cwt(signal, scales, db4); imagesc(t, scales, abs(coefs)); set(gca, YDir, normal); xlabel(時(shí)間 (s)); ylabel(尺度); colorbar;這段代碼的scales參數(shù)決定頻率分辨率的精細(xì)程度1:64在采樣率 1000Hz 的條件下覆蓋了從低頻到高頻的足夠范圍。abs(coefs)取系數(shù)幅值因?yàn)閺?fù)數(shù)小波變換的系數(shù)包含相位信息畫熱圖時(shí)用幅值更直觀。set(gca, YDir, normal)是為了讓坐標(biāo)軸從下到上按尺度遞增排列不然 MATLAB 默認(rèn)的imagesc會(huì)把 Y 軸倒置。熱圖里顏色越亮的區(qū)域表示該時(shí)間點(diǎn)在該尺度下的能量越集中——出現(xiàn)異常時(shí)特定頻帶的亮色區(qū)域會(huì)明顯移動(dòng)或新增。如果信號(hào)里有未知的周期性干擾這張圖能幫你在幾秒內(nèi)鎖定它的位置和頻帶。4.3 高頻細(xì)節(jié)系數(shù)的隱藏周期對(duì) D1 做 FFT 的騷操作有個(gè)實(shí)用技巧是對(duì)高頻細(xì)節(jié)系數(shù)做二次 FFT而不是直接對(duì)原始信號(hào)做 FFT。原理是原始信號(hào)里低頻成分能量很強(qiáng)FFT 結(jié)果中高頻分量的幅度會(huì)被壓制細(xì)節(jié)系數(shù)經(jīng)過(guò)小波分解后已經(jīng)濾掉了低頻主導(dǎo)成分此時(shí)做 FFT 更容易暴露隱藏的高頻周期。這個(gè)操作在機(jī)械故障診斷里尤其好用——設(shè)備在正常運(yùn)行時(shí) D1 層系數(shù)幾乎全是噪聲一旦某部件出現(xiàn)早期磨損D1 層會(huì)突然出現(xiàn)一個(gè)穩(wěn)定的周期頻率成分。% 對(duì)第一層細(xì)節(jié)系數(shù)做 FFT尋找隱藏周期成分 fs 1000; N length(recD1); f (0:N-1) * fs / N; D1_fft abs(fft(recD1)); [max_val, max_idx] max(D1_fft(2:floor(N/2))); % 跳過(guò)直流分量 fprintf(D1 層占優(yōu)頻率%.2f Hz幅度%.3f\n, f(max_idx1), max_val); plot(f(1:N/2), D1_fft(1:N/2)); xlabel(頻率 (Hz)); ylabel(幅度); title(D1 細(xì)節(jié)系數(shù)頻譜);代碼邏輯是提取 D1 重構(gòu)信號(hào)的頻譜峰值找到占優(yōu)頻率點(diǎn)。D1_fft(2:floor(N/2))跳過(guò)首個(gè)直流分量取一半頻譜是因?yàn)?FFT 結(jié)果對(duì)稱。這里的關(guān)鍵是fs必須和真實(shí)采樣率一致我遇到過(guò)采樣率寫錯(cuò)導(dǎo)致頻率軸偏移把故障特征頻率算到完全錯(cuò)誤位置的案例。FFT 點(diǎn)數(shù)不夠時(shí)可以補(bǔ)零到下一個(gè) 2 的冪分辨率不夠的時(shí)候加窗Hamming 或 Blackman抑制頻譜泄漏。這個(gè)方法說(shuō)到底是個(gè)組合拳——先用小波把頻段隔離出來(lái)再用 FFT 把頻段內(nèi)的周期精確量化比單獨(dú)用任何一個(gè)工具都可靠得多。5. 小波分析的常見坑邊界效應(yīng)、層級(jí)上限與基函數(shù)錯(cuò)配5.1 分解層數(shù)超過(guò) log2(N)系數(shù)爆炸與計(jì)算時(shí)間失控現(xiàn)象設(shè)定 7 層分解后wavedec運(yùn)行明顯變慢輸出的l數(shù)組最后一段長(zhǎng)度變得很小重構(gòu)出來(lái)的信號(hào)在兩端出現(xiàn)嚴(yán)重畸變。原因分解層數(shù)受信號(hào)長(zhǎng)度限制。第 k 層系數(shù)長(zhǎng)度約為原始長(zhǎng)度的2^(-k1)倍當(dāng)層數(shù)接近log2(N)時(shí)最后一層系數(shù)只剩幾個(gè)點(diǎn)信息量幾乎為零卻還占著計(jì)算資源。更麻煩的是層數(shù)越多邊界延拓的影響范圍越廣信號(hào)兩端的失真會(huì)沿著重構(gòu)路徑逐級(jí)放大。解決經(jīng)驗(yàn)法則是分解層數(shù)不超過(guò)floor(log2(N)) - 1。比如信號(hào)長(zhǎng)度 5000log2(5000)約 12.3實(shí)際取 5-6 層就夠了。我處理 1000Hz 采樣率的振動(dòng)信號(hào)時(shí)默認(rèn)從 3 層開始根據(jù)目標(biāo)頻段逐步增加。另外注意超過(guò) 5 層的分解中高頻細(xì)節(jié)基本只剩噪聲去噪時(shí)可以先看各層能量占比再?zèng)Q定保留哪些層。% 自檢分解層數(shù)是否合理 N length(signal); max_level floor(log2(N)) - 1; fprintf(信號(hào)長(zhǎng)度 %d安全分解層數(shù)上限%d\n, N, max_level); % 各層能量占比計(jì)算決定保留哪些層 energy zeros(1, max_level); for k 1:max_level d wrcoef(d, c, l, db4, k); energy(k) sum(d.^2); end bar(energy);5.2 基函數(shù)與信號(hào)形態(tài)不匹配EEG 的 alpha 波被切碎現(xiàn)象處理某同學(xué)的腦電數(shù)據(jù)時(shí)用 db4 做 5 層分解alpha 波頻段8-13Hz的重構(gòu)波形呈鋸齒狀原本的節(jié)律性完全丟失。換成 bior3.5 后alpha 波的形態(tài)恢復(fù)正常。原因db4 的尺度函數(shù)和小波函數(shù)都是不對(duì)稱的相位非線性較強(qiáng)。腦電信號(hào)中的 alpha 波是近似正弦的振蕩信號(hào)需要小波基具有良好的線性相位特性否則重構(gòu)時(shí)各頻率分量之間會(huì)出現(xiàn)相位失真。解決先做一次快速小波基對(duì)比實(shí)驗(yàn)把候選基函數(shù)db4、db8、sym5、coif3、bior3.5各分解一次重構(gòu)目標(biāo)頻段后計(jì)算與原始信號(hào)的相關(guān)系數(shù)取相關(guān)系數(shù)最高者。這個(gè)實(shí)驗(yàn) 5 分鐘就能跑完比按經(jīng)驗(yàn)直接拍板靠譜得多。% 小波基對(duì)比腳本用相關(guān)系數(shù)選基函數(shù) wnames {db4, db8, sym5, coif3, bior3.5}; corrs zeros(1, length(wnames)); for i 1:length(wnames) [c, l] wavedec(eeg_signal, 5, wnames{i}); rec wrcoef(a, c, l, wnames{i}, 5); corrs(i) corr(eeg_signal, rec); end [best_val, best_idx] max(corrs); fprintf(最佳基函數(shù)%s相關(guān)系數(shù)%.4f\n, wnames{best_idx}, best_val);5.3 邊界畸變信號(hào)首尾重構(gòu)變形怎么處理現(xiàn)象重構(gòu)后的信號(hào)在開始和結(jié)束的幾百個(gè)采樣點(diǎn)出現(xiàn)大幅擺動(dòng)波形完全偏離原始信號(hào)特征。原因小波分解時(shí)對(duì)信號(hào)兩端做了延拓延拓方式零延拓、對(duì)稱延拓等會(huì)影響邊界附近的系數(shù)計(jì)算。分解層數(shù)越多邊界影響的區(qū)域越寬。解決最直接的辦法是先把邊界區(qū)域砍掉。如果整段信號(hào)的目標(biāo)頻段信息集中在中間區(qū)域直接在重構(gòu)后截取中間部分使用。另外可以在分解前給信號(hào)做鏡像延拓MATLAB 支持的wextend函數(shù)可以擴(kuò)展信號(hào)長(zhǎng)度處理完再截回原長(zhǎng)。次優(yōu)方案是用wavemngr查看當(dāng)前小波的邊界處理方式換成邊界失真小的基函數(shù)。% 處理邊界畸變對(duì)信號(hào)兩端做鏡像延拓后再分解 ext_len 2^5; % 延拓長(zhǎng)度取分解層數(shù)對(duì)應(yīng)的最大影響范圍 signal_ext wextend(1D, sym, signal, ext_len); [c, l] wavedec(signal_ext, 5, db4); rec_ext waverec(c, l, db4); rec rec_ext(ext_len1:end-ext_len); % 截回原始信號(hào)長(zhǎng)度5.4 自動(dòng)閾值不靠譜特殊信號(hào)還是要手動(dòng)調(diào)參現(xiàn)象用ddencmp自動(dòng)生成的閾值處理一段混合信號(hào)發(fā)現(xiàn)去噪后信號(hào)中某些瞬態(tài)沖擊的幅度被嚴(yán)重壓縮。原因ddencmp默認(rèn)的閾值計(jì)算基于全局噪聲水平估計(jì)當(dāng)信號(hào)本身不是均勻噪聲分布時(shí)自動(dòng)閾值會(huì)偏大把有用信號(hào)細(xì)節(jié)也當(dāng)成噪聲處理了。解決遇到這種情況換成逐層閾值策略。先看每層細(xì)節(jié)系數(shù)的直方圖——噪聲層的系數(shù)分布近似高斯、集中在零附近信號(hào)層的系數(shù)直方圖有明顯的重尾。根據(jù)這個(gè)特征手動(dòng)設(shè)定閾值一般噪聲層的閾值取該層標(biāo)準(zhǔn)差的 2-3 倍信號(hào)層的閾值取標(biāo)準(zhǔn)差的 0.5-1 倍。% 逐層手動(dòng)閾值根據(jù)每層系數(shù)的標(biāo)準(zhǔn)差設(shè)定 std_d1 std(cd1); thr_d1 2.5 * std_d1; % 噪聲層取高閾值 std_d3 std(cd3); thr_d3 0.8 * std_d3; % 信號(hào)層取低閾值 cd1_clean wthresh(cd1, s, thr_d1); cd3_clean wthresh(cd3, s, thr_d3);6. 進(jìn)階玩法小波包自適應(yīng)頻段篩選與批處理加速6.1 wpdec 小波包分解能量節(jié)點(diǎn)自動(dòng)定位主成分wavedec只分解低頻分支高頻分支不再細(xì)分這對(duì)某些高頻成分豐富的信號(hào)來(lái)說(shuō)分辨率不夠。改用wpdec做小波包分解高低頻都繼續(xù)拆分然后通過(guò)能量分析自動(dòng)找到主成分所在的頻段節(jié)點(diǎn)特別適合處理未知類型的混合信號(hào)。tree wpdec(signal, 3, db4); % 3 層小波包分解 nodes bestlevt(tree, tplog); % 選擇能量最大的節(jié)點(diǎn)集合 energy sum(wpcoef(tree, nodes).^2, 2); % 計(jì)算各節(jié)點(diǎn)能量 bar(energy); xlabel(節(jié)點(diǎn)編號(hào)); ylabel(能量);bestlevt基于最小熵準(zhǔn)則自動(dòng)選擇最能代表信號(hào)信息的節(jié)點(diǎn)組合。注意小波包分解的計(jì)算量隨層數(shù)指數(shù)增長(zhǎng)3 層是 8 個(gè)節(jié)點(diǎn)5 層就是 32 個(gè)實(shí)時(shí)系統(tǒng)慎用。我當(dāng)時(shí)用這套方法從一段混合傳感器數(shù)據(jù)里自動(dòng)篩出了設(shè)備啟停的瞬態(tài)特征頻段比人工掃圖快了一個(gè)量級(jí)。6.2 wavemngr 注冊(cè)自定義小波造輪子前先檢查消失矩如果自帶小波基都不滿足要求wavemngr可以注冊(cè)自定義小波。wavemngr(add, MyWavelet, mywt, mydwt.m);注冊(cè)后就能在wavedec系列函數(shù)里直接傳MyWavelet。但建議先檢查兩個(gè)核心指標(biāo)消失矩階數(shù)決定能壓制幾階多項(xiàng)式趨勢(shì)和正則性決定重構(gòu)平滑度。消失矩至少要為 2否則低頻逼近會(huì)殘留線性趨勢(shì)導(dǎo)致解橋正則性太差的話重構(gòu)波形會(huì)有毛刺。6.3 parfor 并行加速?gòu)?3 小時(shí)到 20 分鐘批量處理同類數(shù)據(jù)時(shí)小波分解本身是耗時(shí)的但各段數(shù)據(jù)之間完全獨(dú)立天然適合并行。% 初始化 Parallel Computing Toolbox parpool(local, 4); % 批量處理 200 段信號(hào) results cell(1, 200); parfor i 1:200 [c, l] wavedec(data_series{i}, 5, db4); [thr, sorh, keepapp] ddencmp(den, wv, data_series{i}); processed{i} wdencmp(gbl, c, l, db4, 5, thr, sorh, keepapp); end注意parfor循環(huán)體里不能依賴循環(huán)間的共享變量我在第一次用的時(shí)候踩了這個(gè)坑——在循環(huán)內(nèi)引用前一步處理的數(shù)據(jù)MATLAB 直接報(bào)錯(cuò)。改成每個(gè)迭代自己加載數(shù)據(jù)后就正常了。200 組數(shù)據(jù)從串行 3 小時(shí)縮短到 20 分鐘提速接近 9 倍效果相當(dāng)可觀。從那以后我每次批量處理前都會(huì)強(qiáng)制檢查兩個(gè)東西一是數(shù)據(jù)讀取是否在迭代內(nèi)獨(dú)立完成二是循環(huán)體內(nèi)是否引用了外部變量——確認(rèn)沒問(wèn)題后再開并行池這個(gè)習(xí)慣幫我在多個(gè)項(xiàng)目里少走了不少?gòu)澛?。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取