In [ ]:

    

import numpy as np
import scipy as sp
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")

import libtfr
import pointproc


    

In [ ]:

    

def fmsin(N, fnormin, fnormax, period, t0, fnorm0, pm1):
    from numpy import sign, sin, arccos, exp, arange, pi, real
    fnormid = 0.5 * (fnormax+fnormin)
    delta = 0.5 * (fnormax-fnormin)
    phi = - sign(pm1) * arccos((fnorm0 - fnormid)/delta)
    t = arange(0.0, N)
    phase = 2 * pi * fnormid * (t - t0) + delta * period * (sin(2 * pi * (t - t0) / period + phi) - sin(phi))
    return real(exp(1j * phase))

Fs = 1000.
sig = fmsin(17590, 0.15, 0.45, 1024, 256./4, 0.3, -1)

# an inhomogeneous poisson process with rate = exp(sig). Times will be in units of s
p = np.exp(sig - 1)
events = (p > np.random.uniform(size=p.size)).nonzero()[0].astype('d') / Fs
# jitter 
events += np.random.uniform(low=-0.00025, high=0.00025, size=events.size)


    

In [ ]:

    

nfft = sig.size
ntapers = 5
transform = libtfr.mfft_dpss(nfft, 3, ntapers, nfft)


    

In [ ]:

    

# transform a time series
Z = transform.mtpsd(sig)
f, idx = libtfr.fgrid(Fs, transform.nfft)
plt.plot(f, Z[idx])


    

In [ ]:

    

# point process fft - should look like noisy version of signal spectrum
JJ = transform.mtfft_pt(events, 1 / Fs,  0)
SS = np.mean(JJ * JJ.conj(), 1).real
plt.plot(SS)


    

In [ ]:

    

# spectrograms
sns.set_style("dark")
nfft = 256
window = 201
shift = 10
ntapers = 5
nframes = (sig.size - nfft) // shift + 1
D = libtfr.mfft_dpss(nfft, 3, ntapers, window)
Z = D.mtspec(sig, shift)
plt.imshow(Z, cmap="jet")


    

In [ ]:

    

Z, Nsp = D.mtstft_pt(events, 1 / Fs, shift / Fs, 0, sig.size / Fs)
S = (Z * Z.conj()).mean(2).real
plt.imshow(S, cmap="jet")


    

In [ ]: