綜述

作者：aresmiki
鏈接：https://www.zhihu.com/question/23701194/answer/167005497
來源：知乎
著作權(quán)歸作者所有。商業(yè)轉(zhuǎn)載請(qǐng)聯(lián)系作者獲得授權(quán)，非商業(yè)轉(zhuǎn)載請(qǐng)注明出處。
非平穩(wěn)信號(hào)處理應(yīng)該是現(xiàn)在信號(hào)處理技術(shù)最新的也是最熱的研究方向了，信號(hào)處理方法從最早的時(shí)域統(tǒng)計(jì)到傅里葉變換的頻域分析，是人們認(rèn)識(shí)信號(hào)本質(zhì)的一次巨大飛躍，給了分析人員換個(gè)角度看世界的方法，這個(gè)時(shí)期傅里葉變換在調(diào)和分析，諧波分析等領(lǐng)域得到的巨大發(fā)展，當(dāng)然工程人員將傅里葉變換也迅速運(yùn)用的工程運(yùn)用中，得到了巨大發(fā)展，特別是在機(jī)械故障診斷和 功能性磁共振成像fMRI 中。但是現(xiàn)實(shí)卻不是理論那么美好，人們發(fā)現(xiàn)，傅里葉變換作為一個(gè)全局變換，天然的少了另一個(gè)維度，如果將時(shí)間域信號(hào)比作一個(gè)平面中的物體的話，那么頻域信號(hào)也同樣是一個(gè)平面中的物體，只是給我們換了一個(gè)角度而已，而人們總是對(duì)三維世界的物體更具有直觀了解，平面的東西始終是不具體不形象，信號(hào)一樣，工程人員總是想知道信號(hào)有哪些頻率，且這些頻率在何時(shí)產(chǎn)生，而這個(gè)需求就給分析方法提出了一個(gè)要求，必須多一個(gè)維度，也就是頻域變換需要保留時(shí)間信息，聰明的Gabor第一次提出了窗口傅里葉變換，這個(gè)信號(hào)分析帶來了時(shí)間尺度，也第一次有了時(shí)頻分析的概念，當(dāng)然其理論就不介紹了，而其優(yōu)點(diǎn)當(dāng)然就是同時(shí)給了我們時(shí)間和頻率的信息。科學(xué)的魅力就是，永遠(yuǎn)沒有完美的理論，窗口福利葉變換的窗口如何選擇成為了一個(gè)難題，選的太大，時(shí)間分辨率太大，選的太小，頻率分辨率太大，如何能在該大的時(shí)候大，該小的時(shí)候小呢？偉大的Morlet告訴我們，傅里葉變換本身就是全局的，而我們?yōu)槭裁催€非得堅(jiān)守傅里葉變換這個(gè)全局框架下呢？這種變化快的信號(hào)本身就是短時(shí)局部的信號(hào)不斷變換，為什么不是一個(gè)個(gè)局部的去分析信號(hào)呢，這也是小波變換的基本思想，小波將時(shí)頻分析推向了研究高潮，這一時(shí)期小波理論不斷發(fā)展，出現(xiàn)了許多小波，db系列小波，morlet小波，coif系列小波等等，可以看出來，小波研究都是基于小波基函數(shù)的研究，它不像傅里葉變換一樣，基是固定的，而小波基函數(shù)有很多形式。這種靈活性給了小波廣泛的運(yùn)用優(yōu)勢(shì)，但是，科學(xué)就是這樣，一種方法的誕生必然伴隨著缺陷的誕生，小波基的這種靈活性卻給分析人員帶來困擾，如何針對(duì)不同信號(hào)選擇這些基呢？是否有一個(gè)通用基來處理信號(hào)分解任務(wù)呢？還有一個(gè)關(guān)鍵問題是，小波變換受到測(cè)不準(zhǔn)原理的限制，不能無限制的細(xì)分時(shí)間和頻率。當(dāng)然研究人員根據(jù)實(shí)際處理任務(wù)提出了相對(duì)通用的小波基，如圖像處理中的曲波變換和脊波變換，這一時(shí)期小波研究含有另一個(gè)方向，不是一個(gè)小波不能完全符合處理任務(wù)嗎？不是不好選擇小波基嗎？研究人員提出把多個(gè)小波基一起用，必然可以得到更好的結(jié)果，這典型的代表是多小波的發(fā)展。當(dāng)然普通小波變換存在時(shí)移缺陷和頻率混疊的缺陷，而雙樹復(fù)數(shù)小波DTCWT得提出基本解決了該問題。也有針對(duì)小波基難以構(gòu)造的缺陷，提出了提升框架的小波變換，但是小波變換就是小波變換，本質(zhì)還是沒有改變，小波具有的問題，后面發(fā)展出來的方法依然存在這些問題，小波處理方法沒有過時(shí)，但是其諸多問題卻在工程運(yùn)用中越來越明顯，工程運(yùn)用中信號(hào)千變?nèi)f化，小波基選擇始終是擺在分析人員面前的問題，還有就是樓主提到的消噪，小波消噪軟閾值和硬閾值的選擇會(huì)對(duì)信號(hào)消噪能力有巨大影響，選擇不好要么消除不夠，有么太多。更直白的說，你說你消除了噪聲，那么萬一別人問你，你憑什么說你消除的是噪聲而不是有用信號(hào)了，這個(gè)問題是一個(gè)無法回答的問題。當(dāng)然后小波時(shí)代崛起了一批數(shù)據(jù)驅(qū)動(dòng)的非線性非平穩(wěn)信號(hào)處理方法，這些方法的提出就是擺脫小波基選擇困難，測(cè)不準(zhǔn)原理限制的問題，這類方法典型代表就是EMD，EMD方法根據(jù)信號(hào)自身特點(diǎn)不斷分離開來，噪聲和有用信號(hào)安照不同頻帶劃分分離開，當(dāng)然由于EMD方法本身諸多問題，后來發(fā)展出來了EEMD，CEEMD，CEEMDAN。CWEEMDAN等方法，但是其數(shù)學(xué)基礎(chǔ)問題依然存疑，可是，EMD的提出給了人們一個(gè)啟發(fā)，信號(hào)分解就是要讓信號(hào)根據(jù)自己特點(diǎn)分開而不是人工的參與，一切應(yīng)該是自驅(qū)動(dòng)自適應(yīng)才是最好的，受該思想的啟發(fā)，小波研究者也希望將小波推廣到這類方法中，其典型代表SWT和EWT，雖然方法很妙，但是其本質(zhì)問題卻并非數(shù)據(jù)驅(qū)動(dòng)，有興趣可以自行查閱文獻(xiàn)，當(dāng)然繼承該思想的還有VMD、NSP、ITD方法，不得不提一下，這些方法都是非常適合非平穩(wěn)非線性信號(hào)處理，但是分解中的一些參數(shù)設(shè)置卻十分困難。關(guān)于信號(hào)分析最前沿的方法和理論，推薦專欄 信號(hào)處理與機(jī)器學(xué)習(xí) - 知乎專欄 里面講了很多信號(hào)處理方法的核心思想和最新方法。。。。

