快捷導(dǎo)航

Python小波變換去噪的原理解析

更新時間：2021年12月31日 16:37:05 作者：sgggr

這篇文章主要介紹了Python小波變換去噪,對于去噪效果好壞的評價，常用信號的信噪比（SNR）與估計信號同原始信號的均方根誤差（RMSE）來判斷，需要的朋友可以參考下

一.小波去噪的原理

信號產(chǎn)生的小波系數(shù)含有信號的重要信息，將信號經(jīng)小波分解后小波系數(shù)較大，噪聲的小波系數(shù)較小，并且噪聲的小波系數(shù)要小于信號的小波系數(shù)，通過選取一個合適的閥值，大于閥值的小波系數(shù)被認(rèn)為是有信號產(chǎn)生的，應(yīng)予以保留，小于閥值的則認(rèn)為是噪聲產(chǎn)生的，置為零從而達到去噪的目的。
小波閥值去噪的基本問題包括三個方面：小波基的選擇，閥值的選擇，閥值函數(shù)的選擇。

(1) 小波基的選擇：通常我們希望所選取的小波滿足以下條件：正交性、高消失矩、緊支性、對稱性或反對稱性。但事實上具有上述性質(zhì)的小波是不可能存在的，因為小波是對稱或反對稱的只有Haar小波，并且高消失矩與緊支性是一對矛盾，所以在應(yīng)用的時候一般選取具有緊支的小波以及根據(jù)信號的特征來選取較為合適的小波。
(2) 閥值的選擇：直接影響去噪效果的一個重要因素就是閥值的選取，不同的閥值選取將有不同的去噪效果。目前主要有通用閥值（VisuShrink）、SureShrink閥值、Minimax閥值、BayesShrink閥值等。
(3) 閥值函數(shù)的選擇：閥值函數(shù)是修正小波系數(shù)的規(guī)則，不同的反之函數(shù)體現(xiàn)了不同的處理小波系數(shù)的策略。最常用的閥值函數(shù)有兩種：一種是硬閥值函數(shù)，另一種是軟閥值函數(shù)。還有一種介于軟、硬閥值函數(shù)之間的Garrote函數(shù)。

另外，對于去噪效果好壞的評價，常用信號的信噪比（SNR）與估計信號同原始信號的均方根誤差（RMSE）來判斷。

二，在python中使用小波分析進行閾值去噪聲,使用pywt.threshold函數(shù)

#coding=gbk
#使用小波分析進行閾值去噪聲,使用pywt.threshold
 
import pywt
import numpy as np
import pandas as pd 
import matplotlib.pyplot as plt 
import math 
 
data = np.linspace(1, 10, 10)
print(data)
# [ 1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]
# pywt.threshold(data, value, mode, substitute) mode 模式有4種，soft, hard, greater, less; substitute是替換值可以點進函數(shù)里看,data/np.abs(data) * np.maximum(np.abs(data) - value, 0)
 
data_soft = pywt.threshold(data=data, value=6, mode='soft', substitute=12)
print(data_soft)
# [12. 12. 12. 12. 12.  0.  1.  2.  3.  4.] 將小于6 的值設(shè)置為12， 大于等于6 的值全部減去6
 
data_hard = pywt.threshold(data=data, value=6, mode='hard', substitute=12)
print(data_hard)
# [12. 12. 12. 12. 12.  6.  7.  8.  9. 10.] 將小于6 的值設(shè)置為12， 其余的值不變
 
data_greater = pywt.threshold(data, 6, 'greater', 12)
print(data_greater)
# [12. 12. 12. 12. 12.  6.  7.  8.  9. 10.] 將小于6 的值設(shè)置為12，大于等于閾值的值不變化
 
data_less = pywt.threshold(data, 6, 'less', 12)
print(data_less)
# [ 1.  2.  3.  4.  5.  6. 12. 12. 12. 12.] 將大于6 的值設(shè)置為12， 小于等于閾值的值不變

三，在python中使用ecg心電信號進行小波去噪實驗

#-*-coding:utf-8-*-

import matplotlib.pyplot as plt
import pywt
import math
import numpy as np

#get Data
ecg=pywt.data.ecg()  #生成心電信號
index=[]
data=[]
coffs=[]

for i in range(len(ecg)-1):
    X=float(i)
    Y=float(ecg[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  #Threshold for filtering

#Decompose into wavelet components,to the level selected:
coffs=pywt.wavedec(data,'db8',level=maxlev) #將信號進行小波分解

for i in range(1,len(coffs)):
    coffs[i]=pywt.threshold(coffs[i],threshold*max(coeffs[i]))

datarec=pywt.waverec(coffs,'db8')#將信號進行小波重構(gòu)

mintime=0
maxtime=mintime+len(data) 
print(mintime,maxtime)

plt.figure()
plt.subplot(3,1,1)
plt.plot(index[mintime:maxtime], data[mintime:maxtime])
plt.xlabel('time (s)')
plt.ylabel('microvolts (uV)')
plt.title("Raw signal")
plt.subplot(3, 1, 2)
plt.plot(index[mintime:maxtime], datarec[mintime:maxtime])
plt.xlabel('time (s)')
plt.ylabel('microvolts (uV)')
plt.title("De-noised signal using wavelet techniques")
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()