快捷導(dǎo)航

一維信號(hào)小波去噪原理解析及python實(shí)現(xiàn)方式

更新時(shí)間：2023年06月25日 09:58:25 作者：衷科知眠

這篇文章主要介紹了一維信號(hào)小波去噪原理解析及python實(shí)現(xiàn)方式，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。如有錯(cuò)誤或未考慮完全的地方，望不吝賜教

一維信號(hào)小波去噪原理及python實(shí)現(xiàn)

信號(hào)去噪是經(jīng)常用到的信號(hào)預(yù)處理過(guò)程，以達(dá)到在保留原有信號(hào)真是信息的基礎(chǔ)上盡可能低降低或者消除噪聲，獲得更高質(zhì)量的信號(hào)，從而為下一步的處理奠定基礎(chǔ)。

去噪方法可分為時(shí)域方法與頻域方法。時(shí)域方法是指直接在原始信號(hào)上進(jìn)行處理，比如均值濾波器、中值濾波器、EMD分解等方法。頻域方法是指在信號(hào)的變換域進(jìn)行去噪然后再恢復(fù)到時(shí)域得到去噪后的信號(hào)，比如小波變換、傅里葉變換等方法。

無(wú)論是一維信號(hào)還是二維信號(hào)其原理都是一樣的，只不過(guò)實(shí)現(xiàn)上所有不同。這里主要記錄一下我目前需要處理的一維信號(hào)的去噪。

噪聲可分為加性噪聲和乘性噪聲兩種形式，加性噪聲一般假設(shè)其與原始信號(hào)無(wú)關(guān)，假設(shè)含噪信號(hào)為 y(i) ，其可表示為：

y(i) = x(i)+n(i)

其中 x(i) 表示原始的干凈信號(hào)， n(i) 表示噪聲信號(hào)。

而乘性噪聲則假設(shè)噪聲與信號(hào)有關(guān)，假如噪聲和信號(hào)成正比，則含噪信號(hào) y(i) 可表示為：

y(i)=x(i)[1+n(i)]

其中第二個(gè)噪聲項(xiàng)是受原始信號(hào)影響的，原始信號(hào)的值越大，則噪聲越大。

加性噪聲的處理與計(jì)算都比較簡(jiǎn)單，所以一般均假設(shè)噪聲為加性噪聲。

基于小波變換的去噪方法其具體原理這里就不闡述了，其主要思想就是噪聲在經(jīng)過(guò)小波變換后的小波系數(shù)會(huì)遠(yuǎn)遠(yuǎn)小于原始信號(hào)的小波系數(shù)，因此可通過(guò)簡(jiǎn)單的閾值處理就能夠去除絕大部分的噪聲，同時(shí)又不會(huì)對(duì)原始信號(hào)造成很大的損失。對(duì)具體原理感興趣的同學(xué)可以參考其他相關(guān)資料。

python 實(shí)現(xiàn)

使用python去噪，不需要我們自己動(dòng)手實(shí)現(xiàn)相應(yīng)的原理，安裝相應(yīng)的包即可完成。

首先安裝小波變換的python包PyWavelets

pip install PyWavelets

安裝好之后就可以直接使用，我們可通過(guò)定義去噪函數(shù)將去噪過(guò)程打包，假設(shè)去噪函數(shù)名為wavelet_denoising, 其去噪過(guò)程及參數(shù)設(shè)置代碼如下：