VMD（變分模態(tài)分解）（https://zhuanlan.zhihu.com/p/66898788?utm_source=qq）

這一篇寫一下變分模態(tài)分解（原始論文：Variational Mode Decomposition），跟原始論文思路思路一致但有一點(diǎn)點(diǎn)不太一樣，原始論文寫的很好，但我不是通信專業(yè)沒有學(xué)過信號(hào)相關(guān)課程一開始看起來有點(diǎn)費(fèi)勁。模態(tài)分解認(rèn)為信號(hào)是由不同“模態(tài)”的子信號(hào)疊加而成的，而變分模態(tài)分解則認(rèn)為信號(hào)是由不同頻率占優(yōu)的子信號(hào)疊加而成的，其目的是要把信號(hào)分解成不同頻率的子信號(hào)。變分模態(tài)分解的分解結(jié)果如圖所示

基礎(chǔ)

一開始論文看不懂的原因是缺少相關(guān)前置知識(shí)，但一旦順下來就會(huì)感覺其實(shí)沒有那么難，難的是作者的思路很巧妙，先寫下我遇到的這些知識(shí)盲點(diǎn)

第一點(diǎn)是傅立葉變換的微分性質(zhì)，

的傅立葉變換為 ，其導(dǎo)數(shù) 的傅立葉變換為

