Python中 pywt 小波分析库中的一些基本使用方法

浅浅的花香味﹌ 2022-05-10 14:58 517阅读 0赞

尺度函数 ： scaling function （在一些文档中又称为父函数 father wavelet ）  
小波函数 ： wavelet function（在一些文档中又称为母函数 mother wavelet）  
连续的小波变换 ：CWT  
离散的小波变换 ：DWT  
小波变换的基本知识：  
不同的小波基函数，是由同一个基本小波函数经缩放和平移生成的。  
小波变换是将原始图像与小波基函数以及尺度函数进行内积运算,所以一个尺度函数和一个小波基函数就可以确定一个小波变换  
小波变换后低频分量  
基本的小波变换函数

### 二维离散小波变换： ###

（官网上的例子）

import numpy as np
    import pandas as pd 
    import matplotlib.pyplot as plt 
    
    def test_pywt():
        import numpy as np
        import matplotlib.pyplot as plt
        import pywt
        import pywt.data
        
        # Load image
        original = pywt.data.camera()
        
        # Wavelet transform of image, and plot approximation and details
        titles = ['Approximation', ' Horizontal detail',
                  'Vertical detail', 'Diagonal detail']
        coeffs2 = pywt.dwt2(original, 'bior1.3')
        LL, (LH, HL, HH) = coeffs2
        plt.imshow(original)
        plt.colorbar(shrink=0.8)
        fig = plt.figure(figsize=(12, 3))
        for i, a in enumerate([LL, LH, HL, HH]):
            ax = fig.add_subplot(1, 4, i + 1)
            ax.imshow(a, interpolation="nearest", cmap=plt.cm.gray)
            ax.set_title(titles[i], fontsize=10)
            ax.set_xticks([])
            ax.set_yticks([])
        
        fig.tight_layout()
        plt.show()
    # test_pywt()

