信息处理 - python中的尺度图 - 吾爱随笔录

python中的尺度图

信息处理频率小波转换 scipy 英担

2022-01-29 13:11:24

我正在阅读这篇论文来学习 dsp 的基本概念，并且我想重现以下测试信号的尺度图（论文的图 4.2）：

它是由公式的离散化产生的：

\frac{1}{s} \int_{0}^{Ω} f (t) \bar{g (\frac{t - τ}{s})} d t

$\frac{1}{s} \int_{0}^{\Omega} f(t) \overline{g\left(\frac{t-\tau}{s}\right)} d t$

其中是小波： $g(t)$

g (t) = w^{- 1} e^{- π (t / w)^{2}} e^{i 2 π η t / w}

$g(t)=w^{-1} e^{-\pi(t / w)^{2}} e^{i 2 \pi \eta t / w}$

是函数： \ $f(t)$

\begin{array}{l} \sin (2 π ν_{1} t) e^{- π [(t - 0.2) / 0.1]^{10}} \\ + [\sin (2 π ν_{1} t) + 2 \cos (2 π ν_{2} t)] e^{- π [(t - 0.5) / 0.1]^{10}} \\ + [2 \sin (2 π ν_{2} t) - \cos (2 π ν_{3} t)] e^{- π [(t - 0.8) / 0.1]^{10}} \end{array}

$\begin{array}{l} \sin \left(2 \pi \nu_{1} t\right) e^{-\pi[(t-0.2) / 0.1]^{10}} \\ \quad+\left[\sin \left(2 \pi \nu_{1} t\right)+2 \cos \left(2 \pi \nu_{2} t\right)\right] e^{-\pi[(t-0.5) / 0.1]^{10}} \\ \quad+\left[2 \sin \left(2 \pi \nu_{2} t\right)-\cos \left(2 \pi \nu_{3} t\right)\right] e^{-\pi[(t-0.8) / 0.1]^{10}} \end{array}$

我尝试使用该功能重现它们scipy.signal.ctw()。我已将编码为 t 和的函数。尽管如此，我只得到了波浪最后一部分的非零系数，我想我可能有问题。如何使用计算尺度图？ $g(t)$ $\omega$ scipy.signal.ctw()

1个回答

scipy'scwt是原始且容易出错的；下面是通过ssqueezepy.cwt：

代码：

请注意，如果您寻求自己对小波进行编码，则需要先将其带到频域（最好通过傅里叶变换进行分析，然后对其进行采样，而不是通过 FFT），然后像wavelet = Wavelet(my_func); cwt(x, wavelet).

import numpy as np
from ssqueezepy import cwt
from ssqueezepy.visuals import plot, imshow

#%%# Helper fn + params #####################################################
def exp_am(t, offset):
    return np.exp(-pi*((t - offset) / .1)**10)

pi = np.pi
v1, v2, v3 = 64, 128, 32

#%%# Make `x` & plot #########################################################
t = np.linspace(0, 1, 2048, 1)
x = (np.sin(2*pi * v1 * t) * exp_am(t, .2) +
     (np.sin(2*pi * v1 * t) + 2*np.cos(2*pi * v2 * t)) * exp_am(t, .5)  + 
     (2*np.sin(2*pi * v2 * t) - np.cos(2*pi * v3 * t)) * exp_am(t, .8))
plot(x, title="x(t) | t=[0, ..., 1], %s samples" % len(x), show=1)

#%%# Take CWT & plot #########################################################
Wx, scales = cwt(x, 'morlet')
imshow(Wx, yticks=scales, abs=1,
       title="abs(CWT) | Morlet wavelet",
       ylabel="scales", xlabel="samples")

其它你可能感兴趣的问题

上一篇GNU Radio Companion - 生成 4-QAM 星座 - 怎么做？下一篇符号时间偏移模拟