另一點(diǎn)是解析信號(hào)，現(xiàn)實(shí)世界只能采集實(shí)信號(hào)，但實(shí)信號(hào)有很多不好用的性質(zhì)，如存在負(fù)頻率，無法直接得到調(diào)制頻率后的實(shí)信號(hào)等。 設(shè)原始信號(hào)是一個(gè)實(shí)信號(hào)

為了方便表示 為隨 變化的函數(shù)，相當(dāng)于瞬時(shí)頻率，解析信號(hào)是一個(gè)復(fù)信號(hào)，可以通過希爾伯特變換得到 解析信號(hào)的實(shí)部是原本的實(shí)數(shù)信號(hào)，并且經(jīng)過調(diào)頻之后復(fù)數(shù)信號(hào)的實(shí)部仍然是調(diào)頻之后的實(shí)信號(hào)，如對(duì)信號(hào) 增加頻率 ，只需要乘以 即可 由此，其實(shí)部相當(dāng)于在原本頻率 的基礎(chǔ)上增加了頻率 ，如下面的matlab腳本clear;close all;clc; t = 1:0.01:10; %% f1 = sin(20*t).*(t-5).^2; subplot(3,1,1); plot(f1); ylim([-25 25]); %% f2 = sin(50*t).*(t-5).^2; subplot(3,1,2); plot(f2); ylim([-25 25]); %% H = hilbert(f1); f_hat = H.*exp(1i*30.*t); subplot(3,1,3); plot(real(f_hat)); ylim([-25 25]);

matlab的hilbert函數(shù)包括希爾伯特變換和解析函數(shù)轉(zhuǎn)換兩部分，直接得到實(shí)信號(hào)的解析信號(hào)，其中希爾伯特變換

正文

接下來我們看看如何一步一步得到變分模態(tài)分解的思路

原論文通過一個(gè)信號(hào)降噪問題進(jìn)行說明，現(xiàn)需要對(duì)采樣信號(hào)

進(jìn)行降噪重構(gòu)，假設(shè)觀測(cè)信號(hào)是由原始信號(hào)疊加一個(gè)獨(dú)立的高斯噪音 ，需要求 ，又說該等式是一個(gè)不適定問題（ill-posed problem），不滿足識(shí)定問題的三個(gè)條件，所以要用一個(gè)正則化的方法 第一部分是對(duì)原信號(hào)進(jìn)行重構(gòu)，第二部分是為了解決不適定問題的解不唯一，而且不同于機(jī)器學(xué)習(xí)的建模問題一樣 是一個(gè)權(quán)值形式可以直接加權(quán)值 的L1-norm或L2-norm，這里的 是一個(gè)純函數(shù)的形式，其導(dǎo)數(shù)的L2-norm最小化感覺上應(yīng)該是保證了函數(shù) 不會(huì)產(chǎn)生太大的波動(dòng)，這里可不可以跟防止過擬合聯(lián)系起來，接下來說明這個(gè)約束會(huì)使 在頻率域產(chǎn)生什么影響。

下面解這個(gè)式子，首先來看一下直接最小化泛函

能不能解出來 有 ，然后根據(jù)E-L方程的引理 得 ，往下忘了怎么求了。。。。。想起來怎么求再繼續(xù)寫-_-!，這個(gè)直接求的方法與原文沒什么關(guān)系的

由上面的傅立葉變換的微分性質(zhì)知道，用傅立葉變換在復(fù)數(shù)域很方便可以把微分約掉，而且還有一個(gè)叫什么Plancherel傅立葉等距映射的東西，時(shí)域的L2范數(shù)與傅立葉變換到的頻率域的L2范數(shù)等距，故上面的最小化的項(xiàng)可以直接用傅立葉變換轉(zhuǎn)換到頻域

，這里需要注意的是利用上面的那個(gè)什么等距定理有 和 的 是同一個(gè) ，這一點(diǎn)很重要，第二個(gè)是L2范數(shù)大于0，則把其展開為泛函 求最極值，由于這里面只有 直接求偏導(dǎo)即可也不用解微分方程