#coding=gbk
    '''
    Created on 2018年10月9日
    这个模块是为了测试 pywt 库 的相关用法
    @author: Administrator
    '''
    
    import pywt 
    print(pywt.families())  #打印出小波族
    # ['haar', 'db', 'sym', 'coif', 'bior', 'rbio', 'dmey', 'gaus', 'mexh', 'morl', 'cgau', 'shan', 'fbsp', 'cmor']
    
    for family in pywt.families():  #打印出每个小波族的每个小波函数
        print('%s family: '%(family) + ','.join(pywt.wavelist(family)))
    # haar family: haar
    # db family: db1,db2,db3,db4,db5,db6,db7,db8,db9,db10,db11,db12,db13,db14,db15,db16,db17,db18,db19,db20,db21,db22,db23,db24,db25,db26,db27,db28,db29,db30,db31,db32,db33,db34,db35,db36,db37,db38
    # sym family: sym2,sym3,sym4,sym5,sym6,sym7,sym8,sym9,sym10,sym11,sym12,sym13,sym14,sym15,sym16,sym17,sym18,sym19,sym20
    # coif family: coif1,coif2,coif3,coif4,coif5,coif6,coif7,coif8,coif9,coif10,coif11,coif12,coif13,coif14,coif15,coif16,coif17
    # bior family: bior1.1,bior1.3,bior1.5,bior2.2,bior2.4,bior2.6,bior2.8,bior3.1,bior3.3,bior3.5,bior3.7,bior3.9,bior4.4,bior5.5,bior6.8
    # rbio family: rbio1.1,rbio1.3,rbio1.5,rbio2.2,rbio2.4,rbio2.6,rbio2.8,rbio3.1,rbio3.3,rbio3.5,rbio3.7,rbio3.9,rbio4.4,rbio5.5,rbio6.8
    # dmey family: dmey
    # gaus family: gaus1,gaus2,gaus3,gaus4,gaus5,gaus6,gaus7,gaus8
    # mexh family: mexh
    # morl family: morl
    # cgau family: cgau1,cgau2,cgau3,cgau4,cgau5,cgau6,cgau7,cgau8
    # shan family: shan
    # fbsp family: fbsp
    # cmor family: cmor
    
    db3 = pywt.Wavelet('db3')   #创建一个小波对象
    print(db3)
    # Filters length: 6        #滤波器长度
    # Orthogonal:     True    #正交
    # Biorthogonal:   True    #双正交
    # Symmetry:       asymmetric    #对称性，不对称
    # DWT:            True    #离散小波变换
    # CWT:            False    #连续小波变换
    
    def print_array(arr):
        print('[%s]'%','.join(['%.14f'%x for x in arr]))
        
    #离散小波变换的小波滤波系数
    # dec_lo Decomposition filter values 分解滤波值， rec 重构滤波值
    #db3.filter_bank 返回4 个属性
    print(db3.filter_bank == (db3.dec_lo, db3.dec_hi, db3.rec_lo, db3.rec_hi))  #True
    print(db3.dec_len) 
    print(db3.rec_len)  #6 
    
    # DWT 与 IDWT 
    #使用db2 小波函数做dwt
    x = [3, 7, 1, 1, -2, 5, 4, 6]
    cA, cD = pywt.dwt(x, 'db2') #得到近似值和细节系数
    print(cA)   # [5.65685425 7.39923721 0.22414387 3.33677403 7.77817459]
    print(cD)   # [-2.44948974 -1.60368225 -4.44140056 -0.41361256  1.22474487]
    
    #IDWT 
    print(pywt.idwt(cA, cD, 'db2')) # [ 3.  7.  1.  1. -2.  5.  4.  6.]
    
    #传入小波对象，设置模式
    w = pywt.Wavelet('sym3')
    cA, cD = pywt.dwt(x, wavelet=w, mode='constant')
    print(cA)   # [ 4.38354585  3.80302657  7.31813271 -0.58565539  4.09727044  7.81994027]
    print(cD)   # [-1.33068221 -2.78795192 -3.16825651 -0.67715519 -0.09722957 -0.07045258]
    
    print(pywt.Modes.modes)
    # ['zero', 'constant', 'symmetric', 'periodic', 'smooth', 'periodization', 'reflect', 'antisymmetric', 'antireflect']
    
    print(pywt.idwt([1,2,0,1], None, 'db3', 'symmetric'))
    print(pywt.idwt([1,2,0,1], [0,0,0,0], 'db3', 'symmetric'))
    # [ 0.83431373 -0.23479575  0.16178801  0.87734409]
    # [ 0.83431373 -0.23479575  0.16178801  0.87734409]
    
    #小波包 wavelet packets 
    X  = [1, 2, 3, 4, 5, 6, 7, 8]
    wp = pywt.WaveletPacket(data=X, wavelet='db3', mode='symmetric', maxlevel=3)
    print(wp)
    print(wp.data)  #[1 2 3 4 5 6 7 8 9]
    print(repr(wp.path))
    print(wp.level) # 0    #分解级别为0 
    print(wp['ad'].maxlevel)    # 3 
    
    #访问小波包的子节点
    #第一层：
    print(wp['a'].data)
    # [ 4.52111203  1.54666942  2.57019338  5.3986205   8.19182134 11.27067814
    #  12.65348525]        # 当设置分解的 maxlevel 时，分解得到的data
    
    #[ 4.52111203  1.54666942  2.57019338  5.3986205   8.20681003 11.18125264] 设置为2 时
    print(wp['a'].path) # a
    #第2 层
    print(wp['aa'].data)
    # [ 3.63890166  6.00349136  2.89780988  6.80941869 15.41549196]
    print(wp['ad'].data)
    # [ 1.25531439 -0.60300027  0.36403471  0.59368086 -0.53821027]
    print(wp['aa'].path)    # aa
    print(wp['ad'].path)    # ad
    
    #第3 层时：
    print(wp['aaa'].data)
    # [ 6.7736584   5.78857317  5.69392399 10.98672847 19.92241106]
    
    # print(wp['aaaa'].data)  #超过最大层时，会报错
    #获取特定层数的所有节点
    print([node.path for node in wp.get_level(3, 'natural')])   #第3层有8个
    # ['aaa', 'aad', 'ada', 'add', 'daa', 'dad', 'dda', 'ddd']
    
    #依据频带频率进行划分
    print([node.path for node in wp.get_level(3, 'freq')])
    # ['aaa', 'aad', 'add', 'ada', 'dda', 'ddd', 'dad', 'daa']
    
    #从小波包中 重建数据
    X  = [1, 2, 3, 4, 5, 6, 7, 8]
    wp = pywt.WaveletPacket(data=X, wavelet='db1', mode='symmetric', maxlevel=3)
    print(wp['ad'].data)    # [-2,-2]
    new_wp = pywt.WaveletPacket(data=None, wavelet='db1', mode='symmetric')
    new_wp['a'] = wp['a']
    new_wp['aa'] = wp['aa'].data
    new_wp['ad'] = [-2,-2]  # wp['ad'].data
    new_wp['d'] = wp['d']
    print(new_wp.reconstruct(update=False))
    # new_wp['a'] = wp['a']  直接使用高低频也可进行重构
    # new_wp['d'] = wp['d']
    print(new_wp)   #: None
    print(new_wp.reconstruct(update=True))  #更新设置为True时。
    print(new_wp)
    # : [1. 2. 3. 4. 5. 6. 7. 8.]
    
    #获取叶子结点
    print([node.path for node in new_wp.get_leaf_nodes(decompose=False)])
    # ['aa', 'ad', 'd']
    print([node.path for node in new_wp.get_leaf_nodes(decompose=True)])
    # ['aaa', 'aad', 'ada', 'add', 'daa', 'dad', 'dda', 'ddd']
    
    #从小波包树中移除结点
    dummy = wp.get_level(2)
    for i in wp.get_leaf_nodes(False):
        print(i.path, i.data)
    # aa [ 5. 13.]
    # ad [-2. -2.]
    # da [-1. -1.]
    # dd [-1.11022302e-16  0.00000000e+00]
    node = wp['ad']
    print(node) #ad: [-2. -2.]
    del wp['ad'] #删除结点
    for i in wp.get_leaf_nodes(False):
        print(i.path, i.data)
    # aa [ 5. 13.]
    # da [-1. -1.]
    # dd [-1.11022302e-16  0.00000000e+00]    
    
    print(wp.reconstruct()) #进行重建
    # [2. 3. 2. 3. 6. 7. 6. 7.]
    
    wp['ad'].data = node.data   #还原已删除的结点
    print(wp.reconstruct())
    # [1. 2. 3. 4. 5. 6. 7. 8.]
    
    print(wp['a'])
    print(wp.a)
    
        
    
    filename = r'D:\ml_datasets\PHM\c6\c_6_001.csv'
    data = pd.read_csv(filename)
    data = data.iloc[100000:110000, 3]
    cA1, cD1 = pywt.dwt(data, 'db3') #得到近似值和细节系数
    
    wap = pywt.WaveletPacket(data=data, wavelet='db3')
    dataa = wap['a'].data
    print(wap['a'].data)
    print(len(wap['a'].data))
    plt.figure(num='ca')
    plt.plot(dataa)
    plt.figure(num='data')
    plt.plot(dataa)
    plt.show()
    
    
    
    # plt.figure(num='ca')
    # plt.plot(cA1)
    # plt.figure(num='cd')
    # plt.plot(cD1)
    # plt.figure(num='data')
    # plt.plot(data)
    # plt.show()

