香农熵优化小波包拼接的Python/PyWavelet实现

Question

我需要对一个小波包树运行一个Shannon（或任何类似的度量）熵最小化，但是我在使用算法时遇到了问题。有时这被称为在频率空间中优化磁贴。 我可以解决一个比例尺的问题，但我找不到解决与后续比例尺重叠冲突的方法。是否应该简单地用第一个解析层的结果来填充结果比例图中的间隙？

我试过用PyWavelet的pywt.小波包，并在google上搜索解决这个问题的人，但我似乎找不到基于Python的答案。最接近我的是Matlab在小波工具箱中的besttree函数，但我没有那个toobox的许可证，所以我不能简单地将它MCC并在python中运行。

import pywt
wave='db4'
wp=pywt.WaveletPacket(data,wave)
levels=pywt.dwt_max_level(len(data),pywt.Wavelet(wave))+1

在这一点上，我不知道是沿着一个节点并追溯到它的父节点，还是沿着给定的级别进行跟踪更好。

我一直在使用Addison的“插图小波变换手册”作为技术参考。我试图复制第170页图3.41所示的过程，或类似于以下链接： https://www.researchgate.net/figure/TF-tiling-comparison-between-a-a-DWT-and-b-a-sample-WP-decomposition_fig1_4128902

Answer 1

这似乎有效，但对于大型数据集来说，这是相当缓慢的。
另外，我不完全确定使用自然顺序向上建造是正确的，但它似乎确实同意化学信号谱图（即STFFT）。在

另外，在命令中指定'freq'修剪后尝试返回叶节点会重新填充树。在

下面的函数用于返回时频空间中的标度图，类似于化学信号谱图（即图像状阵列）。在

“freq”用于设置图像的y轴范围。在

tile在某种程度上是一种簿记度量，用来查看相对于STFFT使用了多少个数据点

def Shannon(data):
    if len(data)==1:
        S=data
    else:
        E=data**2/len(data)
        P=E/sum(E)
        S=-sum(P*np.log(P))
    return S
def wpscalogram(Data,rate=5e4,thresh=1.,wave='db4',**_):
    wp=pywt.WaveletPacket(Data,wave)
    wp.get_leaf_nodes(decompose=True)
    levels=pywt.dwt_max_level(len(Data),pywt.Wavelet(wave))

    #DO NOT USE "decompose=True" or "get_level(max_level)" FROM THIS POINT
    for level in range(levels,1,-1):
        print level
        nodes=wp.get_level(level,order='natural',decompose=False)
        paths=[n.path for n in nodes]
        n=len(paths)
        for _i in range(0,n,2):
            Cval=np.hstack([wp[paths[_i]].data,wp[paths[_i+1]].data])
            Pval=wp[wp[paths[_i]].parent.path].data

            if Shannon(Cval)>Shannon(Pval)*thresh:
                wp.__delitem__(paths[_i])
                wp.__delitem__(paths[_i+1])
            else:
                wp[wp[paths[_i]].parent.path].data=min(Shannon(Cval),Shannon(Pval))
    leaves=wp.get_leaf_nodes()
    print [len(leaves[i].path) for i in range(len( leaves))]            
    #ONE CAN NOW DECOMPOSE TREE

    tiles=[len(l.data) for l in leaves]
    col=int(np.max(tiles))
    tiles=sum(tiles)
    freq=np.array([0,0])
    for j,l in enumerate(leaves):
        y=l.data
        level=levels-l.level+1
        if len(y)<col:
            x=col*np.arange(1,len(y)+1).astype(float)/(len(y)+1)
            xi=np.linspace(0,col,col)
            yi=griddata(points=x,values=y,xi=xi,method='nearest')
        else:
            yi=y
        if j==0:
            freq[0]=rate*pywt.scale2frequency(wave,level)
            freq[1]=np.copy(freq[0])
            im=np.matlib.repmat(yi,level,1)
        else:
            im=np.vstack([np.matlib.repmat(yi,level,1),im])
            freq[1]+=rate*pywt.scale2frequency(wave,level)
    print freq
    im[im==0]=np.nan
    return im,freq,tiles