Python/PyWavelet implementation of Shannon Entropy optimized WaveletPacket Tiling

Question

I need to run a Shannon (or any such metric) Entropy minimization for a wavelet packet tree, but I'm having problems sussing out the algorithm. Sometimes this is referred to as optimizing the tiles in frequency-space. I can follow solve the problem for a scale, but I can't find a way of resolving overlap conflicts with the subsequent scales. Should the gaps in the resultant scalogram simply be padded with the result from the first resolved layer?

I've tried using PyWavelet's pywt.WaveletPacket, and scouring google for someone that has addressed this problem, but I can't seem to find a Python based answer. The closest I got was Matlab's besttree function in the Wavelet Toolbox, but I don't have a license to that toobox so I can't simply MCC-it and run in python.

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

At that point, I don't know if it's better to proceed along a node and trace back up to its parents, or trace along a given level.

I've been using Addison's "Illustrated Wavelet Transform Handbook" for technical reference. I'm trying to replicate the procedure shown in Figure 3.41 on page 170, or similar to the following link: https://www.researchgate.net/figure/TF-tiling-comparison-between-a-a-DWT-and-b-a-sample-WP-decomposition_fig1_4128902

Answer 1

This seems to work, but it's quite slow for a large data set.
Also, I'm not entirely sure building upwards using natural order is the correct, but it does seem to agree with scipy.signal.spectrogram (i.e. STFFT).

Additionally, trying to return the leaf nodes after pruning specifying 'freq' in the command repopulates the tree.

The function below is intended to return a scalogram in time-frequnecy space similar to scipy.signal.spectrogram (i.e. image-like array).

"freq" is used to set the y-axis range for the image.

tile is somewhat of a book-keeping measure to see how many data points are used relative to a 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