def theta_t(f, t):  # compute complete input ("phase") from instantaneous freq
    theta0 = f[0] * t[0]  # initial phase / input
    theta = [theta0]
    for i in range(1, len(f)):
        theta.append(theta[-1] + (t[i] - t[i - 1]) * f[i])
    return 2 * np.pi * np.array(theta)
    # Equivalent:
    # return 2 * np.pi * np.cumsum([f[0] * t[0]] + list(f[1:] * np.diff(t)))