香农熵优化WaveletPacket拼接的Python/PyWavelet实现
我需要对小波包树运行香农(或任何类似的度量)熵最小化，但我在找出算法时遇到了问题。有时这被称为在频率空间中优化瓦片。我可以跟随解决一个规模的问题，但我找不到一种方法来解决与后续规模的重叠冲突。生成的比例图中的间隙是否应该简单地用第一个解析层的结果填充？
我试过使用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的“图解小波变换手册”作为技术参考。我正在尝试复制170页的图3.41中所示的过程，或类似于以下链接：https://www.researchgate.net/figure/TF-tiling-comparison-between-a-a-DWT-and-b-a-sample-WP-decomposition_fig1_4128902
回答 1
Stack Overflow用户
发布于 2019-05-10 03:47:26
这似乎是有效的，但对于大型数据集来说，它相当慢。
此外，我不完全确定使用自然顺序向上构建是正确的，但它似乎确实符合scipy.signal.spectrogram (即短时傅立叶变换)。
此外，尝试在修剪后返回叶节点，并在命令中指定'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,tileshttps://stackoverflow.com/questions/55923335
复制相似问题