def wavelet_denoising(signal):
    # signal: list
    wavelet_basis = 'db8'  # 此處選擇db8小波基
    w = pywt.Wavelet(wavelet_basis) 
    maxlevel = pywt.dwt_max_level(len(signal), w.dec_len) # 根據(jù)數(shù)據(jù)長(zhǎng)度計(jì)算分解層數(shù)
    threshold = 0.1  # threshold for filtering
    coeffs = pywt.wavedec(signal, wavelet_basis, level=maxlevel)  # 將信號(hào)進(jìn)行小波分解
    # 使用閾值濾除噪聲
    for i in range(1, len(coeffs)):
        coeffs[i] = pywt.threshold(coeffs[i], threshold*max(coeffs[i]))  # 將噪聲濾波
        # pywt.threshold 方法有三種閾值模式：硬閾值、軟閾值、軟硬結(jié)合
        # 默認(rèn)為soft 軟閾值模式，其調(diào)用格式為
        # pywt.threshold(signal, threshold_value, mode='soft', substitute=0)
        # 處理后值為 signal/np.abs(signal) * np.maximum(np.abs(signal) - threshold_value, 0)
        # 具體可參考pywavelets官方文檔的說(shuō)明
    # 重建信號(hào)
    signalrec= pywt.waverec(coeffs, wavelet_basis)
    return signalrec

通過(guò)調(diào)用wavelet_denoising方法，即可實(shí)現(xiàn)對(duì)原始信號(hào)的去噪。

值得注意的是去噪前后信號(hào)的長(zhǎng)度會(huì)有所變化，若原始信號(hào)樣本長(zhǎng)度為奇數(shù)，則去噪后長(zhǎng)度加1，若原始長(zhǎng)度為偶數(shù)，則去噪后長(zhǎng)度不變。

通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)，如果要計(jì)算去噪前后信號(hào)的差，對(duì)于原始數(shù)據(jù)長(zhǎng)度為奇數(shù)時(shí)，應(yīng)該計(jì)算去噪后信號(hào)的第一個(gè)元素到倒數(shù)第二個(gè)元素這個(gè)子序列與原始信號(hào)之間的差值，也即舍去因去噪而多出的最后一個(gè)值。

如果使用去噪后的第二個(gè)元素到最后一個(gè)元素與原始信號(hào)的差，則可能出現(xiàn)偏差較大的情況。

如下圖所示，其中藍(lán)色線為去噪信號(hào)的第一個(gè)值到倒數(shù)第二個(gè)值與原始信號(hào)的差，紅色線為去噪信號(hào)的第二個(gè)值到最后一個(gè)值與原始信號(hào)的差值。

小波去噪代碼

先放代碼，代碼如下：

a = np.loadtxt('C:\\Users\\Administrator\\Desktop\\數(shù)據(jù)文件\\20201008009.txt', dtype=float)# 20201008009為一個(gè)一維列表
index = []
data = []
coffs = []
for i in range(len(a) - 1):
    X = float(i)
    Y = float(a[i])
    index.append(X)
    data.append(Y)
# create wavelet object and define parameters
w = pywt.Wavelet('db8')  # 選用Daubechies8小波
maxlev = pywt.dwt_max_level(len(data), w.dec_len)
print("maximum level is" + str(maxlev))
threshold = 0.1  # Threshold for filtering
# Decompose into wavelet components,to the level selected:
coffs = pywt.wavedec(data, 'db8', level=maxlev)  # 將信號(hào)進(jìn)行小波分解
for i in range(1, len(coffs)):
    coffs[i] = pywt.threshold(coffs[i], threshold * max(coffs[i]))
datarec = pywt.waverec(coffs, 'db8')  # 將信號(hào)進(jìn)行小波重構(gòu)
mintime = 0
maxtime = mintime + len(data)
print(mintime, maxtime)
plt.rcParams['font.sans-serif']=['SimHei'] #使得標(biāo)題輸出為中文
plt.rcParams['axes.unicode_minus'] = False #使得標(biāo)題輸出為中文
plt.figure(1)
#plt.subplot(3, 1, 1)
plt.plot(index[mintime:maxtime], data[mintime:maxtime])
plt.title("原始信號(hào)")
plt.figure(2)
#plt.subplot(3, 1, 2)
plt.plot(index[mintime:maxtime], datarec[mintime:maxtime])
plt.title("小波降噪信號(hào)")
plt.figure(3)
#plt.subplot(3, 1, 3)
plt.plot(index[mintime:maxtime], data[mintime:maxtime] - datarec[mintime:maxtime])
plt.xlabel('time (s)')
plt.ylabel('error (uV)')
plt.tight_layout()
plt.show()