atan2_experiments.py

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

INPUT_OMEGAS = np.array([np.pi / 4, np.pi / 3])
INPUT_PHIS = np.array([np.pi / 6, 2 * np.pi / 3])

def solve(omega1, omega2, phi1, phi2):
    # Let's formulate the equations accordingly
    A = np.array([
        [-np.sin(phi1), 0, np.sin(omega1), 0],
        [-np.cos(phi1), 0, np.cos(omega1), 1],
        [0, -np.sin(phi2), np.sin(omega2), 0],
        [0, -np.cos(phi2), 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 evaluate_feasibility(INPUT_OMEGAS, PHI1, epsilon = 0.2):
    print()
    omega1, omega2 = INPUT_OMEGAS
    phi1, phi2 = PHI1 + epsilon, np.pi - epsilon
    r1, r2, a1, a2 = solve(omega1, omega2, phi1, phi2)
    if r1 > 0 and r2 > 0:
        print(f"r1 = {r1}, r2 = {r2})")
        print(f"Worked for epsilon = {epsilon}")
        # [W, H] = freqz([1, a, b], [0, 0, 1])
        # plt.plot(W, np.angle(H))
        # print(([omega1, phi1], [omega2, phi2]))
        # plt.scatter([omega1, omega2], [phi1, phi2])
        # plt.show()
        omega1, omega2 = INPUT_OMEGAS
        eval_h = lambda omega : np.angle(np.exp(2j * omega) + a1 * np.exp(1j * omega) + a2)
        print(f"Design phases: {phi1}, {phi2}")
        print(f"Obtained phases: {eval_h(omega1)}, {eval_h(omega2)}")
        phi1, phi2 = epsilon, np.pi - INPUT_PHIS[1] - epsilon
        r1, r2, b1, b2 = solve(omega1, omega2, phi1, phi2)
        if r1 > 0 and r2 > 0:
            print("Worked again!")
            eval_h = lambda omega : np.angle(np.exp(2j * omega) + b1 * np.exp(1j * omega) + b2)
            print(f"Denom Design phases: {phi1}, {phi2}")
            print(f"Denom Obtained phases: {eval_h(omega1)}, {eval_h(omega2)}")
            [W, H] = freqz([1, a1, a2], [1, b1, b2])
            plt.plot(W, np.angle(H))
            plt.scatter(INPUT_OMEGAS, INPUT_PHIS)
            plt.show()
        else:
            print(f"DENOM r1 = {r1}, r2 = {r2}")
            print(f"DENOM Failed for epsilon = {epsilon}")

    else:
        print(f"r1 = {r1}, r2 = {r2}")
        print(f"Failed for epsilon = {epsilon}")

for epsilon in np.arange(0.1, 1, 0.1):
    evaluate_feasibility(INPUT_OMEGAS, INPUT_PHIS[0], epsilon = epsilon)