戰(zhàn)：從經(jīng)驗(yàn)?zāi)B(tài)分解到HHT時(shí)頻分析)
簡介Matlab EMD經(jīng)驗(yàn)?zāi)B(tài)分解工具箱面向信號處理、故障診斷與振動(dòng)分析領(lǐng)域的研究人員與工程師用于對非線性、非平穩(wěn)信號做自適應(yīng)分解并提取有物理意義的本征模態(tài)函數(shù)分量。壓縮包共38個(gè)文件其中10個(gè)m腳本承擔(dān)核心算法與在線/局部分解入口17個(gè)c文件實(shí)現(xiàn)底層快速計(jì)算11個(gè)h頭文件提供接口聲明整體僅55KB輕量易部署。資源內(nèi)置標(biāo)準(zhǔn)EMD、復(fù)數(shù)擴(kuò)展EMD與C加速實(shí)現(xiàn)可靈活處理實(shí)信號與復(fù)信號適配科研實(shí)驗(yàn)、課程教學(xué)和工程診斷等場景。目前已有729人學(xué)習(xí)下載適合需要快速搭建EMD分析環(huán)境、二次開發(fā)或集成到現(xiàn)有Matlab流程中的用戶是一份精簡實(shí)用的信號分解工具集。1. EMD工具箱非平穩(wěn)信號分析的入場券matlab EMD工具箱解決的是經(jīng)驗(yàn)?zāi)B(tài)分解的落地問題一套能在 MATLAB 里直接跑通的函數(shù)集把非平穩(wěn)、非線性信號自適應(yīng)地拆成若干個(gè)本征模函數(shù)IMF和一條殘差。我之前幫某實(shí)驗(yàn)室處理設(shè)備啟動(dòng)瞬間的振動(dòng)數(shù)據(jù)前 0.5 秒信號直接做 FFT 根本看不出頻率成分換成 EMD 分解之后特征頻段一下子拆開了。這套資源適合三類人剛接觸經(jīng)驗(yàn)?zāi)B(tài)分解、想跳過論文偽代碼的在校學(xué)生需要做 HHT 時(shí)頻分析的信號處理工程師以及手里攢了一堆實(shí)測波形、想低成本提取特征的分析人員。下面從原理開始逐步講到參數(shù)調(diào)整、批量處理和驗(yàn)證技巧。2. 從原理到函數(shù)EMD拆出了什么怎么讀懂IMF2.1 IMF的定義與篩分循環(huán)三次樣條包絡(luò)在做什么EMD 的核心假設(shè)是任意復(fù)雜信號都能分解為若干本征模函數(shù)。IMF 要滿足兩個(gè)條件一是極值點(diǎn)數(shù)量與過零點(diǎn)數(shù)量相等或相差一個(gè)二是上下包絡(luò)的均值在任意點(diǎn)都接近零。第一條保證波形接近正弦的“局域?qū)ΨQ”第二條保證瞬時(shí)頻率有意義。兩個(gè)條件都是經(jīng)驗(yàn)性的這正是“經(jīng)驗(yàn)?zāi)B(tài)分解”名字的由來。以某跨平臺系統(tǒng)的振動(dòng)信號分析為例我的理解是 EMD 本質(zhì)上是一套自適應(yīng)篩分循環(huán)。具體流程是先找信號的局部極大值和極小值用三次樣條分別插值成上包絡(luò)和下包絡(luò)算出包絡(luò)均值原信號減去均值后得到候選分量判斷它是否滿足 IMF 條件不滿足就繼續(xù)重復(fù)篩分滿足就把它抽出去對殘差繼續(xù)執(zhí)行同樣的過程。整個(gè)過程像一層層剝洋蔥每一層都去掉一個(gè)相對高頻的振蕩成分剩下的殘差通常是單調(diào)趨勢或直流偏置。三次樣條是 EMD 的默認(rèn)插值方式因?yàn)樗亩A光滑性能讓包絡(luò)線在連接處不產(chǎn)生突兀拐點(diǎn)。不過樣條也帶來兩個(gè)副作用一是計(jì)算量大二是端點(diǎn)處的外推不受控。前者在長信號上會直接拖慢速度后者會引發(fā)下文要說的端點(diǎn)效應(yīng)。所以在理解原理的時(shí)候把“包絡(luò)生成——均值相減——迭代判定”當(dāng)作一條主線就夠了細(xì)節(jié)參數(shù)先不用背。選 EMD 而不是其它分解方法主要看三點(diǎn)。FFT 只能把整段信號折算成全局頻譜頻率隨時(shí)間的變化會被平均掉小波分解需要提前選基函數(shù)和分解層數(shù)選錯(cuò)了結(jié)果就差很遠(yuǎn)EMD 不需要預(yù)先設(shè)定基函數(shù)完全按信號自身的極值分布來篩分。代價(jià)是計(jì)算結(jié)果對極值點(diǎn)的分布和插值方式敏感同一段數(shù)據(jù)在不同參數(shù)或不同版本下可能篩出不同數(shù)量的 IMF。這是它的經(jīng)驗(yàn)性決定的不是 bug。2.2 工具箱函數(shù)地圖emd 主函數(shù)與 hht 譜分析我在項(xiàng)目里常用的是官方 Signal Processing Toolbox 自帶的 emd 函數(shù)它的標(biāo)準(zhǔn)調(diào)用方式是[imf, residual] emd(x)。imf 輸出是一個(gè)矩陣每一列是一個(gè) IMF 分量列數(shù)由信號本身的復(fù)雜度決定residual 返回最終殘差。官方版本還支持 Name-Value 形式的可選參數(shù)比如控制篩分終止條件、最大 IMF 數(shù)量、插值方式等具體會在第 4 章展開。配套的 hht 函數(shù)用于做 Hilbert-Huang 變換。它的輸入是 imf 矩陣和采樣率輸出是瞬時(shí)頻率、瞬時(shí)幅值和 Hilbert 譜。很多人會把 EMD 和 HHT 混著叫其實(shí)前者是做分解后者是對各 IMF 做 Hilbert 變換并畫時(shí)頻圖。拿到 imf 矩陣之后先看每個(gè)分量波形是否規(guī)律再用 hht 畫譜才能確認(rèn)分解結(jié)果能不能用于特征提取。除此之外工具箱還帶了自帶繪圖邏輯。比如 emd 在 Display 參數(shù)打開時(shí)會在命令行輸出篩分迭代信息配合tiledlayout可以自己拼一個(gè) IMF 逐行排列的波形圖。這個(gè)圖幾乎是驗(yàn)收分解效果的標(biāo)配看一眼 IMF 數(shù)量、每個(gè)分量的幅值量級、殘差形態(tài)基本就能判斷分解是否可信。關(guān)于這套函數(shù)更完整的用法包括各種返回值和繪圖輔助我在做完源碼包的閱讀筆記之后整理了一張函數(shù)與用途對照表放在后面的章節(jié)里。3. 安裝與第一次分解環(huán)境配置到讀圖全流程3.1 版本選擇與路徑配置別讓老腳本坑了你經(jīng)驗(yàn)?zāi)B(tài)分解在 MATLAB 里的可用性經(jīng)歷了幾個(gè)階段。比較舊的第三方腳本是用純 MATLAB 寫的函數(shù)名通常叫emd參數(shù)順序和官方版本不一樣后來官方在 Signal Processing Toolbox 里內(nèi)置了 emd 和 hht函數(shù)命名更規(guī)范底層也做了優(yōu)化。我手邊用過的 2019、2020、2023 版本都自帶這個(gè)函數(shù)基本不用額外下載安裝包。如果你拿到的是第三方工具箱壓縮包常見做法是解壓后把整個(gè)文件夾加入 MATLAB 搜索路徑。路徑配置只需要在 MATLAB 主界面執(zhí)行addpath我一般還會用savepath把路徑持久化不然重啟后又要重新加。判斷當(dāng)前環(huán)境里用的是官方版本還是第三方版本最直接的方法是which emd -all找到結(jié)果列表里第一個(gè) emd.m 的路徑看它在不在 toolbox 目錄下。這樣處理過的環(huán)境問題多了之后我的習(xí)慣是把這一條檢查寫在每個(gè)腳本開頭注釋掉也行。3.2 第一次分解腳本生成仿真信號并畫出 IMF這里直接用可復(fù)現(xiàn)的仿真信號來跑通流程。我選取了一個(gè)帶調(diào)幅特征的目標(biāo)信號疊加一個(gè)低頻正弦和一個(gè)線性趨勢項(xiàng)模擬真實(shí)測試臺里“轉(zhuǎn)子輕微不對中加基座松動(dòng)”的典型組合%% 生成仿真信號調(diào)幅分量 低頻分量 線性趨勢 fs 1000; % 采樣率 1000 Hz t (0:fs*2-1)/fs; % 2 秒時(shí)間向量 x (1 0.5*cos(2*pi*3*t)) .* cos(2*pi*50*t) ... % 50 Hz 載波受 3 Hz 調(diào)幅 0.3*sin(2*pi*8*t) ... % 8 Hz 低頻正弦 0.1*t; % 線性趨勢項(xiàng) %% 執(zhí)行 EMD 分解 [imf, residual] emd(x, Display, 1); %% 逐行繪制各 IMF 與殘差 figure tiledlayout(size(imf, 2) 1, 1) for k 1:size(imf, 2) nexttile plot(t, imf(:, k)) ylabel(sprintf(IMF%d, k)) end nexttile plot(t, residual) ylabel(殘差) xlabel(時(shí)間/s)這段代碼把采樣率設(shè)成 1000 Hz總長 2 秒第一個(gè)分量用調(diào)幅形式模擬齒輪嚙合頻率受轉(zhuǎn)頻調(diào)制的情況。emd返回的 imf 矩陣按列存放分量所以循環(huán)里取的是imf(:, k)。Display 參數(shù)打開后命令行會輸出每種停止條件對應(yīng)的迭代次數(shù)能直觀看到篩分過程收斂情況。通常第一個(gè) IMF 對應(yīng)最高頻成分往后頻率遞減殘差保留趨勢。如果某個(gè) IMF 的幅值明顯偏離原始信號量級或者殘差里還殘留振蕩成分說明信號本身太復(fù)雜或參數(shù)需要調(diào)整。3.3 讀圖判成敗數(shù)量、殘差與端點(diǎn)形態(tài)分解完成后不要急著做下一步先用眼睛掃三樣?xùn)|西IMF 數(shù)量是否合理、殘差是不是平滑的單調(diào)趨勢、每個(gè)分量的端點(diǎn)是否出現(xiàn)發(fā)散。以這段仿真信號為例理論上有三個(gè)可辨識成分但 IMF 數(shù)量不一定正好等于三因?yàn)檎{(diào)幅分量的包絡(luò)變化可能被拆成兩個(gè)模態(tài)。殘差形態(tài)是最直觀的判據(jù)。正常情況下殘差要么是直線要么是低階多項(xiàng)式形狀如果殘差還帶著明顯周期性振蕩說明高頻分量沒有被充分篩干凈。端點(diǎn)發(fā)散也很好認(rèn)波形兩端會出現(xiàn)幅度突然放大或凹進(jìn)去的“蝴蝶翅膀”形狀。遇到這種情況先檢查是不是插值方式選錯(cuò)再去考慮端點(diǎn)效應(yīng)處理。另一個(gè)容易誤判的地方是幅值量級。原始信號峰峰值大概在 1.5 左右分解出的 IMF 幅值如果動(dòng)輒超過這個(gè)范圍那多半是出現(xiàn)虛假分量了。我一般會用max(abs(imf))批量看一眼每個(gè)分量的最大幅值量級相差 10 倍以上的分量基本可以忽略。讀圖這一步看似基礎(chǔ)卻是整個(gè) EMD 流程里最省時(shí)間的一步能提前發(fā)現(xiàn) 80% 的異常分解。4. 參數(shù)調(diào)整與批量處理讓EMD在業(yè)務(wù)數(shù)據(jù)上落地4.1 五個(gè)核心參數(shù)的默認(rèn)值與調(diào)參方向官方 emd 函數(shù)提供了幾個(gè)常用參數(shù)默認(rèn)值在設(shè)計(jì)上是通用場景的折中。實(shí)際項(xiàng)目里我會根據(jù)信號特點(diǎn)做針對性調(diào)整下表是整理后的參數(shù)說明參數(shù)名默認(rèn)值作用我的調(diào)參習(xí)慣MaxNumIMF自動(dòng)限制最大 IMF 數(shù)量防止篩出過多噪聲分量知道先驗(yàn)成分?jǐn)?shù)時(shí)設(shè)為比期望數(shù)多 12SiftRelativeTolerance0.2篩分過程新建 IMF 的停止容差信號噪聲大時(shí)降到 0.1SiftMaxIterations100單次篩分最大迭代次數(shù)復(fù)雜振動(dòng)信號調(diào)到 300SiftMaxNumExtrema1殘差極值點(diǎn)數(shù)不大于此值時(shí)停止保持默認(rèn)Interpolationspline包絡(luò)插值方式可選 pchip波形有異常尖峰時(shí)換 pchipSiftRelativeTolerance 是最值得調(diào)的參數(shù)。它控制篩分過程中相鄰兩次包絡(luò)均值的變化幅度值越小篩分越嚴(yán)格IMF 越接近理想定義但會顯著增加迭代次數(shù)和計(jì)算時(shí)間。某圖像處理 Demo 的采集信號里混著明顯的白噪聲默認(rèn) 0.2 會把噪聲當(dāng)成分量篩進(jìn) IMF 的低頻段降到 0.1 之后噪聲被更多地留在殘差里。MaxNumIMF 在信號成分?jǐn)?shù)量已知時(shí)很有用。比如振動(dòng)信號里已知是轉(zhuǎn)頻、倍頻加調(diào)制邊帶一共 6 個(gè)特征分量把 MaxNumIMF 設(shè)為 8能避免算法溢出多余的噪聲分量。Interpolation 只在波形有明顯沖擊時(shí)切換成 pchip這種插值不會像三次樣條那樣產(chǎn)生大幅度過沖代價(jià)是包絡(luò)光滑度稍差篩出的 IMF 波形偏銳。這幾個(gè)參數(shù)之間會互相影響不建議一次只調(diào)一個(gè)最好按“先定 MaxNumIMF再調(diào) SiftRelativeTolerance最后換 Interpolation”的順序來。4.2 從 IMF 到時(shí)頻圖hht 的兩行調(diào)用分解完成后下一步通常是做 HHT 譜分析也就是把每個(gè) IMF 做 Hilbert 變換求瞬時(shí)頻率和瞬時(shí)幅值。官方 hht 函數(shù)可以一步到位%% 由 IMF 計(jì)算希爾伯特譜與時(shí)頻譜 [hs, f, t_hht, imf_env, imf_freq] hht(imf, fs, FrequencyLimits, [0 100]); %% 繪制 Hilbert 譜熱力圖 figure imagesc(t_hht, f, hs) axis xy xlabel(時(shí)間/s) ylabel(頻率/Hz) colormap parula colorbar這里 hht 的輸入是 imf 矩陣和采樣率 fs。FrequencyLimits 用來限制輸出的頻率范圍因?yàn)?Hilbert 變換在端點(diǎn)會產(chǎn)生虛假頻率分量限制輸出范圍可以把兩端的高頻噪聲擋在視野外。輸出里 hs 是譜密度矩陣f 是頻率軸t_hht 是時(shí)間軸imf_env 和 imf_freq 分別是瞬時(shí)幅值包絡(luò)和瞬時(shí)頻率矩陣。在多分量信號里hht 譜上能直接看到頻率隨時(shí)間的波動(dòng)。如果某條頻率線在中途發(fā)生斷裂或明顯跳變通常是模態(tài)混疊或端點(diǎn)效應(yīng)引起的這時(shí)要回上一章去調(diào)分解參數(shù)。還有個(gè)常見坑是畫圖時(shí)顏色范圍被個(gè)別極大值拉滿導(dǎo)致微弱分量看不清可以手動(dòng)設(shè)置caxis的下限比如caxis([0 0.5*max(hs(:))])。4.3 長信號的滑窗拆分與結(jié)果拼接實(shí)際采集的信號動(dòng)輒幾十萬點(diǎn)直接喂給 emd 會算得很慢。我一般會先用滑窗把長信號分成多段每段單獨(dú)分解再拼接起來分析。窗口長度取信號中最低特征頻率的 510 倍周期比如關(guān)注頻率是 10 Hz采樣率 1000 Hz窗口至少取 0.5 秒取 1 秒更穩(wěn)。%% 滑窗分解長信號 winLen 1024; % 窗口長度 hopLen 512; % 窗口步長 nWin floor((length(x) - winLen) / hopLen) 1; imfCell cell(nWin, 1); resCell cell(nWin, 1); for idx 1:nWin seg x((idx-1)*hopLen 1 : (idx-1)*hopLen winLen); [imfCell{idx}, resCell{idx}] emd(seg); end窗口重疊設(shè)置為長信號信號分析里的常見做法是重疊 50%保證邊界處不會漏掉關(guān)鍵事件。每個(gè)窗口分解出的 imf 矩陣列數(shù)不一定相同所以用 cell 數(shù)組保存。拼接前先檢查各窗口的列數(shù)分布取眾數(shù)作為統(tǒng)一分量數(shù)再對齊疊加。因?yàn)?EMD 的自適應(yīng)特性相鄰窗口的同一個(gè)物理分量可能被排到不同序號上直接按列拼接反而會串位最常見的處理是按頻率高低重排 IMF再對齊做后續(xù)統(tǒng)計(jì)。批量處理里最耗時(shí)的是循環(huán)體內(nèi)的重復(fù)包絡(luò)插值如果每段窗口都很長建議先把采樣率降下來。降采樣前用resample而不是decimate避免混疊。分段之后還能用parfor并行跑但每個(gè)工作進(jìn)程都要有一份 imf 函數(shù)的訪問權(quán)限路徑配置要在 parpool 啟動(dòng)前完成。5. 避坑與排查端點(diǎn)效應(yīng)、模態(tài)混疊與版本差異5.1 端點(diǎn)效應(yīng)兩端發(fā)散特征頻率被拉低分解結(jié)果在信號兩端出現(xiàn)明顯上翹或凹陷中間部分的波形正常端點(diǎn)處的 IMF 幅度遠(yuǎn)超信號本身幅值。這是 EMD 里最常見的現(xiàn)象我一開始以為是參數(shù)沒調(diào)好反復(fù)調(diào)整 SiftRelativeTolerance 也沒用。原因在于三次樣條包絡(luò)在端點(diǎn)只有單側(cè)有數(shù)據(jù)外推曲線自由度大包絡(luò)均值在端部嚴(yán)重偏離真實(shí)值。這個(gè)誤差會隨著篩分迭代被逐級放大最終表現(xiàn)為兩端發(fā)散。解決思路有兩類。一類是延拓法在分解前把信號首尾做鏡像延拓分解后裁掉延拓段另一類是修正法改用 pchip 插值降低過沖或者直接丟棄兩端各 5% 的采樣點(diǎn)。我現(xiàn)在的習(xí)慣是先判定關(guān)注特征是不是在信號中部如果在中間就直接裁剪端點(diǎn)如果全段都要保留就做鏡像延拓。用官方的wextend函數(shù)做對稱延拓再送進(jìn) emd比手工拼段更可控。5.2 模態(tài)混疊一個(gè) IMF 里混進(jìn)了間歇高頻某個(gè)低頻頻段的 IMF 波形在一個(gè)時(shí)間點(diǎn)突然疊加上一段高頻振蕩波形看起來像被“刺”了一下表現(xiàn)為單個(gè) IMF 內(nèi)既有趨勢又有瞬態(tài)沖擊。原因是 EMD 篩分依賴極值點(diǎn)分布間歇性出現(xiàn)的強(qiáng)沖擊會局部改變極值分布讓原本屬于高頻的成分析出到低頻箱中。這個(gè)問題的底層邏輯是包絡(luò)無法區(qū)分“混合在同一時(shí)段的不同頻率事件”本質(zhì)上是由極值點(diǎn)驅(qū)動(dòng)的自適應(yīng)篩分對間歇信號的天然敏感。解決方向有兩個(gè)層面。第一個(gè)是加白噪聲輔助的總體平均思路把多個(gè)含噪副本的分解結(jié)果做平均讓間歇頻率在各副本間被統(tǒng)計(jì)平均掉這就是 EEMD/CEEMDAN 的基本邏輯第二個(gè)是預(yù)濾波把明顯已知的頻段先用帶通濾波器切開再對各段做 EMD。我用過一段時(shí)間的 CEEMDAN 庫發(fā)現(xiàn)抑制混疊效果好不少但計(jì)算量變大會使參數(shù)迭代更耗時(shí)。5.3 大內(nèi)存與慢循環(huán)批量分解的提速方案處理 60 萬點(diǎn)的長信號時(shí)單個(gè) emd 調(diào)用跑了將近十分鐘內(nèi)存占用也一路飆到好幾 GB任務(wù)管理器里看到內(nèi)存直沖程序直接卡到無響應(yīng)。原因有兩點(diǎn)一是 emd 內(nèi)部每一輪篩分都要對整段信號做三次樣條插值迭代次數(shù)一多耗時(shí)成倍增長二是 imf 矩陣提前分配的空間和實(shí)際篩出的分量數(shù)不一致返回時(shí)可能發(fā)生重新復(fù)制和拼接內(nèi)存開銷更大。解決手段按見效順序排序先降采樣把采樣率降到特征頻率的 5 倍左右即可再分段并行按第 4.3 節(jié)的方式切窗并用parfor跑最后限制迭代次數(shù)把SiftMaxIterations從默認(rèn) 100 降到 50MaxNumIMF設(shè)一個(gè)明確的數(shù)。最有效的是分段因?yàn)?EMD 的計(jì)算復(fù)雜度接近 O(N^2)每段縮短一半總時(shí)間就能降到原來的四分之一以下。從那以后我給所有批量處理都加了“先分段、再并行、后調(diào)參”的固定順序。5.4 官方版本與第三方庫函數(shù)簽名和結(jié)果都有差異同一份腳本在一臺機(jī)器上正常跑換到另一臺機(jī)器就報(bào)參數(shù)錯(cuò)誤或者跑出來的 IMF 數(shù)量明顯不對。原因是不同來源的 emd 實(shí)現(xiàn)差異很大。官方 Signal Processing Toolbox 的函數(shù)簽名是按 Name-Value 方式傳參比如MaxNumIMF而一些第三方工具箱用位置參數(shù)有的還規(guī)定第一個(gè)輸出是 cell 數(shù)組而非矩陣。更隱蔽的是插值邏輯不同第三方庫有的用線性插值做包絡(luò)有的默認(rèn)不做端點(diǎn)延拓結(jié)果自然對不上。解決方法是先在腳本開頭做環(huán)境自檢用which emd查看路徑確認(rèn)是官方函數(shù)還是第三方腳本再用一段固定的基準(zhǔn)信號跑一次分解把期望的 IMF 數(shù)量作為閾值判斷是否兼容。切換環(huán)境時(shí)不要只復(fù)制代碼一定要把工具箱目錄或版本說明一起帶走否則復(fù)現(xiàn)結(jié)果會憑空消耗大量精力去對比差異。這套排查流程在讀者下載到的 EMD 工具箱里同樣適用裝好之后先用自帶示例跑一遍基準(zhǔn)測試再上真實(shí)數(shù)據(jù)比較靠譜。6. 驗(yàn)證IMF質(zhì)量的三個(gè)技巧相關(guān)性、正交性與瞬時(shí)頻率檢查6.1 相關(guān)性檢查IMF 與原始信號分解完之后我先算每個(gè) IMF 與原始信號的皮爾遜相關(guān)系數(shù)。物理意義在于真實(shí)分量應(yīng)該與原始信號共變相關(guān)系數(shù)太低大概率就是虛假振蕩。閾值上振動(dòng)信號里我一般取大于 0.3 才算有效分量低于這個(gè)數(shù)的可以直接剔除或作為噪聲處理。這一步整體上非?？靺s能避免把無效模態(tài)當(dāng)特征送進(jìn)分類器。6.2 正交性檢查一個(gè)數(shù)值衡量分解質(zhì)量第二個(gè)技巧是計(jì)算 IMF 之間的正交性指標(biāo)。EMD 假設(shè)各分量彼此正交但實(shí)際分解總會有少量泄漏所以用歸一化交叉積來衡量“泄漏”程度。我寫一段簡短驗(yàn)證代碼%% 計(jì)算 IMF 之間的正交性指標(biāo) ioSum 0; for m 1:size(imf, 2) for n 1:size(imf, 2) if m ~ n ioSum ioSum abs(sum(imf(:, m) .* imf(:, n))); end end end ioIndex ioSum / sum(x.^2);這個(gè)指標(biāo)越接近 0 表示分量間越正交分解越可信。我在某項(xiàng)目里看到過正交性指標(biāo)超過 0.5 的情況當(dāng)時(shí)立即回去查發(fā)現(xiàn)是信號同時(shí)存在兩個(gè)相近頻率分量EMD 把它們拆成了四五個(gè)交疊模態(tài)。這時(shí)就改用 CEEMDAN 思路重新分解指標(biāo)立刻降到 0.15 以下。6.3 瞬時(shí)頻率曲線看邊界和跳變最后一個(gè)檢查是畫每個(gè) IMF 的瞬時(shí)頻率曲線。瞬時(shí)頻率本身應(yīng)該平滑變化如果某一段曲線突然跳到高頻又跳回說明該分量內(nèi)部混疊或端部失真。我會把瞬時(shí)頻率跟原始信號的關(guān)鍵事件疊加對比確認(rèn)每條曲線都能對應(yīng)到物理意義。從那以后我每次 EMD 跑完都強(qiáng)制走一遍這三項(xiàng)驗(yàn)證不超過十分鐘能篩掉九成以上的無效分解。希望這套流程對你有用也期待你拿到工具箱之后先跑通仿真信號再上真實(shí)數(shù)據(jù)踩過的坑會少很多。本文還有配套的精品資源點(diǎn)擊獲取 MATLAB EMD工具箱實(shí)戰(zhàn)：從經(jīng)驗(yàn)?zāi)B(tài)分解到HHT時(shí)頻分析
簡介Matlab EMD經(jīng)驗(yàn)?zāi)B(tài)分解工具箱面向信號處理、故障診斷與振動(dòng)分析領(lǐng)域的研究人員與工程師用于對非線性、非平穩(wěn)信號做自適應(yīng)分解并提取有物理意義的本征模態(tài)函數(shù)分量。壓縮包共38個(gè)文件其中10個(gè)m腳本承擔(dān)核心算法與在線/局部分解入口17個(gè)c文件實(shí)現(xiàn)底層快速計(jì)算11個(gè)h頭文件提供接口聲明整體僅55KB輕量易部署。資源內(nèi)置標(biāo)準(zhǔn)EMD、復(fù)數(shù)擴(kuò)展EMD與C加速實(shí)現(xiàn)可靈活處理實(shí)信號與復(fù)信號適配科研實(shí)驗(yàn)、課程教學(xué)和工程診斷等場景。目前已有729人學(xué)習(xí)下載適合需要快速搭建EMD分析環(huán)境、二次開發(fā)或集成到現(xiàn)有Matlab流程中的用戶是一份精簡實(shí)用的信號分解工具集。1. EMD工具箱非平穩(wěn)信號分析的入場券matlab EMD工具箱解決的是經(jīng)驗(yàn)?zāi)B(tài)分解的落地問題一套能在 MATLAB 里直接跑通的函數(shù)集把非平穩(wěn)、非線性信號自適應(yīng)地拆成若干個(gè)本征模函數(shù)IMF和一條殘差。我之前幫某實(shí)驗(yàn)室處理設(shè)備啟動(dòng)瞬間的振動(dòng)數(shù)據(jù)前 0.5 秒信號直接做 FFT 根本看不出頻率成分換成 EMD 分解之后特征頻段一下子拆開了。這套資源適合三類人剛接觸經(jīng)驗(yàn)?zāi)B(tài)分解、想跳過論文偽代碼的在校學(xué)生需要做 HHT 時(shí)頻分析的信號處理工程師以及手里攢了一堆實(shí)測波形、想低成本提取特征的分析人員。下面從原理開始逐步講到參數(shù)調(diào)整、批量處理和驗(yàn)證技巧。2. 從原理到函數(shù)EMD拆出了什么怎么讀懂IMF2.1 IMF的定義與篩分循環(huán)三次樣條包絡(luò)在做什么EMD 的核心假設(shè)是任意復(fù)雜信號都能分解為若干本征模函數(shù)。IMF 要滿足兩個(gè)條件一是極值點(diǎn)數(shù)量與過零點(diǎn)數(shù)量相等或相差一個(gè)二是上下包絡(luò)的均值在任意點(diǎn)都接近零。第一條保證波形接近正弦的“局域?qū)ΨQ”第二條保證瞬時(shí)頻率有意義。兩個(gè)條件都是經(jīng)驗(yàn)性的這正是“經(jīng)驗(yàn)?zāi)B(tài)分解”名字的由來。以某跨平臺系統(tǒng)的振動(dòng)信號分析為例我的理解是 EMD 本質(zhì)上是一套自適應(yīng)篩分循環(huán)。具體流程是先找信號的局部極大值和極小值用三次樣條分別插值成上包絡(luò)和下包絡(luò)算出包絡(luò)均值原信號減去均值后得到候選分量判斷它是否滿足 IMF 條件不滿足就繼續(xù)重復(fù)篩分滿足就把它抽出去對殘差繼續(xù)執(zhí)行同樣的過程。整個(gè)過程像一層層剝洋蔥每一層都去掉一個(gè)相對高頻的振蕩成分剩下的殘差通常是單調(diào)趨勢或直流偏置。三次樣條是 EMD 的默認(rèn)插值方式因?yàn)樗亩A光滑性能讓包絡(luò)線在連接處不產(chǎn)生突兀拐點(diǎn)。不過樣條也帶來兩個(gè)副作用一是計(jì)算量大二是端點(diǎn)處的外推不受控。前者在長信號上會直接拖慢速度后者會引發(fā)下文要說的端點(diǎn)效應(yīng)。所以在理解原理的時(shí)候把“包絡(luò)生成——均值相減——迭代判定”當(dāng)作一條主線就夠了細(xì)節(jié)參數(shù)先不用背。選 EMD 而不是其它分解方法主要看三點(diǎn)。FFT 只能把整段信號折算成全局頻譜頻率隨時(shí)間的變化會被平均掉小波分解需要提前選基函數(shù)和分解層數(shù)選錯(cuò)了結(jié)果就差很遠(yuǎn)EMD 不需要預(yù)先設(shè)定基函數(shù)完全按信號自身的極值分布來篩分。代價(jià)是計(jì)算結(jié)果對極值點(diǎn)的分布和插值方式敏感同一段數(shù)據(jù)在不同參數(shù)或不同版本下可能篩出不同數(shù)量的 IMF。這是它的經(jīng)驗(yàn)性決定的不是 bug。2.2 工具箱函數(shù)地圖emd 主函數(shù)與 hht 譜分析我在項(xiàng)目里常用的是官方 Signal Processing Toolbox 自帶的 emd 函數(shù)它的標(biāo)準(zhǔn)調(diào)用方式是[imf, residual] emd(x)。imf 輸出是一個(gè)矩陣每一列是一個(gè) IMF 分量列數(shù)由信號本身的復(fù)雜度決定residual 返回最終殘差。官方版本還支持 Name-Value 形式的可選參數(shù)比如控制篩分終止條件、最大 IMF 數(shù)量、插值方式等具體會在第 4 章展開。配套的 hht 函數(shù)用于做 Hilbert-Huang 變換。它的輸入是 imf 矩陣和采樣率輸出是瞬時(shí)頻率、瞬時(shí)幅值和 Hilbert 譜。很多人會把 EMD 和 HHT 混著叫其實(shí)前者是做分解后者是對各 IMF 做 Hilbert 變換并畫時(shí)頻圖。拿到 imf 矩陣之后先看每個(gè)分量波形是否規(guī)律再用 hht 畫譜才能確認(rèn)分解結(jié)果能不能用于特征提取。除此之外工具箱還帶了自帶繪圖邏輯。比如 emd 在 Display 參數(shù)打開時(shí)會在命令行輸出篩分迭代信息配合tiledlayout可以自己拼一個(gè) IMF 逐行排列的波形圖。這個(gè)圖幾乎是驗(yàn)收分解效果的標(biāo)配看一眼 IMF 數(shù)量、每個(gè)分量的幅值量級、殘差形態(tài)基本就能判斷分解是否可信。關(guān)于這套函數(shù)更完整的用法包括各種返回值和繪圖輔助我在做完源碼包的閱讀筆記之后整理了一張函數(shù)與用途對照表放在后面的章節(jié)里。3. 安裝與第一次分解環(huán)境配置到讀圖全流程3.1 版本選擇與路徑配置別讓老腳本坑了你經(jīng)驗(yàn)?zāi)B(tài)分解在 MATLAB 里的可用性經(jīng)歷了幾個(gè)階段。比較舊的第三方腳本是用純 MATLAB 寫的函數(shù)名通常叫emd參數(shù)順序和官方版本不一樣后來官方在 Signal Processing Toolbox 里內(nèi)置了 emd 和 hht函數(shù)命名更規(guī)范底層也做了優(yōu)化。我手邊用過的 2019、2020、2023 版本都自帶這個(gè)函數(shù)基本不用額外下載安裝包。如果你拿到的是第三方工具箱壓縮包常見做法是解壓后把整個(gè)文件夾加入 MATLAB 搜索路徑。路徑配置只需要在 MATLAB 主界面執(zhí)行addpath我一般還會用savepath把路徑持久化不然重啟后又要重新加。判斷當(dāng)前環(huán)境里用的是官方版本還是第三方版本最直接的方法是which emd -all找到結(jié)果列表里第一個(gè) emd.m 的路徑看它在不在 toolbox 目錄下。這樣處理過的環(huán)境問題多了之后我的習(xí)慣是把這一條檢查寫在每個(gè)腳本開頭注釋掉也行。3.2 第一次分解腳本生成仿真信號并畫出 IMF這里直接用可復(fù)現(xiàn)的仿真信號來跑通流程。我選取了一個(gè)帶調(diào)幅特征的目標(biāo)信號疊加一個(gè)低頻正弦和一個(gè)線性趨勢項(xiàng)模擬真實(shí)測試臺里“轉(zhuǎn)子輕微不對中加基座松動(dòng)”的典型組合%% 生成仿真信號調(diào)幅分量 低頻分量 線性趨勢 fs 1000; % 采樣率 1000 Hz t (0:fs*2-1)/fs; % 2 秒時(shí)間向量 x (1 0.5*cos(2*pi*3*t)) .* cos(2*pi*50*t) ... % 50 Hz 載波受 3 Hz 調(diào)幅 0.3*sin(2*pi*8*t) ... % 8 Hz 低頻正弦 0.1*t; % 線性趨勢項(xiàng) %% 執(zhí)行 EMD 分解 [imf, residual] emd(x, Display, 1); %% 逐行繪制各 IMF 與殘差 figure tiledlayout(size(imf, 2) 1, 1) for k 1:size(imf, 2) nexttile plot(t, imf(:, k)) ylabel(sprintf(IMF%d, k)) end nexttile plot(t, residual) ylabel(殘差) xlabel(時(shí)間/s)這段代碼把采樣率設(shè)成 1000 Hz總長 2 秒第一個(gè)分量用調(diào)幅形式模擬齒輪嚙合頻率受轉(zhuǎn)頻調(diào)制的情況。emd返回的 imf 矩陣按列存放分量所以循環(huán)里取的是imf(:, k)。Display 參數(shù)打開后命令行會輸出每種停止條件對應(yīng)的迭代次數(shù)能直觀看到篩分過程收斂情況。通常第一個(gè) IMF 對應(yīng)最高頻成分往后頻率遞減殘差保留趨勢。如果某個(gè) IMF 的幅值明顯偏離原始信號量級或者殘差里還殘留振蕩成分說明信號本身太復(fù)雜或參數(shù)需要調(diào)整。3.3 讀圖判成敗數(shù)量、殘差與端點(diǎn)形態(tài)分解完成后不要急著做下一步先用眼睛掃三樣?xùn)|西IMF 數(shù)量是否合理、殘差是不是平滑的單調(diào)趨勢、每個(gè)分量的端點(diǎn)是否出現(xiàn)發(fā)散。以這段仿真信號為例理論上有三個(gè)可辨識成分但 IMF 數(shù)量不一定正好等于三因?yàn)檎{(diào)幅分量的包絡(luò)變化可能被拆成兩個(gè)模態(tài)。殘差形態(tài)是最直觀的判據(jù)。正常情況下殘差要么是直線要么是低階多項(xiàng)式形狀如果殘差還帶著明顯周期性振蕩說明高頻分量沒有被充分篩干凈。端點(diǎn)發(fā)散也很好認(rèn)波形兩端會出現(xiàn)幅度突然放大或凹進(jìn)去的“蝴蝶翅膀”形狀。遇到這種情況先檢查是不是插值方式選錯(cuò)再去考慮端點(diǎn)效應(yīng)處理。另一個(gè)容易誤判的地方是幅值量級。原始信號峰峰值大概在 1.5 左右分解出的 IMF 幅值如果動(dòng)輒超過這個(gè)范圍那多半是出現(xiàn)虛假分量了。我一般會用max(abs(imf))批量看一眼每個(gè)分量的最大幅值量級相差 10 倍以上的分量基本可以忽略。讀圖這一步看似基礎(chǔ)卻是整個(gè) EMD 流程里最省時(shí)間的一步能提前發(fā)現(xiàn) 80% 的異常分解。4. 參數(shù)調(diào)整與批量處理讓EMD在業(yè)務(wù)數(shù)據(jù)上落地4.1 五個(gè)核心參數(shù)的默認(rèn)值與調(diào)參方向官方 emd 函數(shù)提供了幾個(gè)常用參數(shù)默認(rèn)值在設(shè)計(jì)上是通用場景的折中。實(shí)際項(xiàng)目里我會根據(jù)信號特點(diǎn)做針對性調(diào)整下表是整理后的參數(shù)說明參數(shù)名默認(rèn)值作用我的調(diào)參習(xí)慣MaxNumIMF自動(dòng)限制最大 IMF 數(shù)量防止篩出過多噪聲分量知道先驗(yàn)成分?jǐn)?shù)時(shí)設(shè)為比期望數(shù)多 12SiftRelativeTolerance0.2篩分過程新建 IMF 的停止容差信號噪聲大時(shí)降到 0.1SiftMaxIterations100單次篩分最大迭代次數(shù)復(fù)雜振動(dòng)信號調(diào)到 300SiftMaxNumExtrema1殘差極值點(diǎn)數(shù)不大于此值時(shí)停止保持默認(rèn)Interpolationspline包絡(luò)插值方式可選 pchip波形有異常尖峰時(shí)換 pchipSiftRelativeTolerance 是最值得調(diào)的參數(shù)。它控制篩分過程中相鄰兩次包絡(luò)均值的變化幅度值越小篩分越嚴(yán)格IMF 越接近理想定義但會顯著增加迭代次數(shù)和計(jì)算時(shí)間。某圖像處理 Demo 的采集信號里混著明顯的白噪聲默認(rèn) 0.2 會把噪聲當(dāng)成分量篩進(jìn) IMF 的低頻段降到 0.1 之后噪聲被更多地留在殘差里。MaxNumIMF 在信號成分?jǐn)?shù)量已知時(shí)很有用。比如振動(dòng)信號里已知是轉(zhuǎn)頻、倍頻加調(diào)制邊帶一共 6 個(gè)特征分量把 MaxNumIMF 設(shè)為 8能避免算法溢出多余的噪聲分量。Interpolation 只在波形有明顯沖擊時(shí)切換成 pchip這種插值不會像三次樣條那樣產(chǎn)生大幅度過沖代價(jià)是包絡(luò)光滑度稍差篩出的 IMF 波形偏銳。這幾個(gè)參數(shù)之間會互相影響不建議一次只調(diào)一個(gè)最好按“先定 MaxNumIMF再調(diào) SiftRelativeTolerance最后換 Interpolation”的順序來。4.2 從 IMF 到時(shí)頻圖hht 的兩行調(diào)用分解完成后下一步通常是做 HHT 譜分析也就是把每個(gè) IMF 做 Hilbert 變換求瞬時(shí)頻率和瞬時(shí)幅值。官方 hht 函數(shù)可以一步到位%% 由 IMF 計(jì)算希爾伯特譜與時(shí)頻譜 [hs, f, t_hht, imf_env, imf_freq] hht(imf, fs, FrequencyLimits, [0 100]); %% 繪制 Hilbert 譜熱力圖 figure imagesc(t_hht, f, hs) axis xy xlabel(時(shí)間/s) ylabel(頻率/Hz) colormap parula colorbar這里 hht 的輸入是 imf 矩陣和采樣率 fs。FrequencyLimits 用來限制輸出的頻率范圍因?yàn)?Hilbert 變換在端點(diǎn)會產(chǎn)生虛假頻率分量限制輸出范圍可以把兩端的高頻噪聲擋在視野外。輸出里 hs 是譜密度矩陣f 是頻率軸t_hht 是時(shí)間軸imf_env 和 imf_freq 分別是瞬時(shí)幅值包絡(luò)和瞬時(shí)頻率矩陣。在多分量信號里hht 譜上能直接看到頻率隨時(shí)間的波動(dòng)。如果某條頻率線在中途發(fā)生斷裂或明顯跳變通常是模態(tài)混疊或端點(diǎn)效應(yīng)引起的這時(shí)要回上一章去調(diào)分解參數(shù)。還有個(gè)常見坑是畫圖時(shí)顏色范圍被個(gè)別極大值拉滿導(dǎo)致微弱分量看不清可以手動(dòng)設(shè)置caxis的下限比如caxis([0 0.5*max(hs(:))])。4.3 長信號的滑窗拆分與結(jié)果拼接實(shí)際采集的信號動(dòng)輒幾十萬點(diǎn)直接喂給 emd 會算得很慢。我一般會先用滑窗把長信號分成多段每段單獨(dú)分解再拼接起來分析。窗口長度取信號中最低特征頻率的 510 倍周期比如關(guān)注頻率是 10 Hz采樣率 1000 Hz窗口至少取 0.5 秒取 1 秒更穩(wěn)。%% 滑窗分解長信號 winLen 1024; % 窗口長度 hopLen 512; % 窗口步長 nWin floor((length(x) - winLen) / hopLen) 1; imfCell cell(nWin, 1); resCell cell(nWin, 1); for idx 1:nWin seg x((idx-1)*hopLen 1 : (idx-1)*hopLen winLen); [imfCell{idx}, resCell{idx}] emd(seg); end窗口重疊設(shè)置為長信號信號分析里的常見做法是重疊 50%保證邊界處不會漏掉關(guān)鍵事件。每個(gè)窗口分解出的 imf 矩陣列數(shù)不一定相同所以用 cell 數(shù)組保存。拼接前先檢查各窗口的列數(shù)分布取眾數(shù)作為統(tǒng)一分量數(shù)再對齊疊加。因?yàn)?EMD 的自適應(yīng)特性相鄰窗口的同一個(gè)物理分量可能被排到不同序號上直接按列拼接反而會串位最常見的處理是按頻率高低重排 IMF再對齊做后續(xù)統(tǒng)計(jì)。批量處理里最耗時(shí)的是循環(huán)體內(nèi)的重復(fù)包絡(luò)插值如果每段窗口都很長建議先把采樣率降下來。降采樣前用resample而不是decimate避免混疊。分段之后還能用parfor并行跑但每個(gè)工作進(jìn)程都要有一份 imf 函數(shù)的訪問權(quán)限路徑配置要在 parpool 啟動(dòng)前完成。5. 避坑與排查端點(diǎn)效應(yīng)、模態(tài)混疊與版本差異5.1 端點(diǎn)效應(yīng)兩端發(fā)散特征頻率被拉低分解結(jié)果在信號兩端出現(xiàn)明顯上翹或凹陷中間部分的波形正常端點(diǎn)處的 IMF 幅度遠(yuǎn)超信號本身幅值。這是 EMD 里最常見的現(xiàn)象我一開始以為是參數(shù)沒調(diào)好反復(fù)調(diào)整 SiftRelativeTolerance 也沒用。原因在于三次樣條包絡(luò)在端點(diǎn)只有單側(cè)有數(shù)據(jù)外推曲線自由度大包絡(luò)均值在端部嚴(yán)重偏離真實(shí)值。這個(gè)誤差會隨著篩分迭代被逐級放大最終表現(xiàn)為兩端發(fā)散。解決思路有兩類。一類是延拓法在分解前把信號首尾做鏡像延拓分解后裁掉延拓段另一類是修正法改用 pchip 插值降低過沖或者直接丟棄兩端各 5% 的采樣點(diǎn)。我現(xiàn)在的習(xí)慣是先判定關(guān)注特征是不是在信號中部如果在中間就直接裁剪端點(diǎn)如果全段都要保留就做鏡像延拓。用官方的wextend函數(shù)做對稱延拓再送進(jìn) emd比手工拼段更可控。5.2 模態(tài)混疊一個(gè) IMF 里混進(jìn)了間歇高頻某個(gè)低頻頻段的 IMF 波形在一個(gè)時(shí)間點(diǎn)突然疊加上一段高頻振蕩波形看起來像被“刺”了一下表現(xiàn)為單個(gè) IMF 內(nèi)既有趨勢又有瞬態(tài)沖擊。原因是 EMD 篩分依賴極值點(diǎn)分布間歇性出現(xiàn)的強(qiáng)沖擊會局部改變極值分布讓原本屬于高頻的成分析出到低頻箱中。這個(gè)問題的底層邏輯是包絡(luò)無法區(qū)分“混合在同一時(shí)段的不同頻率事件”本質(zhì)上是由極值點(diǎn)驅(qū)動(dòng)的自適應(yīng)篩分對間歇信號的天然敏感。解決方向有兩個(gè)層面。第一個(gè)是加白噪聲輔助的總體平均思路把多個(gè)含噪副本的分解結(jié)果做平均讓間歇頻率在各副本間被統(tǒng)計(jì)平均掉這就是 EEMD/CEEMDAN 的基本邏輯第二個(gè)是預(yù)濾波把明顯已知的頻段先用帶通濾波器切開再對各段做 EMD。我用過一段時(shí)間的 CEEMDAN 庫發(fā)現(xiàn)抑制混疊效果好不少但計(jì)算量變大會使參數(shù)迭代更耗時(shí)。5.3 大內(nèi)存與慢循環(huán)批量分解的提速方案處理 60 萬點(diǎn)的長信號時(shí)單個(gè) emd 調(diào)用跑了將近十分鐘內(nèi)存占用也一路飆到好幾 GB任務(wù)管理器里看到內(nèi)存直沖程序直接卡到無響應(yīng)。原因有兩點(diǎn)一是 emd 內(nèi)部每一輪篩分都要對整段信號做三次樣條插值迭代次數(shù)一多耗時(shí)成倍增長二是 imf 矩陣提前分配的空間和實(shí)際篩出的分量數(shù)不一致返回時(shí)可能發(fā)生重新復(fù)制和拼接內(nèi)存開銷更大。解決手段按見效順序排序先降采樣把采樣率降到特征頻率的 5 倍左右即可再分段并行按第 4.3 節(jié)的方式切窗并用parfor跑最后限制迭代次數(shù)把SiftMaxIterations從默認(rèn) 100 降到 50MaxNumIMF設(shè)一個(gè)明確的數(shù)。最有效的是分段因?yàn)?EMD 的計(jì)算復(fù)雜度接近 O(N^2)每段縮短一半總時(shí)間就能降到原來的四分之一以下。從那以后我給所有批量處理都加了“先分段、再并行、后調(diào)參”的固定順序。5.4 官方版本與第三方庫函數(shù)簽名和結(jié)果都有差異同一份腳本在一臺機(jī)器上正常跑換到另一臺機(jī)器就報(bào)參數(shù)錯(cuò)誤或者跑出來的 IMF 數(shù)量明顯不對。原因是不同來源的 emd 實(shí)現(xiàn)差異很大。官方 Signal Processing Toolbox 的函數(shù)簽名是按 Name-Value 方式傳參比如MaxNumIMF而一些第三方工具箱用位置參數(shù)有的還規(guī)定第一個(gè)輸出是 cell 數(shù)組而非矩陣。更隱蔽的是插值邏輯不同第三方庫有的用線性插值做包絡(luò)有的默認(rèn)不做端點(diǎn)延拓結(jié)果自然對不上。解決方法是先在腳本開頭做環(huán)境自檢用which emd查看路徑確認(rèn)是官方函數(shù)還是第三方腳本再用一段固定的基準(zhǔn)信號跑一次分解把期望的 IMF 數(shù)量作為閾值判斷是否兼容。切換環(huán)境時(shí)不要只復(fù)制代碼一定要把工具箱目錄或版本說明一起帶走否則復(fù)現(xiàn)結(jié)果會憑空消耗大量精力去對比差異。這套排查流程在讀者下載到的 EMD 工具箱里同樣適用裝好之后先用自帶示例跑一遍基準(zhǔn)測試再上真實(shí)數(shù)據(jù)比較靠譜。6. 驗(yàn)證IMF質(zhì)量的三個(gè)技巧相關(guān)性、正交性與瞬時(shí)頻率檢查6.1 相關(guān)性檢查IMF 與原始信號分解完之后我先算每個(gè) IMF 與原始信號的皮爾遜相關(guān)系數(shù)。物理意義在于真實(shí)分量應(yīng)該與原始信號共變相關(guān)系數(shù)太低大概率就是虛假振蕩。閾值上振動(dòng)信號里我一般取大于 0.3 才算有效分量低于這個(gè)數(shù)的可以直接剔除或作為噪聲處理。這一步整體上非?？靺s能避免把無效模態(tài)當(dāng)特征送進(jìn)分類器。6.2 正交性檢查一個(gè)數(shù)值衡量分解質(zhì)量第二個(gè)技巧是計(jì)算 IMF 之間的正交性指標(biāo)。EMD 假設(shè)各分量彼此正交但實(shí)際分解總會有少量泄漏所以用歸一化交叉積來衡量“泄漏”程度。我寫一段簡短驗(yàn)證代碼%% 計(jì)算 IMF 之間的正交性指標(biāo) ioSum 0; for m 1:size(imf, 2) for n 1:size(imf, 2) if m ~ n ioSum ioSum abs(sum(imf(:, m) .* imf(:, n))); end end end ioIndex ioSum / sum(x.^2);這個(gè)指標(biāo)越接近 0 表示分量間越正交分解越可信。我在某項(xiàng)目里看到過正交性指標(biāo)超過 0.5 的情況當(dāng)時(shí)立即回去查發(fā)現(xiàn)是信號同時(shí)存在兩個(gè)相近頻率分量EMD 把它們拆成了四五個(gè)交疊模態(tài)。這時(shí)就改用 CEEMDAN 思路重新分解指標(biāo)立刻降到 0.15 以下。6.3 瞬時(shí)頻率曲線看邊界和跳變最后一個(gè)檢查是畫每個(gè) IMF 的瞬時(shí)頻率曲線。瞬時(shí)頻率本身應(yīng)該平滑變化如果某一段曲線突然跳到高頻又跳回說明該分量內(nèi)部混疊或端部失真。我會把瞬時(shí)頻率跟原始信號的關(guān)鍵事件疊加對比確認(rèn)每條曲線都能對應(yīng)到物理意義。從那以后我每次 EMD 跑完都強(qiáng)制走一遍這三項(xiàng)驗(yàn)證不超過十分鐘能篩掉九成以上的無效分解。希望這套流程對你有用也期待你拿到工具箱之后先跑通仿真信號再上真實(shí)數(shù)據(jù)踩過的坑會少很多。本文還有配套的精品資源點(diǎn)擊獲取