溫馨提示×

溫馨提示×

您好,登錄后才能下訂單哦!

密碼登錄×
登錄注冊×
其他方式登錄
點擊 登錄注冊 即表示同意《億速云用戶服務條款》

詳解python實現(xiàn)小波變換的一個簡單例子

發(fā)布時間:2020-08-25 08:47:15 來源:腳本之家 閱讀:289 作者:weixin_33888907 欄目:開發(fā)技術(shù)

最近工作需要,看了一下小波變換方面的東西,用python實現(xiàn)了一個簡單的小波變換類,將來可以用在工作中。

簡單說幾句原理,小波變換類似于傅里葉變換,都是把函數(shù)用一組正交基函數(shù)展開,選取不同的基函數(shù)給出不同的變換。例如傅里葉變換,選擇的是sin和cos,或者exp(ikx)這種復指數(shù)函數(shù);而小波變換,選取基函數(shù)的方式更加靈活,可以根據(jù)要處理的數(shù)據(jù)的特點(比如某一段上信息量比較多),在不同尺度上采用不同的頻寬來對已知信號進行分解,從而盡可能保留多一點信息,同時又避免了原始傅里葉變換的大計算量。以下計算采用的是haar基,它把函數(shù)分為2段(A1和B1,但第一次不分),對第一段內(nèi)相鄰的2個采樣點進行變換(只考慮A1),變換矩陣為

sqrt(0.5)       sqrt(0.5)

sqrt(0.5)        -sqrt(0.5)

變換完之后,再把第一段(A1)分為兩段,同樣對相鄰的點進行變換,直到無法再分。

下面直接上代碼

Wavelet.py

import math
 
class wave:
  def __init__(self):
    M_SQRT1_2 = math.sqrt(0.5)
    self.h2 = [M_SQRT1_2, M_SQRT1_2]
    self.g1 = [M_SQRT1_2, -M_SQRT1_2]
    self.h3 = [M_SQRT1_2, M_SQRT1_2]
    self.g2 = [M_SQRT1_2, -M_SQRT1_2]
    self.nc = 2
    self.offset = 0
 
  def __del__(self):
    return
 
class Wavelet:
  def __init__(self, n):
    self._haar_centered_Init()
    self._scratch = []
    for i in range(0,n):
      self._scratch.append(0.0)
    return
    
  def __del__(self):
    return
    
  def transform_inverse(self, list, stride):
    self._wavelet_transform(list, stride, -1)
    return
    
  def transform_forward(self, list, stride):
    self._wavelet_transform(list, stride, 1)
    return
    
  def _haarInit(self):
    self._wave = wave()
    self._wave.offset = 0
    return
    
  def _haar_centered_Init(self):
    self._wave = wave()
    self._wave.offset = 1
    return
    
  def _wavelet_transform(self, list, stride, dir):
    n = len(list)
    if (len(self._scratch) < n):
      print("not enough workspace provided")
      exit()
    if (not self._ispower2(n)):
      print("the list size is not a power of 2")
      exit()
    
    if (n < 2):
      return
 
    if (dir == 1): # 正變換
      i = n
      while(i >= 2):
        self._step(list, stride, i, dir)
        i = i>>1
       
    if (dir == -1):  # 逆變換
      i = 2
      while(i <= n):
        self._step(list, stride, i, dir)
        i = i << 1
    return
    
  def _ispower2(self, n):
    power = math.log(n,2)
    intpow = int(power)
    intn = math.pow(2,intpow)
    if (abs(n - intn) > 1e-6):
      return False
    else:
      return True
      
  def _step(self, list, stride, n, dir):
    for i in range(0, len(self._scratch)):
      self._scratch[i] = 0.0
    
    nmod = self._wave.nc * n
    nmod -= self._wave.offset
    n1 = n - 1
    nh = n >> 1
    
    if (dir == 1): # 正變換
      ii = 0
      i = 0
      while (i < n):
        h = 0
        g = 0
        ni = i + nmod
        for k in range(0, self._wave.nc):
          jf = n1 & (ni + k)
          h += self._wave.h2[k] * list[stride*jf]
          g += self._wave.g1[k] * list[stride*jf]
        self._scratch[ii] += h
        self._scratch[ii + nh] += g
        i += 2
        ii += 1
    
    if (dir == -1):  # 逆變換
      ii = 0
      i = 0
      while (i < n):
        ai = list[stride*ii]
        ai1 = list[stride*(ii+nh)]
        ni = i + nmod
        for k in range(0, self._wave.nc):
          jf = n1 & (ni + k)
          self._scratch[jf] += self._wave.h3[k] * ai + self._wave.g2[k] * ai1
        i += 2
        ii += 1
        
    for i in range(0, n):
      list[stride*i] = self._scratch[i]

測試代碼如下:

test.py

import math
import Wavelet
 
waveletn = 256
waveletnc = 20  #保留的分量數(shù)
wavelettest = Wavelet.Wavelet(waveletn)
waveletorigindata = []
waveletdata = []
for i in range(0, waveletn):
  waveletorigindata.append(math.sin(i)*math.exp(-math.pow((i-100)/50,2))+1)
  waveletdata.append(waveletorigindata[-1])
  
Wavelet.wavelettest.transform_forward(waveletdata, 1)
newdata = sorted(waveletdata, key = lambda ele: abs(ele), reverse=True)
for i in range(waveletnc, waveletn):  # 篩選出前 waveletnc個分量保留
  for j in range(0, waveletn):
    if (abs(newdata[i] - waveletdata[j]) < 1e-6):
      waveletdata[j] = 0.0
      break
  
Wavelet.wavelettest.transform_inverse(waveletdata, 1)
waveleterr = 0.0
for i in range(0, waveletn):
  print(waveletorigindata[i], ",", waveletdata[i])
  waveleterr += abs(waveletorigindata[i] - waveletdata[i])/abs(waveletorigindata[i])
print("error: ", waveleterr/waveletn)

當waveletnc = 20時,可得到下圖,誤差大約為2.1

詳解python實現(xiàn)小波變換的一個簡單例子

當waveletnc = 100時,則為下圖,誤差大約為0.04

詳解python實現(xiàn)小波變換的一個簡單例子

當waveletnc = 200時,得到下圖,誤差大約為0.0005

詳解python實現(xiàn)小波變換的一個簡單例子

 以上就是本文的全部內(nèi)容,希望對大家的學習有所幫助,也希望大家多多支持億速云。

向AI問一下細節(jié)

免責聲明:本站發(fā)布的內(nèi)容(圖片、視頻和文字)以原創(chuàng)、轉(zhuǎn)載和分享為主,文章觀點不代表本網(wǎng)站立場,如果涉及侵權(quán)請聯(lián)系站長郵箱:is@yisu.com進行舉報,并提供相關(guān)證據(jù),一經(jīng)查實,將立刻刪除涉嫌侵權(quán)內(nèi)容。

AI