可以看到得到的

相當(dāng)于對(duì)觀測(cè)信號(hào)在 在頻率段進(jìn)行濾波，過濾掉了高頻部分，這說明加了該導(dǎo)數(shù)的L1正則約束與上面的直觀感覺是一樣的，過濾掉了高頻部分，減弱 的波動(dòng)

再往下進(jìn)行，模態(tài)分解需要把原始信號(hào)分解成多個(gè)子信號(hào)的和，我們?yōu)榱撕驮膶?duì)應(yīng)修改一下符號(hào)表示，

表示觀測(cè)的采樣信號(hào)， 表示分解得到的基函數(shù)，則上面的約束對(duì)象變?yōu)? 同樣先轉(zhuǎn)化為頻率域再求極值 泛函

每個(gè)基函數(shù)基于其他的基函數(shù)更新，相當(dāng)于每個(gè)基函數(shù)是原信號(hào)剩余部分的低通濾波，每次迭代都是保留剩余信號(hào)的低頻率部分。

到現(xiàn)在為止我們發(fā)現(xiàn)每個(gè)基函數(shù)都會(huì)趨向于每次的剩余信號(hào)分量的低頻部分，這與我們?cè)嫉募僭O(shè)“每個(gè)基函數(shù)都有不同的頻率分量”是相悖的，但根據(jù)上面的低通濾波的性質(zhì)，每個(gè)基函數(shù)進(jìn)行特定頻率的濾波應(yīng)該就能解決這個(gè)問題了，那么上面的式子就簡單的變?yōu)?

其中， 為每個(gè)基函數(shù) 的中心頻率，該式就是變分模態(tài)分解的基函數(shù)的更新公式，我們來看一下這個(gè)式子應(yīng)該如何得到，以便于找到中心頻率的更新公式

由上面的一步一步的演化發(fā)現(xiàn)，基函數(shù)對(duì)剩余信號(hào)的低通濾波是由導(dǎo)數(shù)的L2正則最小化帶來的，要得到基函數(shù)的中心頻率約束也要從這個(gè)地方入手，由上面的推導(dǎo)可知每個(gè)基函數(shù)都會(huì)被約束到

頻率附近，那么我們把基函數(shù)的頻率增加各自的中心頻率 得到 ，并保證 即可，則相當(dāng)于對(duì)每個(gè)基函數(shù)乘以了一個(gè) ，這里需要對(duì)頻率進(jìn)行變換，我們沿用文章開頭的解析信號(hào)的性質(zhì)，認(rèn)為 是一個(gè)復(fù)數(shù)的解析信號(hào)，同樣也需要把觀測(cè)信號(hào) 預(yù)先轉(zhuǎn)化為解析信號(hào)，下文默認(rèn)都是解析信號(hào)，則每一個(gè)基函數(shù)轉(zhuǎn)換頻率后的導(dǎo)數(shù)變?yōu)?p>

傅立葉變換得 先來看看傅立葉變換為什么會(huì)得到這個(gè)，而不是 反傅立葉變換 感覺應(yīng)該不是這樣證的，我數(shù)學(xué)不是很好，不怎么會(huì)這個(gè)變量變換

至此，約束變?yōu)?

傅立葉變換 根據(jù)論文原文把 然后求最小就可以得到上面 的更新公式

最后，為了保證每個(gè)點(diǎn)處的重構(gòu)信號(hào)與原信號(hào)盡可能相似，增加了每個(gè)點(diǎn)處的重構(gòu)約束，其實(shí)這一項(xiàng)并不是必需的，最終的約束對(duì)象為

，然后拉格朗日乘子法帶進(jìn)去，但其需要滿足下式才有意義 這個(gè)式子還是符合Parseval定理，故整理一下

最后整理一下更新公式：

其中 使用梯度下降更新

代碼

