簡(jiǎn)介這份資源面向具備一定Python基礎(chǔ)、希望上手非線性信號(hào)處理的開(kāi)發(fā)者與算法學(xué)習(xí)者提供變分模態(tài)分解VMD的可運(yùn)行代碼與配套測(cè)試用例可用于聲音分析、故障診斷、金融時(shí)序等場(chǎng)景的信號(hào)分解實(shí)踐。壓縮包共2個(gè)文件包含1個(gè)py腳本與1個(gè)txt說(shuō)明整體約1KB體量輕巧py文件承載VMD分解函數(shù)與測(cè)試邏輯txt則給出依賴庫(kù)的安裝指引便于快速搭建numpy、scipy環(huán)境。目前已有5452人學(xué)習(xí)下載熱度較高。通過(guò)研讀代碼讀者可理解模態(tài)個(gè)數(shù)K、迭代次數(shù)與尺度因子α等參數(shù)如何影響分解效果掌握將非平穩(wěn)信號(hào)拆解為多個(gè)簡(jiǎn)諧模態(tài)、并檢驗(yàn)?zāi)B(tài)正交性與頻譜特性的完整流程也能借鑒其測(cè)試思路把VMD作為預(yù)處理環(huán)節(jié)接入機(jī)器學(xué)習(xí)或后端數(shù)據(jù)分析管道提升時(shí)序特征提取與預(yù)測(cè)精度。1. 變分模態(tài)分解到底解決了什么信號(hào)難題如果你做過(guò)旋轉(zhuǎn)機(jī)械故障診斷、腦電節(jié)律提取或者電力諧波分析大概率被同一個(gè)問(wèn)題折磨過(guò)一段非平穩(wěn)信號(hào)里混著好幾個(gè)中心頻率不同的分量想拆開(kāi)看用傅里葉變換會(huì)糊成一片用經(jīng)驗(yàn)?zāi)B(tài)分解EMD又經(jīng)常模態(tài)混疊端點(diǎn)還發(fā)散。變分模態(tài)分解VMD就是沖著這個(gè)痛點(diǎn)來(lái)的。它把「分解」這件事從遞歸篩分改寫成變分求解先假設(shè)信號(hào)由 K 個(gè)有限帶寬的模態(tài)疊加而成再讓每個(gè)模態(tài)圍繞自己的中心頻率最緊湊同時(shí)所有模態(tài)加起來(lái)能重構(gòu)原信號(hào)。這套思路的好處是抗噪、抗混疊模態(tài)數(shù)和中心頻率都能提前約束結(jié)果比 EMD 穩(wěn)定得多。這篇筆記不講空理論直接給能跑的 python 代碼、測(cè)試用例以及我調(diào)參時(shí)踩過(guò)的坑適合已經(jīng)會(huì)寫 python、想把這套方法落到自己數(shù)據(jù)上的工程師。2. VMD 的數(shù)學(xué)骨架與 python 實(shí)現(xiàn)選型2.1 從維納濾波到 ADMMVMD 為什么能拆得干凈VMD 的核心是把每個(gè)模態(tài)定義成一個(gè)調(diào)幅調(diào)頻信號(hào)寫成 $u_k(t) A_k(t)\cos(\phi_k(t))$然后對(duì)每個(gè)模態(tài)做希爾伯特變換得到解析信號(hào)再乘以 $e^{-j\omega_k t}$ 把頻譜搬到基帶。搬完之后用 L2 范數(shù)平方去度量這個(gè)基帶信號(hào)的帶寬帶寬越小說(shuō)明這個(gè)模態(tài)越「純」。于是整個(gè)分解變成一個(gè)帶約束的優(yōu)化問(wèn)題所有模態(tài)之和等于原信號(hào)。約束不好直接解就引入二次懲罰項(xiàng)和拉格朗日乘子轉(zhuǎn)成無(wú)約束問(wèn)題再用交替方向乘子法ADMM迭代。每一輪迭代里模態(tài)、中心頻率、拉格朗日乘子輪流更新模態(tài)在頻域有閉式解中心頻率是模態(tài)頻譜的重心。這就是為什么 VMD 比 EMD 穩(wěn)它不是靠極值點(diǎn)篩出來(lái)的而是靠頻域迭代收斂出來(lái)的。理解這一點(diǎn)對(duì)調(diào)參很關(guān)鍵。懲罰因子 α 控制的是帶寬約束的強(qiáng)弱α 越大每個(gè)模態(tài)被逼得越窄容易把有用分量切碎α 越小模態(tài)越寬容易混疊。模態(tài)數(shù) K 直接決定拆成幾份K 給小了會(huì)欠分解兩個(gè)分量擠在一個(gè)模態(tài)里K 給大了會(huì)過(guò)分解同一個(gè)物理分量被劈成兩半中心頻率還會(huì)靠得很近。這兩個(gè)參數(shù)是 VMD 的命門后面會(huì)專門講怎么定。2.2 用 vmdpy 跑通第一段信號(hào)最小可復(fù)現(xiàn)腳本python 里現(xiàn)成的 VMD 實(shí)現(xiàn)不多常見(jiàn)做法是用vmdpy這個(gè)包它把 ADMM 迭代封裝好了接口干凈。先裝依賴再跑一段合成信號(hào)驗(yàn)證。pip install vmdpy numpy matplotlibimport numpy as np import matplotlib.pyplot as plt from vmdpy import VMD # 構(gòu)造合成信號(hào)三個(gè)不同頻率分量 噪聲 fs 1000 # 采樣率 Hz t np.arange(0, 1, 1/fs) # 1 秒時(shí)間軸 f1, f2, f3 30, 80, 150 # 三個(gè)分量中心頻率 sig (np.cos(2*np.pi*f1*t) 0.6*np.cos(2*np.pi*f2*t) 0.4*np.cos(2*np.pi*f3*t)) sig 0.1 * np.random.randn(len(t)) # 加一點(diǎn)高斯噪聲 # VMD 參數(shù) alpha 2000 # 帶寬約束中等強(qiáng)度 tau 0 # 噪聲容限0 表示無(wú)松弛 K 3 # 模態(tài)數(shù)已知三個(gè)分量 DC 0 # 不含直流分量 init 1 # 中心頻率均勻初始化 tol 1e-7 # 收斂容差 u, u_hat, omega VMD(sig, alpha, tau, K, DC, init, tol) # u: 各模態(tài)時(shí)域波形形狀 (K, N) # u_hat: 各模態(tài)頻域譜 # omega: 每次迭代的中心頻率最后一行是收斂結(jié)果 print(收斂后中心頻率(Hz):, omega[-1] * fs) plt.figure(figsize(10, 6)) for i in range(K): plt.subplot(K, 1, i1) plt.plot(t, u[i]) plt.title(fMode {i1}) plt.tight_layout() plt.show()這段代碼的邏輯是先造一個(gè)頻率成分已知的信號(hào)方便對(duì)照分解結(jié)果對(duì)不對(duì)。VMD返回三個(gè)東西u是時(shí)域模態(tài)u_hat是頻域omega記錄每輪迭代的中心頻率。最后一行omega[-1] * fs把歸一化頻率還原成 Hz正常應(yīng)該接近 30、80、150。參數(shù)上alpha2000是個(gè)常用起點(diǎn)K3是因?yàn)槲覀兪孪戎婪至繑?shù)init1表示中心頻率均勻鋪開(kāi)初始化比全零初始化收斂更穩(wěn)。跑完如果三個(gè)模態(tài)的頻譜峰值分別落在 30、80、150 附近說(shuō)明這套流程通了。2.3 參數(shù)怎么定K 和 alpha 的實(shí)操判據(jù)新手最容易卡在 K 和 alpha 上。我的經(jīng)驗(yàn)是先定 K 再調(diào) alpha。定 K 有個(gè)便宜辦法對(duì)信號(hào)做 FFT看主峰個(gè)數(shù)主峰幾個(gè) K 就取幾個(gè)再往上加一兩個(gè)試。更嚴(yán)謹(jǐn)?shù)淖龇ㄊ强捶纸夂蟾髂B(tài)的中心頻率有沒(méi)有靠得太近如果兩個(gè)模態(tài)中心頻率差不到主頻的 10%基本就是過(guò)分解了。另一個(gè)判據(jù)是重構(gòu)誤差把所有模態(tài)加起來(lái)和原信號(hào)比誤差應(yīng)該很小如果 K 加大后誤差沒(méi)明顯下降說(shuō)明加多了。alpha 的取值和采樣率、信號(hào)帶寬有關(guān)。采樣率高、分量窄的時(shí)候 alpha 要往大調(diào)常見(jiàn)范圍是 1000 到 5000。判斷標(biāo)準(zhǔn)是看模態(tài)頻譜的帶寬如果模態(tài)譜拖得很寬、和鄰居重疊alpha 調(diào)大如果模態(tài)譜被切得只剩一根尖峰、時(shí)域波形失真alpha 調(diào)小。我一般會(huì)寫個(gè)小循環(huán)掃幾個(gè) alpha畫重構(gòu)誤差和模態(tài)帶寬的曲線挑拐點(diǎn)。for a in [500, 1000, 2000, 4000, 8000]: u_tmp, _, om_tmp VMD(sig, a, 0, 3, 0, 1, 1e-7) recon u_tmp.sum(axis0) err np.linalg.norm(sig - recon) / np.linalg.norm(sig) print(falpha{a}, 重構(gòu)相對(duì)誤差{err:.4f}, 中心頻率{om_tmp[-1]*fs})這段掃描代碼能幫你快速看出 alpha 對(duì)結(jié)果的影響。注意重構(gòu)誤差不是越小越好過(guò)小的 alpha 會(huì)讓模態(tài)變寬、誤差也小但分解沒(méi)意義。要結(jié)合中心頻率是否穩(wěn)定一起看。3. 測(cè)試用例怎么寫從合成信號(hào)到真實(shí)數(shù)據(jù)的驗(yàn)證3.1 合成信號(hào)測(cè)試用例已知答案的回歸驗(yàn)證給 VMD 寫測(cè)試用例第一層是合成信號(hào)因?yàn)榇鸢敢阎墚?dāng)回歸測(cè)試用。核心斷言有三條分解出的模態(tài)數(shù)等于 K、中心頻率和設(shè)定值偏差在容差內(nèi)、重構(gòu)誤差低于閾值。用 pytest 組織方便每次改參數(shù)后重跑。import numpy as np import pytest from vmdpy import VMD def make_signal(fs1000, dur1.0, freqs(30, 80, 150), noise0.1): t np.arange(0, dur, 1/fs) sig sum(np.cos(2*np.pi*f*t) for f in freqs) rng np.random.default_rng(42) # 固定隨機(jī)種子保證可復(fù)現(xiàn) sig sig noise * rng.standard_normal(len(t)) return t, sig def test_vmd_recovers_known_frequencies(): fs 1000 freqs (30, 80, 150) t, sig make_signal(fsfs, freqsfreqs) u, _, omega VMD(sig, 2000, 0, 3, 0, 1, 1e-7) est np.sort(omega[-1] * fs) # 中心頻率偏差不超過(guò) 5 Hz assert np.allclose(est, sorted(freqs), atol5), f頻率偏差過(guò)大: {est} def test_vmd_reconstruction_error(): t, sig make_signal() u, _, _ VMD(sig, 2000, 0, 3, 0, 1, 1e-7) recon u.sum(axis0) err np.linalg.norm(sig - recon) / np.linalg.norm(sig) assert err 0.05, f重構(gòu)誤差過(guò)大: {err}這里兩個(gè)用例分別盯頻率恢復(fù)和重構(gòu)精度。default_rng(42)固定種子是關(guān)鍵否則每次噪聲不同斷言閾值不好定。atol5是頻率容差采樣率 1000 Hz 下 5 Hz 已經(jīng)夠松如果實(shí)際偏差超過(guò)這個(gè)數(shù)說(shuō)明 K 或 alpha 選錯(cuò)了。重構(gòu)誤差閾值 0.05 是經(jīng)驗(yàn)值噪聲大時(shí)可以放寬到 0.1。3.2 邊界與異常用例短信號(hào)、單頻、強(qiáng)噪合成信號(hào)跑通不代表代碼健壯。真實(shí)數(shù)據(jù)里經(jīng)常遇到信號(hào)很短、只有一個(gè)頻率成分、信噪比很低的情況這些都要單獨(dú)寫用例。短信號(hào)的問(wèn)題是 ADMM 迭代次數(shù)不夠中心頻率還沒(méi)收斂就停了單頻信號(hào)如果 K 設(shè)成 3會(huì)強(qiáng)行拆出三個(gè)模態(tài)其中兩個(gè)是噪聲這時(shí)候要斷言「多余模態(tài)的能量占比很低」。def test_vmd_short_signal(): # 只有 100 個(gè)采樣點(diǎn)驗(yàn)證不報(bào)錯(cuò)且能返回 K 個(gè)模態(tài) t, sig make_signal(fs1000, dur0.1) u, _, _ VMD(sig, 2000, 0, 3, 0, 1, 1e-7) assert u.shape[0] 3 def test_vmd_single_tone_extra_modes_low_energy(): # 單頻信號(hào)K3多出的模態(tài)能量應(yīng)很低 fs 1000 t np.arange(0, 1, 1/fs) sig np.cos(2*np.pi*50*t) 0.05*np.random.randn(len(t)) u, _, _ VMD(sig, 2000, 0, 3, 0, 1, 1e-7) energy np.sum(u**2, axis1) energy_ratio np.sort(energy) / energy.sum() # 最小的兩個(gè)模態(tài)能量占比之和應(yīng)小于 0.2 assert energy_ratio[:2].sum() 0.2, f多余模態(tài)能量過(guò)高: {energy_ratio}短信號(hào)用例只斷言形狀因?yàn)?100 點(diǎn)下頻率分辨率只有 10 Hz硬卡頻率會(huì)誤報(bào)。單頻用例用能量占比判斷過(guò)分解比看頻率更穩(wěn)。這兩個(gè)用例能擋住大部分「參數(shù)亂設(shè)也能跑但結(jié)果沒(méi)意義」的情況。3.3 真實(shí)數(shù)據(jù)接入以軸承振動(dòng)信號(hào)為例真實(shí)數(shù)據(jù)接入時(shí)第一件事是確認(rèn)采樣率和單位。軸承振動(dòng)信號(hào)常見(jiàn)采樣率是 12 kHz 或 25.6 kHz直接丟進(jìn) VMD 之前最好先做歸一化否則 alpha 的量綱對(duì)不上。我一般先把信號(hào)減均值、除標(biāo)準(zhǔn)差再送進(jìn) VMD。另外真實(shí)信號(hào)往往有趨勢(shì)項(xiàng)低頻漂移會(huì)干擾第一個(gè)模態(tài)可以先高通濾波或者去趨勢(shì)。from scipy.signal import detrend def preprocess(x): x detrend(x) # 去線性趨勢(shì) x (x - x.mean()) / x.std() # 標(biāo)準(zhǔn)化 return x # 假設(shè) raw 是從 csv 讀入的一維振動(dòng)信號(hào) # raw np.loadtxt(bearing.csv, delimiter,)[:, 1] # sig preprocess(raw) # u, _, omega VMD(sig, 2000, 0, 5, 0, 1, 1e-7)預(yù)處理這步別省。我見(jiàn)過(guò)直接把原始振動(dòng)信號(hào)丟進(jìn)去的結(jié)果第一個(gè)模態(tài)全是趨勢(shì)項(xiàng)后面幾個(gè)模態(tài)擠在一起調(diào)半天參數(shù)都沒(méi)用。標(biāo)準(zhǔn)化之后 alpha 的取值也更好遷移同一套參數(shù)換個(gè)數(shù)據(jù)集不至于完全失效。4. 避坑與排查VMD 落地時(shí)最容易翻車的五件事4.1 模態(tài)數(shù) K 設(shè)大了中心頻率擠在一起現(xiàn)象分解出 5 個(gè)模態(tài)但打印中心頻率發(fā)現(xiàn)有兩個(gè)只差幾 Hz時(shí)域波形看著像同一個(gè)分量被劈開(kāi)。原因K 超過(guò)信號(hào)實(shí)際分量數(shù)ADMM 會(huì)把一個(gè)寬帶模態(tài)拆成兩個(gè)窄的。解決先做 FFT 數(shù)主峰K 從主峰數(shù)開(kāi)始試或者畫中心頻率隨 K 變化的曲線找頻率不再明顯分開(kāi)的拐點(diǎn)。我一般會(huì)跑 K2 到 8看中心頻率分布出現(xiàn)兩個(gè)模態(tài)頻率差小于主頻 10% 就停。4.2 alpha 太小導(dǎo)致模態(tài)混疊現(xiàn)象模態(tài)時(shí)域波形里能看到明顯的拍頻頻譜上兩個(gè)峰連在一起。原因alpha 太小帶寬約束太弱ADMM 允許模態(tài)變寬去擬合更多成分。解決把 alpha 往上調(diào)常見(jiàn)從 2000 起每次翻倍試。注意 alpha 太大又會(huì)讓模態(tài)譜變成尖峰、時(shí)域失真所以要兩邊夾。判斷標(biāo)準(zhǔn)是模態(tài)頻譜的 3dB 帶寬應(yīng)該和該分量的物理帶寬量級(jí)一致。4.3 信號(hào)沒(méi)預(yù)處理低頻趨勢(shì)吃掉第一個(gè)模態(tài)現(xiàn)象第一個(gè)模態(tài)幾乎是一條緩慢變化的曲線后面模態(tài)才是有用的振蕩分量。原因原始信號(hào)有直流或線性趨勢(shì)VMD 把它當(dāng)成一個(gè)低頻模態(tài)。解決進(jìn) VMD 前做去趨勢(shì)和標(biāo)準(zhǔn)化。如果趨勢(shì)是非線性的先用高通濾波或者經(jīng)驗(yàn)方法去掉。這一步不做后面參數(shù)怎么調(diào)都是白費(fèi)。4.4 短信號(hào)迭代不收斂結(jié)果每次不一樣現(xiàn)象信號(hào)長(zhǎng)度只有幾百點(diǎn)跑兩次 VMD 結(jié)果中心頻率差很多。原因ADMM 迭代次數(shù)上限到了還沒(méi)收斂而初始化對(duì)短信號(hào)影響大。解決短信號(hào)要么加長(zhǎng)數(shù)據(jù)窗要么固定init參數(shù)并放寬tol同時(shí)把最大迭代次數(shù)調(diào)大。如果還是不穩(wěn)考慮對(duì)信號(hào)做鏡像延拓再分解取中間段結(jié)果。4.5 把 VMD 當(dāng)萬(wàn)能分解器什么信號(hào)都上現(xiàn)象對(duì)沖擊成分為主的信號(hào)做 VMD分解結(jié)果沒(méi)有物理意義。原因VMD 假設(shè)模態(tài)是窄帶準(zhǔn)平穩(wěn)的沖擊信號(hào)是寬帶瞬態(tài)不符合假設(shè)。解決先判斷信號(hào)類型窄帶多分量用 VMD寬帶瞬態(tài)用其他方法或者先做包絡(luò)分析。工具沒(méi)有好壞只有合不合適硬套只會(huì)得到一堆沒(méi)法解釋的模態(tài)。5. 進(jìn)階技巧用中心頻率穩(wěn)定性反推最優(yōu) K調(diào) K 最煩的是沒(méi)有客觀標(biāo)準(zhǔn)。我后來(lái)固定用一個(gè)技巧對(duì)同一個(gè)信號(hào)K 從 2 掃到 10每個(gè) K 跑多次改變初始化或加微小擾動(dòng)記錄每次收斂后的中心頻率看哪些頻率在所有 K 下都穩(wěn)定出現(xiàn)。真正屬于信號(hào)的頻率會(huì)反復(fù)出現(xiàn)過(guò)分解產(chǎn)生的虛假頻率會(huì)隨 K 跳來(lái)跳去。這個(gè)方法比單看一次分解結(jié)果靠譜得多。def stable_frequencies(sig, fs, k_rangerange(2, 11), repeats5): from collections import defaultdict freq_count defaultdict(int) for K in k_range: for r in range(repeats): # 每次加極小擾動(dòng)模擬初始化差異 noise 1e-6 * np.random.randn(len(sig)) _, _, omega VMD(sig noise, 2000, 0, K, 0, 1, 1e-7) for f in omega[-1] * fs: # 把頻率歸到 5 Hz 的桶里統(tǒng)計(jì) bucket round(f / 5) * 5 freq_count[bucket] 1 # 出現(xiàn)次數(shù)多的頻率就是穩(wěn)定分量 stable sorted([f for f, c in freq_count.items() if c repeats], keylambda x: -freq_count[x]) return stable, freq_count # stable, counts stable_frequencies(sig, fs) # print(穩(wěn)定中心頻率:, stable)這段代碼把頻率按 5 Hz 分桶統(tǒng)計(jì)出現(xiàn)次數(shù)達(dá)到重復(fù)次數(shù)的頻率才算穩(wěn)定。repeats5是每個(gè) K 跑 5 次k_range覆蓋 2 到 10。跑完看stable列表前幾個(gè)就是信號(hào)的真實(shí)分量頻率K 取這個(gè)個(gè)數(shù)就行。這個(gè)方法對(duì)噪聲有一定容忍度因?yàn)樵肼暜a(chǎn)生的虛假頻率不會(huì)穩(wěn)定復(fù)現(xiàn)。參數(shù)上分桶寬度 5 Hz 要按采樣率調(diào)采樣率低的時(shí)候桶要寬一點(diǎn)否則同一個(gè)頻率會(huì)被分到相鄰?fù)袄?。repeats越大越穩(wěn)但越慢一般 5 次夠用。如果信號(hào)很長(zhǎng)可以只取一段做這個(gè)掃描定好 K 再對(duì)全長(zhǎng)信號(hào)分解。我自己現(xiàn)在的習(xí)慣是拿到新信號(hào)先跑這個(gè)穩(wěn)定性掃描把 K 定下來(lái)再單獨(dú)調(diào) alpha最后寫兩個(gè) pytest 用例把參數(shù)鎖住。這樣換數(shù)據(jù)、換同事接手結(jié)果都能復(fù)現(xiàn)。VMD 不難難的是參數(shù)有依據(jù)、結(jié)果可驗(yàn)證把這兩件事做扎實(shí)它就能穩(wěn)定產(chǎn)出可解釋的模態(tài)。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取