小波降噪与重构例子 python
原理講解
傅里葉變換
關(guān)于傅里葉變換的基本概念在此我們就不再贅述了,
下面我們主要將傅里葉變換的不足。即我們知道傅里葉變化可以分析信號的頻譜,那么為什么還要提出小波變換?答案“對非平穩(wěn)過程,傅里葉變換有局限性”。看如下一個(gè)簡單的信號:
做完FFT(傅里葉變換)后,可以在頻譜上看到清晰的四條線,信號包含四個(gè)頻率部分。一切沒有問題,但是,如果是非平穩(wěn)信號呢?
如上圖,最上邊的是頻率始終不變的平穩(wěn)信號。而下邊兩個(gè)則是頻率隨著時(shí)間改變的非平穩(wěn)信號,它們同樣包含和最上信號相同頻率的四個(gè)成分。做FFT后,我們發(fā)現(xiàn)這三個(gè)時(shí)域上有巨大差異的信號,頻譜(幅值譜)卻非常一致。尤其是下邊兩個(gè)非平穩(wěn)信號,我們從頻譜上無法區(qū)分它們,因?yàn)樗鼈儼乃膫€(gè)頻率的信號的成分確實(shí)是一樣的,只是出現(xiàn)的先后順序不同。
可見,傅里葉變換處理非平穩(wěn)信號有先天缺陷。它只能獲取一段信號總體上包含哪些頻率的成分,但是對各成分出現(xiàn)的時(shí)刻并無所知。因此,時(shí)域相差很大的兩個(gè)信號,可能頻譜一樣。
然而平穩(wěn)信號大多是制造出來的,自然界的大量信號幾乎都是非平穩(wěn)的。
短時(shí)傅里葉變換
一個(gè)簡單可行的方法就是——加窗。 “把整個(gè)時(shí)域過程分解成無數(shù)個(gè)等長的小過程,每個(gè)小過程近似平穩(wěn),再傅里葉變換,就知道在哪個(gè)時(shí)間點(diǎn)上出現(xiàn)了什么頻率了。”這就是短時(shí)傅里葉變換。
看圖:
時(shí)域上分成一段一段做FFT,不就知道頻率成分隨著時(shí)間的變化情況了嗎?
使用STFT存在一個(gè)問題,我們應(yīng)該用多寬的窗函數(shù)?
窗太寬太窄都有問題:
小波變換
那么你可能會想到,讓窗口大小變起來,多做幾次STFT不就可以了嗎?!沒錯(cuò),小波變換就有著這樣的思路。
但事實(shí)上小波并不是這么做的(有人認(rèn)為“小波變換就是根據(jù)算法,加不等長的窗,對每一小部分進(jìn)行傅里葉變換”,這是不準(zhǔn)確的。小波變換并沒有采用窗的思想,更沒有做傅里葉變換。)
至于為什么不采用可變窗的STFT呢,我認(rèn)為是因?yàn)檫@樣做冗余會太嚴(yán)重,STFT做不到正交化,這也是它的一大缺陷。
于是小波變換的出發(fā)點(diǎn)和STFT還是不同的。STFT是給信號加窗,分段做FFT;而小波直接把傅里葉變換的基給換了——將無限長的三角函數(shù)基換成了有限長的會衰減的小波基。這樣不僅能夠獲取頻率,還可以定位到時(shí)間了~
解釋:
來我們回顧一下傅里葉變換吧,沒弄清傅里葉變換為什么能得到信號各個(gè)頻率成分的同學(xué)也可以再借我的圖理解下。
傅里葉變換把無限長的三角函數(shù)作為基函數(shù):
這個(gè)基函數(shù)會伸縮,平移(其實(shí)是兩個(gè)正交基的分解)。縮得窄,對應(yīng)高頻;伸得寬,對應(yīng)低頻。然后這個(gè)基函數(shù)不斷和信號做相乘。某一個(gè)尺度(寬窄)下乘出來的結(jié)果,就可以理解成信號所包含的當(dāng)前尺度對應(yīng)頻率成分有多少。于是,基函數(shù)會在某些尺度下,與信號相乘得到一個(gè)很大的值,因?yàn)榇藭r(shí)二者有一種重合關(guān)系。那么我們就知道信號包含該頻率的成分的多少。
看,這兩種尺度能乘出一個(gè)大的值(相關(guān)度高),所以信號包含較多的這兩個(gè)頻率成分,在頻譜上這兩個(gè)頻率會出現(xiàn)兩個(gè)峰。
以上,就是粗淺意義上傅里葉變換的原理。
如前邊所說,小波做的改變就在于,將無限長的三角函數(shù)基換成了有限長的會衰減的小波基。
從公式可以看出,不同于傅里葉變換,變量只有頻率ω,小波變換有兩個(gè)變量:尺度a(scale)和平移量 τ(translation)。尺度a控制小波函數(shù)的伸縮,平移量 τ控制小波函數(shù)的平移。尺度就對應(yīng)于頻率(反比),平移量 τ就對應(yīng)于時(shí)間。
當(dāng)伸縮、平移到這么一種重合情況時(shí),也會相乘得到一個(gè)大的值。這時(shí)候和傅里葉變換不同的是,這不僅可以知道信號有這樣頻率的成分,而且知道它在時(shí)域上存在的具體位置。
而當(dāng)我們在每個(gè)尺度下都平移著和信號乘過一遍后,我們就知道信號在每個(gè)位置都包含哪些頻率成分。
原理部分參考文獻(xiàn)
原理部分參考來源
小波變換原理
感謝大佬,感謝大佬。
代碼部分
簡單例子,python版
讀取數(shù)據(jù)
#小波降噪與重構(gòu)畫圖
#獲取重構(gòu)后的趨勢部分和噪聲
def choose(data, w):data=data[:,-1]#只對電價(jià)數(shù)據(jù)進(jìn)行去噪重構(gòu)mode = pywt.Modes.smoothw = pywt.Wavelet(w) # 選取小波函數(shù)a = dataca = [] # 近似分量cd = [] # 細(xì)節(jié)分量for i in range(1):#1階變換(a, d) = pywt.dwt(a, w, mode) # 進(jìn)行1階離散小波變換 dwt為分解ca.append(a)cd.append(d)rec_a = []rec_d = []for i, coeff in enumerate(ca):coeff_list = [coeff, None] +[None] * i#填充長度rec_a.append(pywt.waverec(coeff_list, w)) # waverec重構(gòu)rec_a=np.array(rec_a[-1])#-1取最后一個(gè)基波rec_a=rec_a.reshape(-1,1)for i, coeff in enumerate(cd):coeff_list = [None, coeff] + [None] * irec_d.append(pywt.waverec(coeff_list, w))#取出所有噪聲波,1階變換,有3個(gè)噪聲rec_d=np.array(rec_d)return rec_a,rec_d#取出趨勢,噪聲 trainingrec_a, trainingrec_d=choose(txt, 'sym5')#訓(xùn)練集電價(jià)趨勢 和噪聲trainingrec_a為趨勢
trainingrec_d為噪聲
算例2 多階變換 代碼
import numpy as np import matplotlib.pyplot as plt import matplotlib; matplotlib.use('TkAgg') from pylab import * mpl.rcParams['font.sans-serif'] = ['SimHei'] mpl.rcParams['axes.unicode_minus'] = False import pywt #導(dǎo)入小波包#==============生成數(shù)據(jù)============== X=np.random.uniform(90,100,100) #生成范圍90-100的數(shù),數(shù)據(jù)100個(gè)#============小波繪圖============== def plot_signal_decomp(data, w,title):""":param data: 要進(jìn)行小波分解重構(gòu)的數(shù)據(jù):param w: 小波基名稱:param title: 繪圖名稱:return: 圖"""mode = pywt.Modes.smoothw = pywt.Wavelet(w) # 選取小波函數(shù)a = dataca = [] # 近似分量 ,趨勢cd = [] # 細(xì)節(jié)分量,噪聲n=3 #定義階數(shù),這里定義3階for i in range(n):#3階(a, d) = pywt.dwt(a, w, mode) # 進(jìn)行n階離散小波變換ca.append(a)#存放趨勢cd.append(d)#存儲噪聲rec_a = []#存放重構(gòu)后的趨勢rec_d = []#存放重構(gòu)后的噪聲for i, coeff in enumerate(ca):coeff_list = [coeff, None] + [None] * irec_a.append(pywt.waverec(coeff_list, w)) # 重構(gòu)for i, coeff in enumerate(cd):coeff_list = [None, coeff] + [None] * irec_d.append(pywt.waverec(coeff_list, w))print('趨勢個(gè)數(shù)len(rec_a):',len(rec_a))print('噪聲個(gè)數(shù)len(rec_d):',len(rec_d))#趨勢只要最后一個(gè),噪聲全要。 加原始數(shù)據(jù)個(gè)數(shù)。所以子圖個(gè)數(shù)=n+1+2fig = plt.figure()for i in range(n+2):if i==n+1:plt.subplot(n+2, 1, i+1)plt.plot(data, label='原數(shù)據(jù)曲線')plt.legend()elif i==n:plt.subplot(n+2, 1, i+1)plt.plot(rec_a[-1], 'g', label='數(shù)據(jù)曲線趨勢')plt.legend()else:plt.subplot(n+2, 1, i+1)plt.plot(rec_d[i], 'r', label='數(shù)據(jù)曲線噪聲%d'%(i+1))plt.legend()plt.show()plot_signal_decomp(X, 'sym5','sym5')#這里選擇sym5小波,小波還有#=======獲取重構(gòu)后的數(shù)據(jù)======== def choose(data, w):""":param data: 要進(jìn)行小波分解重構(gòu)的數(shù)據(jù):param w: 小波基:return: 趨勢 和噪聲"""mode = pywt.Modes.smoothw = pywt.Wavelet(w) # 選取小波函數(shù)a = dataca = [] # 近似分量cd = [] # 細(xì)節(jié)分量n = 3 # 定義階數(shù),這里定義3階for i in range(n): # 3階(a, d) = pywt.dwt(a, w, mode) # 進(jìn)行n階離散小波變換ca.append(a) # 存放趨勢cd.append(d) # 存儲噪聲rec_a = [] # 存放重構(gòu)后的趨勢rec_d = [] # 存放重構(gòu)后的噪聲for i, coeff in enumerate(ca):coeff_list = [coeff, None] + [None] * irec_a.append(pywt.waverec(coeff_list, w)) # 重構(gòu)for i, coeff in enumerate(cd):coeff_list = [None, coeff] + [None] * irec_d.append(pywt.waverec(coeff_list, w))#返回最后一個(gè)趨勢、所有噪聲 return rec_a[-1],rec_dtrainingrec_a, trainingrec_d=choose(X, 'sym5')#所有趨勢 和噪聲 print('趨勢\n',trainingrec_a) print('噪聲1\n',trainingrec_d[0]) print('噪聲2\n',trainingrec_d[1]) print('噪聲3\n',trainingrec_d[2])最后得到的是 趨勢和3部分噪聲。。重構(gòu)一般是:舍棄一部分噪聲,如(保留2部分噪聲)。重構(gòu)結(jié)果=噪聲1+噪聲2+趨勢。。或者只要趨勢。。 原式數(shù)據(jù)=噪聲1+噪聲2+噪聲3+趨勢。若不舍棄噪聲,小波分解沒有意義。。
總結(jié)
以上是生活随笔為你收集整理的小波降噪与重构例子 python的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: cambodia怎么读?
- 下一篇: kmean python实现