pywt 库中 upcoef 函数的使用

pywt.upcoef(part='a'or'd', coeffs=cA, cD, wavelet, level, take)

data = np.array([8,9,10,11,1,2,3,4,5,6,7]).reshape(-1,)
    print('origin data:',data)
    (cA, cD) = pywt.dwt(data, 'haar')
    print('cA para:',cA)
    # cA para: [12.02081528 14.8492424   2.12132034  4.94974747  7.77817459  9.89949494]
    print('cD para:',cD)
    print('take length:len(cA)',len(cA))
    # take length:len(cA) 6
    print('take length:len(cD)',len(cD))
    
    print(pywt.upcoef('a', cA, 'haar',take=len(cA)) + pywt.upcoef('d', cD, 'haar',take=len(cD)))
    #截取原始 数据中间的take 个数值
    print('take length:len(data)',len(data))
    print(pywt.upcoef('a', cA, 'haar',take=len(data)) + pywt.upcoef('d', cD, 'haar',take=len(data)))
    print(pywt.upcoef('a', cA, 'haar',take=1) + pywt.upcoef('d', cD, 'haar',take=1))
    print(pywt.upcoef('a', cA, 'haar',take=2) + pywt.upcoef('d', cD, 'haar',take=2))
    print(pywt.upcoef('a', cA, 'haar',take=3))
    # [1.5 1.5 3.5]
    print(pywt.upcoef('d', cD, 'haar',take=3))
    # [-0.5  0.5 -0.5]
    print(pywt.upcoef('a', cA, 'haar',take=3) + pywt.upcoef('d', cD, 'haar',take=3))
    # [1. 2. 3.]
    # 从最中间的位置截取指定长度的数据；如果为1，就是最中间的2；如果为2，就是最中间的两个数，就是2与3；以此类推；

