Реализация на Python / PyWavelet оптимизированного для Шеннона энтропии тайлинга WaveletPacket

Question

Мне нужно запустить минимизацию энтропии Шеннона (или любую другую метрику) для дерева вейвлет-пакетов, но у меня возникают проблемы с выяснением алгоритма.Иногда это называется оптимизацией тайлов в частотном пространстве.Я могу последовать решению проблемы для шкалы, но я не могу найти способ разрешения конфликтов перекрытия с последующими шкалами.Должны ли пробелы в результирующей скалограмме просто дополняться результатом первого разрешенного слоя?

Я пытался использовать pywt.WaveletPacket из PyWavelet и отыскивать Google для кого-то, кто решил эту проблему, но я не могу 'Кажется, я не нашел ответ на основе Python.Самым близким из всех, что я получил, была функция лучшего дерева Matlab в Wavelet Toolbox, но у меня нет лицензии на этот тобокс, поэтому я не могу просто MCC-it и запустить его на python.

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

На этом этапе, Я не знаю, лучше ли идти по узлу и проследить до его родителей или проследить по заданному уровню.

Я использовал «Справочник по иллюстрированному вейвлет-преобразованию» Аддисона для технической справки,Я пытаюсь повторить процедуру, показанную на рис. 3.41 на стр. 170, или похожую на следующую ссылку: https://www.researchgate.net/figure/TF-tiling-comparison-between-a-a-DWT-and-b-a-sample-WP-decomposition_fig1_4128902

Answer 1

Кажется, это работает, но это довольно медленно для большого набора данных.
Кроме того, я не совсем уверен, что построение вверх с использованием естественного порядка является правильным, но, похоже, оно действительно согласуется с scipy.signal.spectrogram (т.е. STFFT).

Кроме того, попытка вернуть листовые узлы после обрезки, указав 'freq' в команде, заполняет дерево.

Функция, приведенная ниже, предназначена для возврата скалограммы в частотно-временном пространстве, аналогичной scipy.signal.spectrogram (то есть массиву, подобному изображению).

«freq» используется для установки диапазона оси Y для изображения.

плитка является своего рода мерой бухгалтерского учета, чтобы увидеть, сколько точек данных используется относительно 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