%matplotlib inline from matplotlib import pyplot as plt import numpy as np from scipy.signal import hilbertT = 1000 fs = 1./T t = np.linspace(0, 1, 1000,endpoint=True) f_1 = 10 f_2 = 50 f_3 = 100 mode_1 = (2 * t) ** 2 mode_2 = np.sin(2 * np.pi * f_1 * t) mode_3 = np.sin(2 * np.pi * f_2 * t) mode_4 = np.sin(2 * np.pi * f_3 * t) f = mode_1 + mode_2 + mode_3 + mode_4 + 0.5 * np.random.randn(1000)plt.figure(figsize=(6,3), dpi=150) plt.plot(f, linewidth=1)

class VMD:def __init__(self, K, alpha, tau, tol=1e-7, maxIters=200, eps=1e-9):""":param K: 模態(tài)數(shù):param alpha: 每個(gè)模態(tài)初始中心約束強(qiáng)度:param tau: 對(duì)偶項(xiàng)的梯度下降學(xué)習(xí)率:param tol: 終止閾值:param maxIters: 最大迭代次數(shù):param eps: eps"""self.K =Kself.alpha = alphaself.tau = tauself.tol = tolself.maxIters = maxItersself.eps = epsdef __call__(self, f):T = f.shape[0]t = np.linspace(1, T, T) / Tomega = t - 1. / T# 轉(zhuǎn)換為解析信號(hào)f = hilbert(f)f_hat = np.fft.fft(f)u_hat = np.zeros((self.K, T), dtype=np.complex)omega_K = np.zeros((self.K,))lambda_hat = np.zeros((T,), dtype=np.complex)# 用以判斷u_hat_pre = np.zeros((self.K, T), dtype=np.complex)u_D = self.tol + self.eps# 迭代n = 0while n < self.maxIters and u_D > self.tol:for k in range(self.K):# u_hatsum_u_hat = np.sum(u_hat, axis=0) - u_hat[k, :]res = f_hat - sum_u_hatu_hat[k, :] = (res + lambda_hat / 2) / (1 + self.alpha * (omega - omega_K[k]) ** 2)# omegau_hat_k_2 = np.abs(u_hat[k, :]) ** 2omega_K[k] = np.sum(omega * u_hat_k_2) / np.sum(u_hat_k_2)# lambda_hatsum_u_hat = np.sum(u_hat, axis=0)res = f_hat - sum_u_hatlambda_hat -= self.tau * resn += 1u_D = np.sum(np.abs(u_hat - u_hat_pre) ** 2)u_hat_pre[::] = u_hat[::]# 重構(gòu)，反傅立葉之后取實(shí)部u = np.real(np.fft.ifft(u_hat, axis=-1))omega_K = omega_K * Tidx = np.argsort(omega_K)omega_K = omega_K[idx]u = u[idx, :]return u, omega_KK = 4 alpha = 2000 tau = 1e-6 vmd = VMD(K, alpha, tau) u, omega_K = vmd(f) omega_K # array([0.85049797, 10.08516203, 50.0835613, 100.13259275])) plt.figure(figsize=(5,7), dpi=200) plt.subplot(4,1,1) plt.plot(mode_1, linewidth=0.5, linestyle='--') plt.plot(u[0,:], linewidth=0.2, c='r')plt.subplot(4,1,2) plt.plot(mode_2, linewidth=0.5, linestyle='--') plt.plot(u[1,:], linewidth=0.2, c='r')plt.subplot(4,1,3) plt.plot(mode_3, linewidth=0.5, linestyle='--') plt.plot(u[2,:], linewidth=0.2, c='r')plt.subplot(4,1,4) plt.plot(mode_4, linewidth=0.5, linestyle='--') plt.plot(u[3,:], linewidth=0.2, c='r') # [<matplotlib.lines.Line2D at 0x7fac532f4dd8>]

---

可以看到有比較強(qiáng)的端點(diǎn)效應(yīng)，端點(diǎn)處會(huì)有重疊，文章原始論文中采用的方法是對(duì)稱拼接的方法，把原信號(hào)復(fù)制一份然后拼成兩半，一半對(duì)稱放前面，一般對(duì)稱放后面