### 项目推荐： ###

**[2000多G的计算机各行业电子资源分享（持续更新）][2000_G]**

**[2020年微信小程序全栈项目之喵喵交友【附课件和源码】][2020]**

**[Spring Boot开发小而美的个人博客【附课件和源码】][Spring Boot]**

**[Java微服务实战296集大型视频-谷粒商城【附代码和课件】][Java_296_-]**

**[Java开发微服务畅购商城实战【全357集大项目】-附代码和课件][Java_357_-]**

**[最全最详细数据结构与算法视频-【附课件和源码】][-]**

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQwNTg3NTc1_size_16_color_FFFFFF_t_70_pic_center]

[2000_G]: https://mp.weixin.qq.com/s/sP4JgGWkCzpgwKr9sAV2_Q
[2020]: https://mp.weixin.qq.com/s?__biz=MzIyNTI3NDQ4NQ==&mid=2247487704&idx=1&sn=5f4b2127c4d49fd07ae072a0721424a2&chksm=e8036fc2df74e6d489c4aa9b06f917ef7cee6027f13e150fca53cf79d5d188d4ccc1af49e098&scene=21#wechat_redirect
[Spring Boot]: https://mp.weixin.qq.com/s?__biz=MzIyNTI3NDQ4NQ==&mid=2247487798&idx=2&sn=ac0293b996521b872a9dba5fbb3e65e6&chksm=e8036e2cdf74e73aba104a9a994a5b2e31483e8dcbe0f1d9936f6d5173b887e1560f59d2819c&scene=21#wechat_redirect
[Java_296_-]: https://mp.weixin.qq.com/s?__biz=MzIyNTI3NDQ4NQ==&mid=2247487674&idx=1&sn=7aff0bdf2bb727303f3d3618995aef21&chksm=e8036fa0df74e6b6d872c7e6ece179c524ed463a4a6b74c96875475c9a3d5ddb903427dd993b&scene=21#wechat_redirect
[Java_357_-]: https://mp.weixin.qq.com/s?__biz=MzIyNTI3NDQ4NQ==&mid=2247486376&idx=1&sn=d1fef270c463ea8ac663f6fbfedd70a0&chksm=e80374b2df74fda4d3bafba878a106a19e18c5fcda266008f4f37975847a21bc612ffcd5ff39&scene=21#wechat_redirect
[-]: https://mp.weixin.qq.com/s?__biz=MzIyNTI3NDQ4NQ==&mid=2247487750&idx=1&sn=747bccbb5f5ea6b58915198de40da777&chksm=e8036e1cdf74e70ae97a5e8e265b49d7236d904d291203309159d07ba1724033062c0e370843&scene=21#wechat_redirect
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzQwNTg3NTc1_size_16_color_FFFFFF_t_70_pic_center]: /images/20220505/09ce5aa92d474d5ea0b4e325a9a16ffb.png