戰(zhàn):wavelet.zip工具箱的MATLAB去噪與信號(hào)分析指南)
簡(jiǎn)介這份MATLAB小波變換擴(kuò)展工具包面向信號(hào)處理、圖像分析、故障診斷等領(lǐng)域的開發(fā)者與研究人員用于彌補(bǔ)標(biāo)準(zhǔn)小波工具箱在快速變換、灰度可視化、自定義濾波與信號(hào)重構(gòu)方面的不足。壓縮包共1447個(gè)文件大小14.08MB核心為1108個(gè)m腳本同時(shí)配備183個(gè)mat數(shù)據(jù)文件、68個(gè)png圖像、26個(gè)fig圖形窗口以及html說(shuō)明文檔便于直接加載測(cè)試和可視化驗(yàn)證。包內(nèi)wavefast實(shí)現(xiàn)快速小波變換以降低計(jì)算開銷wave2gray將小波系數(shù)轉(zhuǎn)為灰度圖像wavefilter支持去噪與特征提取wavecopy用于系數(shù)復(fù)制與實(shí)驗(yàn)對(duì)比waveback完成信號(hào)重構(gòu)組合起來(lái)可搭建完整的小波分析流程。此外資源覆蓋故障診斷、金融時(shí)間序列、地震信號(hào)處理、醫(yī)學(xué)圖像分析等典型場(chǎng)景目前已有1717人學(xué)習(xí)下載適合需要擴(kuò)展MATLAB小波能力、開展實(shí)驗(yàn)對(duì)比或工程落地的中高級(jí)用戶。1. 小波變換離線工具箱wavelet.zip 到底能幫你省多少事小波變換在信號(hào)處理里的地位不用多言但真正上手時(shí)MATLAB自帶的小波工具箱函數(shù)分散、參數(shù)繁多尤其做振動(dòng)信號(hào)去噪時(shí)每次都要查一遍wden、wavedec、thselect的用法閾值調(diào)起來(lái)像玄學(xué)。所以當(dāng)我在某技術(shù)社區(qū)里看到這份wavelet.zip小波工具箱時(shí)第一反應(yīng)是“這不就是把常用函數(shù)封裝好了嘛”。解壓后掃了一眼它把多級(jí)分解、重構(gòu)、去噪、小波包都整理成了統(tǒng)一接口還帶一個(gè)demo腳本確實(shí)能省半天時(shí)間。這篇筆記會(huì)講清楚它適合什么場(chǎng)景、怎么安裝、參數(shù)怎么調(diào)、坑在哪里正在做故障診斷、振動(dòng)分析、圖像去噪的從業(yè)者可以直接照著做剛接觸小波變換的初學(xué)者也能按步驟跑通。2. 小波變換的核心概念與wavelet.zip的功能地圖這一章先建立對(duì)小波變換的直覺(jué)再對(duì)照zip包里的函數(shù)清單知道每段代碼在干什么。理解這三個(gè)關(guān)鍵參數(shù)——小波基、分解層數(shù)、閾值規(guī)則后面調(diào)參就有了方向。2.1 從傅里葉到小波時(shí)間-頻率局部化意味著什么傅里葉變換把信號(hào)拆成正弦波的疊加得到全局頻譜但丟掉了時(shí)間位置信息。兩個(gè)不同時(shí)刻出現(xiàn)的突變信號(hào)可能頻譜幾乎一樣傅里葉看不出差異。短時(shí)傅里葉變換用固定窗去截信號(hào)時(shí)間分辨率有了但頻率分辨率被窗長(zhǎng)鎖死窗短了高頻能分辨低頻就模糊窗長(zhǎng)了低頻清楚了突變點(diǎn)又被抹平。小波變換用“可伸縮的窗”——母小波通過(guò)平移和縮放產(chǎn)生一系列子波與信號(hào)做內(nèi)積在低頻用長(zhǎng)窗獲得細(xì)頻率分辨率在高頻用短窗捕捉突變時(shí)刻這就是時(shí)間-頻率局部化。離散小波變換DWT通過(guò)濾波器組實(shí)現(xiàn)信號(hào)經(jīng)過(guò)低通和高通濾波器分別得到近似分量和細(xì)節(jié)分量再對(duì)近似分量繼續(xù)分解形成多級(jí)結(jié)構(gòu)。每一步輸出兩個(gè)部分近似系數(shù)cA和細(xì)節(jié)系數(shù)cD。如果用MATLAB底層函數(shù)查看濾波器形狀可以這樣[LoD, HiD, LoR, HiR] wfilters(db4); % 低通分解濾波器的沖擊響應(yīng) subplot(2,1,1); stem(LoD); title(db4 低通分解濾波器 LoD); subplot(2,1,2); stem(HiD); title(db4 高通分解濾波器 HiD);LoD是低通分解濾波器系數(shù)HiD是高通分解濾波器LoR、HiR是重構(gòu)對(duì)應(yīng)的低通和高通濾波器??催@些系數(shù)的形狀就能直觀感受到短的小波比如db2突變更集中長(zhǎng)的小波比如db8過(guò)渡更平滑。這個(gè)特性直接決定了小波基的選型。實(shí)際做DWT時(shí)不需要手動(dòng)算卷積MATLAB的wavedec已經(jīng)封裝好了。wavelet.zip里的my_multilevel_dwt就是在wavedec外面加了一層輸入校驗(yàn)和邊界處理調(diào)用起來(lái)更省心。2.2 wavelet.zip 的文件結(jié)構(gòu)每個(gè)函數(shù)負(fù)責(zé)什么我解壓后把wavelet.zip里的文件列了個(gè)清單典型工具包會(huì)包含以下文件不同版本可能有增減但核心功能基本一致??梢杂胐ir命令查看d dir(wavelet_tool/*.m); for i 1:length(d) disp(d(i).name); end我見(jiàn)到的版本通常包含這些文件作用核心接口wavelet_filters.m獲取指定小波基的濾波器系數(shù)[LoD,HiD,LoR,HiR] wavelet_filters(db4)my_multilevel_dwt.m多級(jí)離散小波分解[C,L] my_multilevel_dwt(x,db4,5)my_waverec.m多級(jí)重構(gòu)從系數(shù)還原信號(hào)[xr,err] my_waverec(C,L,db4,5)wavelet_denoise.m小波去噪主函數(shù)[xd,C,L] wavelet_denoise(x,db4,5,heursure)my_wavelet_packet.m小波包分解[wp,wpt] my_wavelet_packet(x,db4,3)example_demo.m示例腳本生成含噪信號(hào)并演示去噪直接運(yùn)行注意wavelet_denoise返回的C和L是去噪時(shí)用到的原始分解系數(shù)不是去噪后的系數(shù)。如果想基于去噪后的信號(hào)做特征提取需要先得到去噪信號(hào)再重新分解。2.3 小波基選型dbN、symN、coifN 到底有什么區(qū)別小波基直接影響濾波器長(zhǎng)度、正交性、對(duì)稱性和消失矩。dbN系列Daubechies小波是經(jīng)典正交小波N是消失矩階數(shù)N越大濾波器越長(zhǎng)頻帶劃分越精細(xì)但計(jì)算量和邊界失真也增加。symN系列是dbN的改進(jìn)對(duì)稱性更好重構(gòu)信號(hào)相位失真更小適合需要保真的信號(hào)處理。coifN系列在對(duì)稱性和消失矩之間折中常用于圖像處理。選型有個(gè)通用經(jīng)驗(yàn)先跑一個(gè)demo用相同分解層數(shù)對(duì)比不同小波基的重構(gòu)誤差和去噪效果。大部分情況下db4和sym4是可用且穩(wěn)定的起點(diǎn)。當(dāng)信號(hào)以瞬態(tài)突變?yōu)橹鞫桃恍┑男〔ㄈ鏳b2更敏感當(dāng)信號(hào)是平滑準(zhǔn)周期信號(hào)長(zhǎng)一些的小波如db8或sym6能更好分離頻帶。提示不要被“N越大越好”誤導(dǎo)過(guò)長(zhǎng)的濾波器會(huì)把噪聲也當(dāng)成有效分量計(jì)算量還翻倍。驗(yàn)證小波基可用性的一個(gè)直接辦法是用wavedec和waverec做分解-重構(gòu)看歸一化誤差load(noisy_signal.mat, x); for N 1:4 wname sprintf(db%d, N); [C, L] wavedec(x, 5, wname); xr waverec(C, L, wname); err(N) norm(x - xr) / norm(x); end bar(1:4, err); xlabel(db N); ylabel(歸一化重構(gòu)誤差);只要小波基名正確、系數(shù)沒(méi)動(dòng)err基本在1e-12量級(jí)。如果某個(gè)N對(duì)應(yīng)的誤差明顯變大說(shuō)明該小波基和信號(hào)長(zhǎng)度搭配有問(wèn)題或者你的調(diào)用代碼有誤。這就是用wavelet.zip做封裝時(shí)應(yīng)該保留的驗(yàn)證手段。3. 把wavelet.zip裝進(jìn)MATLAB路徑設(shè)置與第一個(gè)demo跑通這一章的目標(biāo)是讓讀者在五分鐘內(nèi)跑通整個(gè)工具包。我們分四步解壓、驗(yàn)證路徑、運(yùn)行demo、替換自己的數(shù)據(jù)。3.1 解壓與路徑設(shè)置別讓MATLAB找不到函數(shù)打開MATLAB把當(dāng)前工作目錄切到wavelet.zip所在文件夾然后執(zhí)行解壓和添加路徑。常見(jiàn)做法是直接在命令窗口里用unzipunzip(wavelet.zip, wavelet_tool); cd(wavelet_tool); addpath(genpath(pwd)); savepath;第一行把壓縮包解壓到當(dāng)前目錄下的wavelet_tool文件夾第二行進(jìn)入該文件夾第三行用genpath遞歸地將該文件夾及所有子目錄添加到MATLAB搜索路徑這樣private子目錄或示例數(shù)據(jù)目錄里的函數(shù)也能被找到第四行把路徑保存到pathdef.m下次啟動(dòng)MATLAB仍在。之后驗(yàn)證函數(shù)是否可見(jiàn)which wavelet_denoise which my_multilevel_dwt如果輸出路徑說(shuō)明設(shè)置成功。如果提示找不到先檢查文件名是否帶了隱藏?cái)U(kuò)展名比如wavelet_denoise.m.txt在Windows資源管理器里改后綴名時(shí)容易犯這個(gè)毛病。3.2 運(yùn)行example_demo.m先看標(biāo)準(zhǔn)輸出工具包自帶的example_demo.m會(huì)生成一段含噪正弦信號(hào)疊加隨機(jī)噪聲和兩個(gè)間斷突變?cè)僬{(diào)用去噪函數(shù)繪制原始、去噪和殘差三條曲線。直接運(yùn)行example_demo正常會(huì)彈出圖窗左邊是原始含噪信號(hào)右邊是去噪信號(hào)。如果沒(méi)報(bào)錯(cuò)說(shuō)明核心鏈路是通的。我一般會(huì)先看殘差幅值如果殘差里仍有明顯的規(guī)則波形說(shuō)明閾值太寬松如果殘差全是白噪聲說(shuō)明去噪強(qiáng)度合適。有些版本的example_demo.m內(nèi)部可能是這樣生成信號(hào)的Fs 1000; t 0:1/Fs:1; x sin(2*pi*50*t) 0.5*sin(2*pi*150*t); x([300, 700]) x([300, 700]) 1; % 模擬突變 x x 0.3*randn(size(t));生成之后再調(diào)用wavelet_denoise。你可以把這段代碼復(fù)制進(jìn)自己的腳本改改頻率和噪聲強(qiáng)度快速理解參數(shù)影響。3.3 用自己的信號(hào)調(diào)用wavelet_denoise假設(shè)你有一段采樣率為1kHz的振動(dòng)信號(hào)x存儲(chǔ)為列向量直接調(diào)用load(my_vibration.mat, x); Fs 1000; t (0:length(x)-1)/Fs; [xd, C, L] wavelet_denoise(x, sym4, 5, heursure); plot(t, x); hold on; plot(t, xd, r); legend(原始, 去噪);wavelet_denoise的第一個(gè)輸出是去噪信號(hào)第二、第三分別是分解系數(shù)和長(zhǎng)度記錄供后續(xù)重構(gòu)或分析使用。這里小波基選了sym4層數(shù)是5閾值規(guī)則是heursure。這四個(gè)參數(shù)是核心下一章詳細(xì)展開。如果你的信號(hào)長(zhǎng)度不是2的冪完全沒(méi)問(wèn)題wavedec本身不要求長(zhǎng)度等于2的冪但分解層數(shù)受長(zhǎng)度限制太多時(shí)會(huì)報(bào)錯(cuò)下一章講。3.4 運(yùn)行過(guò)程的常見(jiàn)警告與初步排查首次運(yùn)行時(shí)如果MATLAB版本較舊heursure可能不被支持因?yàn)樵撘?guī)則在新版本中才引入。解決方法是把參數(shù)改為heursure對(duì)應(yīng)的舊式寫法或者先查看doc thselect確認(rèn)。如果demo運(yùn)行后沒(méi)有圖窗檢查是否有figure命令被注釋了或當(dāng)前是無(wú)圖形桌面的系統(tǒng)環(huán)境。桌面版MATLAB一般不會(huì)遇到這個(gè)問(wèn)題但Linux遠(yuǎn)端用-nodisplay模式時(shí)需要注意。還有一個(gè)容易被忽略的地方如果wavelet_denoise函數(shù)內(nèi)部使用的wden依賴Wavelet Toolbox而你的MATLAB沒(méi)有安裝該工具箱會(huì)報(bào)Undefined function wden。wavelet.zip只是封裝不能替代工具箱本身所以安裝前確認(rèn)一下license(test, Wavelet_Toolbox)返回1。4. 核心函數(shù)逐個(gè)拆解分解、重構(gòu)、去噪與小波包的參數(shù)這一章把每個(gè)接口的輸入輸出、內(nèi)部邏輯、參數(shù)含義講透。只有理解了這些才能針對(duì)自己的信號(hào)調(diào)參。4.1 my_multilevel_dwt.m多級(jí)分解的接線口這個(gè)函數(shù)負(fù)責(zé)把一維信號(hào)分解成多層系數(shù)。調(diào)用方式[C, L] my_multilevel_dwt(x, db4, 5);輸入x是待分解信號(hào)wname是小波基名level是分解層數(shù)。輸出C是系數(shù)向量L是各層長(zhǎng)度記錄。C的前L(end)個(gè)元素是最低層近似系數(shù)之后依次是各層細(xì)節(jié)系數(shù)。想提某層細(xì)節(jié)可以用wrcoefd5 wrcoef(d, C, L, db4, 5); plot(d5); title(第5層細(xì)節(jié)系數(shù));my_multilevel_dwt內(nèi)部等效于wavedec但加上了輸入校驗(yàn)行向量會(huì)轉(zhuǎn)成列向量level超過(guò)可能的最大層數(shù)時(shí)自動(dòng)截?cái)嗖⑻崾尽_@樣避免一些誤操作。參數(shù)說(shuō)明level不是越大越好。分解層數(shù)受信號(hào)長(zhǎng)度限制最小可分解長(zhǎng)度滿足N 2^level更長(zhǎng)的小波基會(huì)拉高這個(gè)限制。層數(shù)越多頻率分辨率越高但每一層系數(shù)樣本數(shù)越少統(tǒng)計(jì)意義下降。對(duì)于1000個(gè)點(diǎn)的信號(hào)用db4試到6層基本就是極限。4.2 my_waverec.m重構(gòu)系數(shù)并驗(yàn)證完整性重構(gòu)是分解的逆操作[xr, err] my_waverec(C, L, db4, 5);輸出xr是重構(gòu)信號(hào)err是重構(gòu)誤差。理想情況下分解-重構(gòu)后err應(yīng)該在浮點(diǎn)誤差量級(jí)。驗(yàn)證代碼[C, L] my_multilevel_dwt(x, db4, 5); [xr, err] my_waverec(C, L, db4, 5); fprintf(重構(gòu)誤差 %e\n, err);只要小波基一致、系數(shù)沒(méi)被修改err通常在1e-12左右。如果誤差達(dá)到1%以上先檢查C或L是否被改動(dòng)過(guò)再檢查小波基名是否一致。這個(gè)驗(yàn)證方法也是判斷wavelet.zip封裝可靠性的好工具。4.3 wavelet_denoise.m閾值規(guī)則與去噪強(qiáng)度去噪函數(shù)是核心。調(diào)用[xd, C, L] wavelet_denoise(x, db4, 5, sqtwolog);內(nèi)部采用軟閾值或硬閾值閾值規(guī)則支持四種rigrsure無(wú)偏風(fēng)險(xiǎn)估計(jì)、heursure啟發(fā)式、sqtwolog固定閾值、minimaxi極小極大。四種規(guī)則適用場(chǎng)景不同我通常用快速實(shí)驗(yàn)比較對(duì)同一段噪聲信號(hào)分別用四種規(guī)則去噪計(jì)算去噪后殘差的標(biāo)準(zhǔn)差標(biāo)準(zhǔn)差接近噪聲水平的最合適。對(duì)比實(shí)驗(yàn)可以這樣寫rules {rigrsure, heursure, sqtwolog, minimaxi}; for i 1:4 [xd, ~, ~] wavelet_denoise(x, sym4, 5, rules{i}); subplot(2,2,i); plot(t, xd); title(rules{i}); end注意wavelet_denoise返回的C和L是去噪前的系數(shù)不是去噪后的系數(shù)。如果你想用去噪后的系數(shù)做后續(xù)特征提取需要先用my_multilevel_dwt對(duì)去噪信號(hào)重新分解。去噪強(qiáng)度主要受閾值規(guī)則影響。sqtwolog最激進(jìn)適合噪聲強(qiáng)度大且平穩(wěn)的情況rigrsure最保守適合信號(hào)本身有較多高頻成分的情況heursure介于兩者之間是默認(rèn)推薦。如果規(guī)則不合適去噪后要么噪聲殘留多要么有效信號(hào)被削平這時(shí)需要改用分層閾值。4.4 my_wavelet_packet.m小波包分解與頻帶能量多級(jí)DWT只對(duì)低頻繼續(xù)分解高頻細(xì)節(jié)的頻帶寬度固定小波包則對(duì)高頻也做分解得到更細(xì)致的頻帶劃分。調(diào)用[wp, wpt] my_wavelet_packet(x, db4, 3);輸出wp是小波包樹對(duì)象wpt是各節(jié)點(diǎn)系數(shù)。做頻帶特征提取時(shí)通常需要計(jì)算每個(gè)葉子節(jié)點(diǎn)的能量n 2^3; E zeros(n, 1); for i 1:n % 示意按索引從 wpt 中取出第i個(gè)葉子節(jié)點(diǎn)系數(shù) seg wpt{i}; E(i) sum(seg.^2); end E_norm E / sum(E); bar(1:n, E_norm); xlabel(葉子節(jié)點(diǎn)編號(hào)); ylabel(歸一化能量);這里2^3是分解到第3層時(shí)的葉子數(shù)。如果采樣率是2000Hz第1個(gè)葉子大致對(duì)應(yīng)0~125Hz第2個(gè)125~250Hz以此類推。想要聚焦某個(gè)故障特征頻率就增加分解層數(shù)但層數(shù)增加后每個(gè)葉子系數(shù)長(zhǎng)度減半能量估計(jì)的方差變大。所以選擇層數(shù)要在頻率分辨率與統(tǒng)計(jì)穩(wěn)定性之間權(quán)衡。5. 避坑wavelet.zip 使用中的5個(gè)典型問(wèn)題與排查這里從實(shí)際問(wèn)題角度回顧。這些坑不只在wavelet.zip里存在用小波工具箱都會(huì)遇到。我按“現(xiàn)象→原因→解決”來(lái)寫。5.1 函數(shù)變成“Undefined function or variable”怎么辦現(xiàn)象命令窗口輸入wavelet_denoise報(bào)錯(cuò)Undefined function or variable wavelet_denoise。原因絕大多數(shù)時(shí)候是路徑?jīng)]添加成功或者解壓后把文件夾移動(dòng)過(guò)也有可能是文件名誤帶了隱藏?cái)U(kuò)展名比如wavelet_denoise.m.txtWindows下容易犯這個(gè)毛病。解決重新執(zhí)行addpath(genpath(wavelet_tool))然后用which wavelet_denoise確認(rèn)。如果還是找不到檢查文件擴(kuò)展名并在MATLAB里用fileparts打印文件名確認(rèn)實(shí)際狀態(tài)。如果確認(rèn)文件在但which仍是空可以試試rehash toolboxcache刷新緩存。5.2 信號(hào)兩端出現(xiàn)明顯畸變邊界延拓踩坑現(xiàn)象去噪后的信號(hào)兩端有大幅度振蕩中間正常。原因DWT的卷積過(guò)程在信號(hào)邊界做延拓MATLAB默認(rèn)模式是zpd零填充對(duì)非平穩(wěn)信號(hào)會(huì)在端點(diǎn)產(chǎn)生突然跳變導(dǎo)致重構(gòu)后端點(diǎn)失真。分解層數(shù)較多時(shí)尤其明顯。解決改用dwtmode(ppd)或dwtmode(sym)分別對(duì)應(yīng)周期延拓和平滑延拓。我一般會(huì)先用小波分解少量數(shù)據(jù)觀察邊界再?zèng)Q定延拓模式。注意dwtmode是全局設(shè)置調(diào)完記得改回來(lái)否則別的腳本會(huì)受影響。對(duì)邊界的驗(yàn)證方法很簡(jiǎn)單dwtmode(sym); % 改為對(duì)稱延拓 [xd, ~, ~] wavelet_denoise(x, db4, 5, heursure); plot(t, xd); % 觀察兩端是否改善5.3 閾值規(guī)則不合適去噪后信號(hào)被“壓扁”或噪聲殘留現(xiàn)象用sqtwolog去噪后信號(hào)過(guò)于平滑峰值被明顯削弱用rigrsure去噪后仍有明顯噪聲殘留。原因固定閾值在噪聲方差估計(jì)過(guò)高時(shí)把有效信號(hào)當(dāng)噪聲濾掉無(wú)偏風(fēng)險(xiǎn)估計(jì)在信號(hào)本身有較多高頻成分時(shí)會(huì)保留更多細(xì)節(jié)代價(jià)是噪聲也保留。解決先對(duì)噪聲段做方差估計(jì)再用分層閾值替代全局閾值。可以這樣計(jì)算每層噪聲強(qiáng)度[C, L] wavedec(x, 5, db4); noise wnoisest(C, L, 1); % 估計(jì)第一層細(xì)節(jié)的噪聲方差 thr thselect(x, heursure); [xd, C, L] wdencmp(gbl, x, db4, 5, thr, s, 1);如果覺(jué)得效果還不夠就自己算分層閾值把每層細(xì)節(jié)系數(shù)分別過(guò)thselect再做軟閾值處理。5.4 分解層數(shù)超過(guò)極限近似分量全部變成零現(xiàn)象調(diào)用my_multilevel_dwt時(shí)提示分解層數(shù)太大或返回的近似分量全為零。原因信號(hào)長(zhǎng)度不夠。對(duì)于長(zhǎng)度1000的信號(hào)最大層數(shù)可能只有6~8如果用db4這種稍長(zhǎng)的濾波器限制更嚴(yán)格。超過(guò)極限后濾波器長(zhǎng)度大于數(shù)據(jù)段系數(shù)失去意義。解決用wmaxlev查看允許的最大層數(shù)maxLev wmaxlev(length(x), db4); fprintf(最大允許層數(shù): %d\n, maxLev);再取小于等于該值的層數(shù)。如果必須增加層數(shù)就先對(duì)信號(hào)做分段或補(bǔ)零把每段長(zhǎng)度補(bǔ)到2的冪次再設(shè)置層數(shù)。注意補(bǔ)零會(huì)影響邊界補(bǔ)完要結(jié)合5.2的邊界處理。5.5 小波基與信號(hào)不匹配導(dǎo)致重構(gòu)誤差大現(xiàn)象分解-重構(gòu)后誤差遠(yuǎn)大于1e-12信號(hào)波形整體變形。原因小波基不具備正交性或者重構(gòu)濾波器與分解濾波器不匹配。有些開源代碼把LoD和HiD順序?qū)懛椿蛘叻纸庥胐b4重構(gòu)時(shí)寫成dB4大小寫錯(cuò)誤導(dǎo)致濾波器完全不同。解決固定小波基名重構(gòu)時(shí)通過(guò)my_waverec的err輸出檢查。如果誤差大先核對(duì)小波基名再檢查系數(shù)向量是否被修改。血淚經(jīng)驗(yàn)我曾經(jīng)做過(guò)一次實(shí)驗(yàn)分解時(shí)用db4重構(gòu)時(shí)寫成dB4大小寫不一致導(dǎo)致重構(gòu)誤差直接到0.5后來(lái)用strcmp核對(duì)字符串才找到原因?,F(xiàn)在每次跑完分解重構(gòu)我都會(huì)輸出err看一眼不只看波形。6. 進(jìn)階技巧小波包頻帶能量分析與批量處理最后分享一個(gè)具體技巧用小波包計(jì)算頻帶能量占比用于某設(shè)備狀態(tài)識(shí)別。這個(gè)技巧在振動(dòng)信號(hào)診斷里很實(shí)用。假設(shè)有一段8000點(diǎn)的采樣信號(hào)采樣率2000Hz想觀察0~1000Hz內(nèi)頻帶能量比例分解層數(shù)選3層[~, wpt] my_wavelet_packet(x, db4, 3); n 2^3; E zeros(n, 1); for i 1:n seg wpt{i}; % 第i個(gè)葉子節(jié)點(diǎn)系數(shù) E(i) sum(seg.^2); end E_norm E / sum(E); bar(1:n, E_norm); xlabel(葉子節(jié)點(diǎn)); ylabel(歸一化能量);這里wpt{i}的索引方式取決于wavelet.zip的節(jié)點(diǎn)排列約定我一般會(huì)先打開源碼確認(rèn)順序。采樣率2000Hz時(shí)第一層葉子大致對(duì)應(yīng)0~125Hz第二層125~250Hz以此類推。如果想聚焦某個(gè)故障特征頻率就加大分解層數(shù)但每個(gè)葉子系數(shù)長(zhǎng)度減半能量統(tǒng)計(jì)穩(wěn)定性下降需要在分辨率與穩(wěn)定性之間權(quán)衡。批量處理時(shí)我會(huì)寫一個(gè)循環(huán)對(duì)多段信號(hào)提取特征矩陣numSeg length(segList); featMat zeros(numSeg, n); for k 1:numSeg [~, wpt] my_wavelet_packet(segList{k}, sym4, 3); for i 1:n featMat(k, i) sum(wpt{i}.^2); end featMat(k, :) featMat(k, :) / sum(featMat(k, :)); end最后把這個(gè)特征矩陣輸入分類器就能做簡(jiǎn)單的狀態(tài)識(shí)別。需要注意不同段的長(zhǎng)度要一致否則小波包系數(shù)長(zhǎng)度不同能量比較沒(méi)有意義。我的習(xí)慣是每次拿到新信號(hào)先不做去噪直接做三層小波包分解畫出能量柱狀圖再?zèng)Q定該用哪個(gè)頻帶。這個(gè)習(xí)慣幫我避開了很多“先入為主”的調(diào)參錯(cuò)誤。從那以后我每次跑wavelet.zip里的函數(shù)都會(huì)強(qiáng)制先走一遍邊界檢查和誤差校驗(yàn)。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取