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