%matplotlib inline from matplotlib import pyplot as plt import numpy as np from scipy.signal import hilbertT = 1000 fs = 1./T t = np.linspace(0, 1, 1000,endpoint=True)f_1 = 10 f_2 = 50 f_3 = 100 mode_1 = np.sin(2 * np.pi * f_1 * t) mode_2 = np.sin(2 * np.pi * f_2 * t) mode_3 = np.sin(2 * np.pi * f_3 * t) f = np.concatenate((mode_1[:301], mode_2[301:701], mode_3[701:])) + 0.1 * np.random.randn(1000)plt.figure(figsize=(6,3), dpi=150) plt.plot(f, linewidth=0.5) # [<matplotlib.lines.Line2D at 0x7fc2134b8630>]

class VMD:def __init__(self, K, alpha, tau, tol=1e-7, maxIters=200, eps=1e-9):""":param K: 模態(tài)數(shù):param alpha: 每個(gè)模態(tài)初始中心約束強(qiáng)度:param tau: 對(duì)偶項(xiàng)的梯度下降學(xué)習(xí)率:param tol: 終止閾值:param maxIters: 最大迭代次數(shù):param eps: eps"""self.K =Kself.alpha = alphaself.tau = tauself.tol = tolself.maxIters = maxItersself.eps = epsdef __call__(self, f):N = f.shape[0]# 對(duì)稱拼接f = np.concatenate((f[:N//2][::-1], f, f[N//2:][::-1]))T = f.shape[0]t = np.linspace(1, T, T) / Tomega = t - 1. / T# 轉(zhuǎn)換為解析信號(hào)f = hilbert(f)f_hat = np.fft.fft(f)u_hat = np.zeros((self.K, T), dtype=np.complex)omega_K = np.zeros((self.K,))lambda_hat = np.zeros((T,), dtype=np.complex)# 用以判斷u_hat_pre = np.zeros((self.K, T), dtype=np.complex)u_D = self.tol + self.eps# 迭代n = 0while n < self.maxIters and u_D > self.tol:for k in range(self.K):# u_hatsum_u_hat = np.sum(u_hat, axis=0) - u_hat[k, :]res = f_hat - sum_u_hatu_hat[k, :] = (res + lambda_hat / 2) / (1 + self.alpha * (omega - omega_K[k]) ** 2)# omegau_hat_k_2 = np.abs(u_hat[k, :]) ** 2omega_K[k] = np.sum(omega * u_hat_k_2) / np.sum(u_hat_k_2)# lambda_hatsum_u_hat = np.sum(u_hat, axis=0)res = f_hat - sum_u_hatlambda_hat -= self.tau * resn += 1u_D = np.sum(np.abs(u_hat - u_hat_pre) ** 2)u_hat_pre[::] = u_hat[::]# 重構(gòu)，反傅立葉之后取實(shí)部u = np.real(np.fft.ifft(u_hat, axis=-1))u = u[:, N//2 : N//2 + N]omega_K = omega_K * T / 2idx = np.argsort(omega_K)omega_K = omega_K[idx]u = u[idx, :]return u, omega_KK = 3 alpha = 2000 tau = 1e-6 vmd = VMD(K, alpha, tau)u, omega_K = vmd(f) omega_K # array([ 9.68579292, 50.05232833, 100.12321047])plt.figure(figsize=(5,7), dpi=200) plt.subplot(3,1,1) plt.plot(mode_1, linewidth=0.5, linestyle='--') plt.plot(u[0,:], linewidth=0.2, c='r')plt.subplot(3,1,2) plt.plot(mode_2, linewidth=0.5, linestyle='--') plt.plot(u[1,:], linewidth=0.2, c='r')plt.subplot(3,1,3) plt.plot(mode_3, linewidth=0.5, linestyle='--') plt.plot(u[2,:], linewidth=0.2, c='r') # [<matplotlib.lines.Line2D at 0x7fc2134075c0>]

好像結(jié)果要好一點(diǎn)，VMD的一個(gè)缺點(diǎn)是K的值對(duì)結(jié)果有很大影響，但這個(gè)迭代過程，怎么開始怎么迭代怎么結(jié)束都可以自己控制，感覺可以按照自己的需求來定制怎么動(dòng)態(tài)決定K

總結(jié)

以上是生活随笔為你收集整理的基于小波变换的信号降噪处理及仿真研究_信号处理方法推荐--1（转载自用，侵删）...的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。