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

A_c = 1

f_1 = 500
T_1 = 1/f_1
omega_1 = 2*np.pi*f_1

f_c = 20000
T_c = 1/f_c
omega_c = 2*np.pi*f_c

def m(t):
    return 2*np.cos(2*omega_1*t) + 3*np.cos(3*omega_1*t)

def s(t):
    return A_c*m(t)*np.cos(omega_c*t)

f_s = 40*f_c
T_s = 1/f_s
t = np.arange(0,T_1,T_s)

plt.plot(t,s(t))
plt.plot(t,m(t))
plt.show()

S = sp.fft.fft(s(t))[35:45]
f = t*f_s/T_1
f = f[35:45]
for i in range(len(f)):
    print(f[i], S[i].real)

plt.stem(f,S.real/1600)
plt.show()