libtfr example for the new libtfr API (develop branch)

libtfr_test.py

import sys
import matplotlib.pyplot as plt
import numpy as np
import libtfr

def fmsin(N, fnormin=0.05, fnormax=0.45, period=None, t0=None, fnorm0=0.25, pm1=1):
    """
    Signal with sinusoidal frequency modulation.

    generates a frequency modulation with a sinusoidal frequency.
    This sinusoidal modulation is designed such that the instantaneous
    frequency at time T0 is equal to FNORM0, and the ambiguity between
    increasing or decreasing frequency is solved by PM1.

    N       : number of points.
    FNORMIN : smallest normalized frequency          (default: 0.05)
    FNORMAX : highest normalized frequency           (default: 0.45)
    PERIOD  : period of the sinusoidal fm            (default: N   )
    T0      : time reference for the phase           (default: N/2 )
    FNORM0  : normalized frequency at time T0        (default: 0.25)
    PM1     : frequency direction at T0 (-1 or +1)       (default: +1  )

    Returns:
    Y       : signal
    IFLAW   : its instantaneous frequency law

    Example:
    z,i=fmsin(140,0.05,0.45,100,20,0.3,-1.0)

    Original MATLAB code F. Auger, July 1995.
    (note: Licensed under GPL; see main LICENSE file)
    """

    if period==None:
        period = N
    if t0==None:
        t0 = N/2
    pm1 = np.sign(pm1)

    fnormid=0.5*(fnormax+fnormin);
    delta  =0.5*(fnormax-fnormin);
    phi    =-pm1*np.arccos((fnorm0-fnormid)/delta);
    time   =np.arange(1,N)-t0;
    phase  =2*np.pi*fnormid*time+delta*period*(np.sin(2*np.pi*time/period+phi)-np.sin(phi));
    y      =np.exp(1j*phase)
    iflaw  =fnormid+delta*np.cos(2*np.pi*time/period+phi);

    return y,iflaw

if __name__=="__main__":
    N = 256
    NW = 3.5
    step = 10
    k = 6
    tm = 6.0
    Np = 201

    nloop = 1 if len(sys.argv) == 1 else int(sys.argv[1])

    # generate a nice dynamic signal
    siglen = 17590
    ss,iff = fmsin(siglen, .15, 0.45, 1024, 256/4,0.3,-1)
    signal = ss.real + np.random.randn(ss.size)*.1

    # generate a transform object with size equal to signal length and 5 tapers
    D = libtfr.mfft_dpss(Np, 3, 5)
    # complex windowed FFTs, one per taper
    Z = D.mtfft(signal)
    # power spectrum density estimate using adaptive method to average tapers
    P = D.mtpsd(signal)

    # generate a transform object with size 512 samples and 5 tapers for short-time analysis
    D = libtfr.mfft_dpss(512, 3, 5)
    # complex STFT with frame shift of 10 samples
    Z = D.mtstft(signal, 10)
    # spectrogram with frame shift of 10 samples
    P = D.mtspec(signal, 10)
    plt.imshow(P, interpolation='nearest', cmap='viridis', aspect='auto')
    plt.show()
    #print (P)

    # generate a transform object with 200-sample hanning window padded to 256 samples
    D = libtfr.mfft_precalc(256, np.hanning(200))
    # spectrogram with frame shift of 10 samples
    P = D.mtspec(signal, 10)
    plt.imshow(P, interpolation='nearest', cmap='viridis', aspect='auto')
    plt.show()