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

Question

我正在测试 wigner ville 分布，看看它是否适用于估计带有噪声的信号的原始振幅。

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

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

看起来它正在运行：

但是，我真的不明白这里发生了什么，所以我想改变频率轴，并可能打印关于绝对频率而不是标准化的？不过，我不太了解 Wigner Ville Distribution，所以如果我看错了，我将不胜感激！

Answer 1

这里是一些使用绝对频率绘制 Wigner-Ville 变换 (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()

结果如下：

因为信号是真实的，所以在正频率和负频率中有 activity，就像 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：

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

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

Plotting Wigner-Ville distribution with pytftb on python 3

python

signal-processing

time-frequency