pythonsignal-processingtime-frequency

Plotting Wigner-Ville distribution with pytftb on python 3


I am testing wigner ville distribution to see if it works for the estimation of original amplitude of a signal with noise.

The pytftb provides a wigner ville function that works well with their examples. I use it like this:

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

And it looks like it's working:

enter image description here

However, I can't really understand what's happening here, so i'd like to alter the frequency axis, and probably print regarding to absolute frequency rather than normalized? I don't really understand Wigner Ville Distribution though, so if am looking at it wrong i'd appreciate some help!


Solution

  • Here is some code to plot the Wigner-Ville Transform (WVT) using the absolute frequency. I use the python library: https://pypi.org/project/tftb/, version 0.1.1

    First I create a signal:

    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()
    

    Here is a plot of the signal signal

    For reference, I will also construct a spectrogram before the 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.
    

    Now plotting the heatmaps:

    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()
    

    Here is the result:

    spectrogram and Wigner-Ville transform

    Because the signal is real, there is activity in the positive and in the negative frequencies, just like for the FFT. In addition, there is the interference between the 2 terms, at frequency 0 (i.e., in the middle between the 2). This is also what you saw in your graph. What you had above the middle lines are the negative frequencies and what is at the very top is frequency 0.

    For an analytical signal, the image is simpler. Here I construct the signal using:

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

    And here is the resulting spectrogram and WVT: spectrogram and Wigner-Ville transform of analytic signal

    Let me know if you need any additional clarifications.