import numpy as np
from scipy.signal import lfilter
import matplotlib.pyplot as plt

# Construct the input signal
idx = np.arange(512)
w0 = [2*np.pi/(80 + x/8) for x in idx]
x0 = np.sin(w0*idx)
# Calculate the transforms
h = hilbert(x0)
I = x0 # normally it should be np.roll(x0, 3)
Q = lfilter([0.25, 0, 0.75, 0, -0.25, 0, -0.75], [1], x0)
# Calculate amplitude arrays
A1 = np.abs(h)
A2 = np.sqrt(I**2 + Q**2)
# Calcuclate phase arrays
ph1 = np.angle(h)
ph2 = -np.arctan2(Q, I)
# Calculate instant. frequency
w1 = np.diff(np.unwrap(ph1))
w2 = np.diff(np.unwrap(ph2))

# Plot the results
fig = plt.figure()
ax = plt.subplot(311)
plt.plot(x0, lw=3, color='b', label='signal')
plt.plot(A1*np.cos(ph1), lw=2, color='g', label='x1')
plt.plot(A1, lw=2, ls='dotted', color='g', label='A1')
plt.plot(A2*np.cos(ph2), lw=1, color='k', label='x2')
plt.plot(A2, lw=1, ls='dotted', color='k', label='A2')
plt.legend(bbox_to_anchor=(1, 1), loc=2, borderaxespad=0.)
plt.subplot(312, sharex=ax)
plt.plot(ph1, lw=2, color='g', label='ph1')
plt.plot(ph2, lw=1, color='k', label='ph2')
plt.legend(bbox_to_anchor=(1, 1), loc=2, borderaxespad=0.)
plt.subplot(313, sharex=ax)
plt.plot(w0, lw=3, color='b', label='w0')
plt.plot(w1, lw=2, color='g', label='w1')
plt.plot(w2, lw=1, color='k', label='w2')
plt.legend(bbox_to_anchor=(1, 1), loc=2, borderaxespad=0.)
fig.show()