"""
ECE 09351 - Digital Signal Processing
Lab 1 Question 3

Translated from MATLAB by Aidan Sharpe
Last Modified: January 17, 2025
"""

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt


def cft(g, f):
    result = np.zeros(len(f), dtype=complex)

    for i, ff in enumerate(f):
        result[i] = complex_quad(lambda t: g(t)*np.exp(-2j*np.pi*ff*t), -5, 5)

    return result

def complex_quad(g, a, b):
    t = np.linspace(a, b, 2501)
    x = g(t)
    return sp.integrate.simpson(y=x, x=t)


def main():
    g = lambda t: np.sinc(0.5*t/np.pi)

    f_s = 1
    T_s = 1/f_s
    t_0 = 1000
    t = np.arange(-t_0, t_0, T_s)
    f = np.linspace(-f_s/2, f_s/2, len(t), endpoint=False)
    omega = 2*np.pi*f

    G = np.fft.fftshift(np.fft.fft(g(t)) * np.exp(-2j*np.pi*f*t_0) /f_s)
    G_analytic = 2*np.pi*(np.heaviside(omega+0.5, 1) - np.heaviside(omega-0.5, 1))

    plt.plot(f, G_analytic, color="red")
    plt.scatter(f, G)
    plt.show()
    


if __name__ == "__main__":
    main()