import numpy as np
import scipy.signal
import matplotlib.pyplot as mp
import matplotlib.mlab as mm

# Create a function for the instantaneous frequency ...
dt = 0.01
t = np.arange(0,20,dt)
f = (0.2 + 1.3*np.exp(-((t-10)/4)**2))

# ... and envelope
env = 2-f/2

#  Plot the corresponding function values
df = 2*np.pi*f*dt
fsum = np.cumsum(df)
x = np.sin(fsum)*env

#  Using the Hilbert transform, find the envelope and zero crossings
envelope = abs(scipy.signal.hilbert(x))
phase = np.angle(scipy.signal.hilbert(x))
zeroCrossing = mm.find(np.diff(np.sign(np.cos(phase)))==2)

# Plot the results
mp.plot(x)
mp.hold('on')
mp.plot(envelope, 'r')
mp.plot((0, 2000), (0,0), 'k--')
mp.plot(zeroCrossing, np.zeros(zeroCrossing.size), 'ko')
mp.ylim((-2, 2.2))
mp.savefig('hilbert.eps')
mp.show()