python - 时变频率的频谱图与比例图

Question

我正在比较 FFT 与 CWT 的特定信号。

我不清楚如何从 CWT 的相应比例图中读取相应的频率和幅度。此外，我的印象是 CWT 很不精确？

频谱图在预测精确频率方面似乎相当不错，但对于 CWT，我尝试了许多不同的小波，结果是相同的。

我监督了什么？这不是解决这个问题的合适工具吗？

下面，您可以找到我的示例源代码和相应的绘图。

import matplotlib.pyplot as plt
import numpy as np
from numpy import pi as π
from scipy.signal import spectrogram
import pywt

f_s = 200              # Sampling rate = number of measurements per second in [Hz]
t   = np.arange(-10,10, 1 / f_s) # Time between [-10s,10s].
T1  = np.tanh(t)/2  + 1.0 # Period in [s]
T2  = 0.125               # Period in [s]
f1  = 1 / T1              # Frequency in [Hz]
f2  = 1 / T2              # Frequency in [Hz] 

N = len(t)
x = 13 * np.sin(2 * π * f1 * t) + 42 * np.sin(2 * π * f2 * t)

fig, (ax1, ax2, ax3, ax4) = plt.subplots(4,1, sharex = True, figsize = (10,8))

# Signal
ax1.plot(t, x)
ax1.grid(True)
ax1.set_ylabel("$x(t)$")
ax1.set_title("Signal x(t)")

# Frequency change
ax2.plot(t, f1)
ax2.grid(True)
ax2.set_ylabel("$f_1$ in [Hz]")
ax2.set_title("Change of frequency $f_1(t)$")

# Moving fourier transform, i.e. spectrogram
Δt = 4 # window length in [s]
Nw = np.int(2**np.round(np.log2(Δt * f_s))) # Number of datapoints within window
f, t_, Sxx = spectrogram(x, f_s, window='hanning', nperseg=Nw, noverlap = Nw - 100, detrend=False, scaling='spectrum')
Δf  =  f[1] - f[0]
Δt_ = t_[1] - t_[0]
t2  = t_ + t[0] - Δt_

im = ax3.pcolormesh(t2, f - Δf/2, np.sqrt(2*Sxx), cmap = "inferno_r")#, alpha = 0.5)
ax3.grid(True)
ax3.set_ylabel("Frequency in [Hz]")
ax3.set_ylim(0, 10)
ax3.set_xlim(np.min(t2),np.max(t2))
ax3.set_title("Spectrogram using FFT and Hanning window")

# Wavelet transform, i.e. scaleogram
cwtmatr, freqs = pywt.cwt(x, np.arange(1, 512), "gaus1", sampling_period = 1 / f_s)
im2 = ax4.pcolormesh(t, freqs, cwtmatr, vmin=0, cmap = "inferno" )  
ax4.set_ylim(0,10)
ax4.set_ylabel("Frequency in [Hz]")
ax4.set_xlabel("Time in [s]")
ax4.set_title("Scaleogram using wavelet GAUS1")

# plt.savefig("./fourplot.pdf")

plt.show()

标签： pythonsignal-processingfftwaveletpywavelets

Answer 1

您的示例波形在任何给定时间仅由几个纯音组成，因此频谱图看起来非常干净且可读。

相比之下，小波图看起来“凌乱”，因为您必须对可能不同尺度（以及频率）的高斯小波求和，以近似原始信号中的每个分量纯音。

短时 FFT 和小波变换都属于时频变换类别，但具有不同的核。FFT的内核是纯音，但是您提供的小波变换的内核是高斯小波。FFT 纯音内核将完全对应于您在问题中显示的波形类型，但小波内核不会。

我怀疑你得到的不同小波的结果在数值上完全相同。肉眼看起来可能与您绘制它的方式相同。您似乎将小波用于错误的目的。小波在分析信号方面比绘制信号更有用。小波分析是独一无二的，因为小波组合中的每个数据点同时编码频率、相位和窗口信息。这允许设计位于时间序列和频率分析之间连续统一体的算法，并且可以非常强大。

至于您关于不同小波给出相同结果的说法，这显然不适用于所有小波：我修改了您的代码，它会生成跟随它的图像。当然，GAUS2 和 MEXH 似乎产生了相似的图（放大到非常接近，你会发现它们有细微的不同），但那是因为二阶高斯小波看起来类似于墨西哥帽小波。

import matplotlib.pyplot as plt
import numpy as np
from numpy import pi as π
from scipy.signal import spectrogram, wavelets
import pywt
import random


f_s = 200              # Sampling rate = number of measurements per second in [Hz]
t   = np.arange(-10,10, 1 / f_s) # Time between [-10s,10s].
T1  = np.tanh(t)/2  + 1.0 # Period in [s]
T2  = 0.125               # Period in [s]
f1  = 1 / T1              # Frequency in [Hz]
f2  = 1 / T2              # Frequency in [Hz] 

N = len(t)
x = 13 * np.sin(2 * π * f1 * t) + 42 * np.sin(2 * π * f2 * t)

fig, (ax1, ax2, ax3, ax4) = plt.subplots(4,1, sharex = True, figsize = (10,8))


wvoptions=iter(['gaus1','gaus2','mexh','morl'])

axes=[ax1,ax2,ax3,ax4]


for ax in axes:
    # Wavelet transform, i.e. scaleogram
    try:
        choice=next(wvoptions)
        cwtmatr, freqs = pywt.cwt(x, np.arange(1, 512), choice, sampling_period = 1 / f_s)
        im = ax.pcolormesh(t, freqs, cwtmatr, vmin=0, cmap = "inferno" )  
        ax.set_ylim(0,10)
        ax.set_ylabel("Frequency in [Hz]")
        ax.set_xlabel("Time in [s]")
        ax.set_title(f"Scaleogram using wavelet {choice.upper()}")
    except:
        pass
# plt.savefig("./fourplot.pdf")

plt.tight_layout()
plt.show()