gen_tsp2.py

#!/usr/bin/python3
import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack


def gen_tsp2(T=1200, fs=44.1, bw=22.05, bs=0, ta=120, tb=50, plot=False):
    '''
    GEN_TSP2: function to generate and plot TSP signals

    TSP stands for "time-stretched pulse". A TSP is a pulse whose magnitude
    spectrum is flat, as in the case of the spectra of an impulse and of
    white noise. The difference is in the phase spectrum, which has a quadratic
    shape in the TSP case, in contrast with the linear phase observed in the
    impulse case and the random phase observed in the white noise case.

    The group delay spectrum, defined as "minus" the derivative of the phase
    spectrum with respect to the frequency, is therefore linear. This means
    that, in a TSP, each frequency arrives with a delay that changes linearly
    with the value of that frequency.

    The sintax of the function gen_tsp2 is as follows:

    gen_tsp2(T: float, fs: float, bw: float, bs: float, ta: float, tb: float)
        => (x: np.ndarray, t: np.ndarray)

    x  -> TSP generated by gen_tsp2
    t  -> time vector associated with x

    T  -> duration of the TSP. The default is 1200ms.
    fs -> sampling frequency (in kHz). The default is 44.1kHz.
    bw -> bandwidth (in kHz) of the TSP. The default is the Nyquist frequency,
          i.e. half of the sampling frequency.
    bs -> bandshift (in kHz) of the TSP. The default is 0kHz.
    ta -> starting time of the TSP (in ms). The default is T/10 ms.
    tb -> group delay growing rate (in ms/kHz). The default is 50ms/kHz.
          Note that, if a TSP with monotonically growing (or decreasing)
          frequency, is desired, the condition T >= ta + tb*(bs+bw) must
          be satisfied.

          Obs. The duration of the TSP usually becomes longer than T in order
          to make the number of samples in x equal to a power of 2. This speeds
          up the computation of the discrete Fourier transform used in the
          procedure.

          (Function written by Hani Yehia 31/10/1999)
    '''

    # Modify ta and tb so that the TSP phase at the Nyquist Frequency
    # (i.e. fs/2) is a multiple of pi. (This is necessary to make the TSP in the
    # time domain a real signal)
    ta = round(ta*fs / fs)
    tb = round(tb*fs**2/8) * 8 / fs**2

    # Set the FFT length as a power of 2
    fft_length = int(2**(np.ceil(np.log2(fs * T))))
    T = fft_length / fs

    # Set frequency step and angular frequency values
    step_f = fs / fft_length
    w1 = 2*np.pi * np.arange(0, fs/2 + step_f, step_f)

    # Make the magnitude spectrum of the TSP
    if (bw == fs / 2):
        Mag = np.ones((fft_length))
        Mag1 = Mag[:len(w1)]
    else:
        power = 12
        Mag1 = np.exp(- ((w1 / (2*np.pi) - (bw/2 + bs)) / (bw/2)) ** power)
        Mag2 = Mag1[len(w1)-1:-1:2]
        Mag = np.hstack((Mag1, Mag2))

    # Make the group delay spectrum of the TSP
    group_delay1 = ta + tb * w1/(2*np.pi)
    group_delay = np.hstack((group_delay1, -group_delay1[-1:-1:2]))

    # Make the phase spectrum of the TSP
    Ph1 = -(ta + tb * w1/(4*np.pi)) * w1
    Ph = np.hstack((Ph1, -Ph1[-2:0:-1]))

    # Make the complex spectrum of the TSP
    X = Mag * np.exp(1j * Ph)

    # Make the TSP in the time domain
    x = np.real(scipy.fftpack.ifft(X))  # imag(ifft(X)) ~= 0 due to finite numerical precision
    x = x / np.max(np.abs(x)) * .98  # normalization

    # Set time step and time vector
    step_t = T / fft_length
    t = np.arange(0, T, step_t)

    if plot:
        # Plot a figure showing the main characteristics of the TSP generated
        fig = plt.figure('TSP')
        fig.suptitle('TSP', fontsize=20)
        plt.tight_layout()

        ax = fig.add_subplot(2, 2, 1)
        ax.plot(w1[::int(np.ceil(fft_length / 256))] / (2*np.pi),
                Mag1[::int(np.ceil(fft_length / 256))])
        ax.set_xlim([0, fs/2])
        # ax.set_xticks([np.arange(0, 30, np.round(fs / 8))])
        ax.set_ylim([0, 1.25])
        # ax.set_yticks([np.arange(0, 1, .1)])
        ax.set_xlabel('Frequency (kHz)')
        ax.set_ylabel('Normalized Amplitude')
        ax.set_title('TSP Magnitude Spectrum')
        ax.text(x=.1, y=.95, s='Passband = {:.4f}--kHz'.format(bs))

        ax = fig.add_subplot(2, 2, 2)
        ax.plot(w1[::int(np.ceil(fft_length / 256))] / (2*np.pi),
                group_delay[::int(np.ceil(fft_length / 256))])
        ax.set_xlim([0, fs/2])
        ax.set_ylim(0, group_delay[-1] * 1.2)
        ax.set_xlabel('Frequency (kHz)')
        ax.set_ylabel('Group Delay (ms)')
        ax.set_title('TSP Group Delay Spectrum')
        ax.text(x=.1, y=.9, s='t0 = {:.4f}ms'.format(ta))
        ax.text(x=.1, y=.8, s='t1 = {:.4f}ms/kHz'.format(tb))

        ax = fig.add_subplot(2, 1, 2)
        ax.plot(t, x)
        ax.set_xlim([0, T])
        ax.set_ylim([1.1*min(x), 1.1*max(x)])
        ax.set_xlabel('Time (ms)')
        ax.set_ylabel('Amplitude')
        ax.set_title('TSP Signal')

        plt.show()

    return x, t