bilinear_experiment.py

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

def solve(omega1, omega2, phi1, phi2):
    # Let's formulate the equations accordingly
    A = np.array([
        [-np.sin(phi1 + omega1), 0, np.sin(omega1), 0],
        [-np.cos(phi1 + omega1), 0, np.cos(omega1), 1],
        [0, -np.sin(phi2 + omega2), np.sin(omega2), 0],
        [0, -np.cos(phi2 + omega2), np.cos(omega2), 1],
    ])
    b_ = np.array([-np.sin(2 * omega1),
                   -np.cos(2 * omega1),
                   -np.sin(2 * omega2),
                   -np.cos(2 * omega2)])
    b_ = np.linalg.solve(A, b_)
    return b_

def plot_freqz(a1, a2):
    # Omega = np.arange(0, 10, 0.01)
    # plt.plot(Omega, np.angle(-Omega**2 + a2 + 1j * a1 * Omega))
    # plt.title('Continuous')
    # plt.show()
    #[W, H] = freqz([1,(2*a2*T*T-8) / (a2*T*T+2*a1*T+4),(a2*T*T-2*a1*T + 4) / (a2*T*T+2*a1*T+4)], T*T*np.array([1, 2, 1]))
    [W, H] = freqz([1,(2*a2-8) / (a2+2*a1+4),(a2-2*a1 + 4) / (a2+2*a1+4)], T*T*np.array([1, 2, 1]))
    plt.plot(W, np.angle(H))
    plt.title('Discrete')
    plt.show()
    # zeros = np.roots([1,(2*a2-8) / (a2+2*a1+4),(a2-2*a1 + 4) / (a2+2*a1+4)])
    # t = np.linspace(0, 2 * np.pi, 1000)
    # plt.scatter(np.cos(t), np.sin(t))
    # plt.scatter(np.real(zeros), np.imag(zeros))
    # v = plt.axis()
    # plt.axis([-2, 2, -2, 2])
    # plt.show()

# for i in range(10):
#     a1 = np.abs(np.random.randn())
#     a2 = 10*np.abs(np.random.randn())
#     plot_freqz(a1, a2)

INPUT_OMEGAS = np.array([np.pi / 4, np.pi / 3])
INPUT_PHIS = np.array([np.pi / 6, 2 * np.pi / 3])
epsilon = 0.3
omega1, omega2 = INPUT_OMEGAS
phi1, phi2 = INPUT_PHIS[0] + epsilon, np.pi - epsilon
r1, r2, a1, a2 = solve(omega1, omega2, phi1, phi2)
assert r1 > 0 and r2 > 0
[W, H] = freqz([1, a1, a2], [1, 2, 1])
plt.plot(W, np.angle(H))
print(([omega1, phi1], [omega2, phi2]))
plt.scatter([omega1, omega2], [phi1, phi2])
plt.show()