0

我需要为小波包树运行香农（或任何此类度量）熵最小化，但我在找出算法时遇到问题。有时这被称为优化频率空间中的瓦片。我可以按照比例解决问题，但我找不到解决与后续比例重叠冲突的方法。生成的尺度图中的空白是否应该简单地用第一个解析层的结果填充？

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

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

那时，我不知道是沿着一个节点继续并追溯到它的父节点，还是沿着给定的级别追踪是否更好。

我一直在使用 Addison 的“Illustrated Wavelet Transform Handbook”作为技术参考。我正在尝试复制第 170 页上的图 3.41 中显示的过程，或类似于以下链接： https ://www.researchgate.net/figure/TF-tiling-comparison-between-aa-DWT-and-ba -sample-WP-decomposition_fig1_4128902

4

1 回答 1

0

这似乎可行，但对于大型数据集来说速度很慢。
另外，我不完全确定使用自然顺序向上构建是正确的，但它似乎与 scipy.signal.spectrogram （即 STFFT）一致。

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

下面的函数旨在返回类似于 scipy.signal.spectrogram 的时频空间中的尺度图（即类似图像的数组）。

“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
于 2019-05-09T19:47:26.103 回答