在 python 3 上使用 pytftb 绘制 Wigner-Ville 分布

Question

我正在测试维格纳维尔分布，看看它是否适用于估计带有噪声的信号的原始幅度。

pytftb 提供了一个 wigner ville 函数，可以很好地配合他们的示例。我这样使用它：

tfr = WignerVilleDistribution(prepack[0])
tfr.run()
tfr.plot(show_tf=True)

看起来它正在工作：

但是，我无法真正理解这里发生的情况，所以我想更改频率轴，并且可能打印绝对频率而不是标准化频率？我不太了解维格纳维尔分布，所以如果我看错了，我会很感激一些帮助！

Answer 1

下面是一些使用绝对频率绘制维格纳-维尔变换 (WVT) 的代码。我使用python库：https://pypi.org/project/tftb/，版本0.1.1

首先我创建一个信号：

import tftb
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sig


T = 2  # signal duration
dt = 1/500  # sample interval/spacing
freq_s = 1/dt  # sampling frequency
N = T / dt  # number of samples
ts = np.arange(N) * dt  # times

#  constructing a chirp multiplied by a Gaussian
t0 = T/2
freq = np.linspace(10, 30, int(N))
sigma = 0.1
signal = np.cos((ts-t0) * 2 * np.pi * freq) * np.exp(-(ts-t0)**2/(2*sigma**2))/np.sqrt(sigma)

# adding some noise
signal += np.random.randn(len(signal))*0.5


#  plotting the signal
plt.figure()
plt.plot(ts, signal)
plt.show()

这是信号图

作为参考，我还将在 WVT 之前构建一个频谱图：

# first looking at the power of the short time fourier transform (SFTF):
nperseg = 2**6  # window size of the STFT
f_stft, t_stft, Zxx = sig.stft(signal, freq_s, nperseg=nperseg,
                           noverlap=nperseg-1, return_onesided=False)

# shifting the frequency axis for better representation
Zxx = np.fft.fftshift(Zxx, axes=0)
f_stft = np.fft.fftshift(f_stft)

# Doing the WVT
wvd = tftb.processing.WignerVilleDistribution(signal, timestamps=ts)
tfr_wvd, t_wvd, f_wvd = wvd.run()
# here t_wvd is the same as our ts, and f_wvd are the "normalized frequencies"
# so we will not use them and construct our own.

现在绘制热图：

f, axx = plt.subplots(2, 1)

df1 = f_stft[1] - f_stft[0]  # the frequency step
im = axx[0].imshow(np.real(Zxx * np.conj(Zxx)), aspect='auto',
          interpolation=None, origin='lower',
          extent=(ts[0] - dt/2, ts[-1] + dt/2,
                  f_stft[0] - df1/2, f_stft[-1] + df1/2))
axx[0].set_ylabel('frequency [Hz]')
plt.colorbar(im, ax=axx[0])
axx[0].set_title('spectrogram')


# because of how they implemented WVT, the maximum frequency is half of
# the sampling Nyquist frequency, so 125 Hz instead of 250 Hz, and the sampling
# is 2 * dt instead of dt
f_wvd = np.fft.fftshift(np.fft.fftfreq(tfr_wvd.shape[0], d=2 * dt))
df_wvd = f_wvd[1]-f_wvd[0]  # the frequency step in the WVT
im = axx[1].imshow(np.fft.fftshift(tfr_wvd, axes=0), aspect='auto', origin='lower',
       extent=(ts[0] - dt/2, ts[-1] + dt/2,
               f_wvd[0]-df_wvd/2, f_wvd[-1]+df_wvd/2))
axx[1].set_xlabel('time [s]')
axx[1].set_ylabel('frequency [Hz]')
plt.colorbar(im, ax=axx[1])
axx[1].set_title('Wigner-Ville Transform')
plt.show()

结果如下：

因为信号是真实的，所以正频率和负频率都有活动，就像 FFT 一样。此外，这两个项之间存在频率为 0 的干扰（即，位于两个项之间的中间）。这也是您在图表中看到的。中线上方的是负频率，最顶部的是频率 0。

对于分析信号，图像更简单。在这里，我使用以下方法构建信号：

signal = np.exp(1j * (ts-t0) * 2 * np.pi * freq) * np.exp(-(ts-t0)**2/(2*sigma**2))/np.sqrt(sigma)

这是生成的频谱图和 WVT：

如果您需要任何其他说明，请告诉我。