Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
import numpy as np
import scipy
from scipy import signal
import pylab as pl
if __name__ == "__main__":
N = 1024*8
A = 2**14
An = 400
phase = np.pi/5
fs = 250E6
t = [(i/fs) for i in range(N)]
window = signal.flattop(N) #np.hamming(N)
# f_current = 300E6
for f_current in np.arange(50E6,500E6,10E6):
if f_current % 125E6 == 0:
continue
print("")
print("%g Hz"%(f_current))
noise1 = np.random.normal(0,1,N)
noise2 = np.random.normal(0,1,N)
ref1 = [A*np.sin(2*np.pi*f_current*i) for i in t]
ref1 = [(An*noise1[i] + ref1[i])*window[i] for i in range(N)]
comp = [.8*A*np.sin(2*np.pi*f_current*i + phase) for i in t]
comp = [(An*noise2[i] + comp[i])*window[i] for i in range(N)]
ref1_f = np.fft.rfft(ref1)
comp_f = np.fft.rfft(comp)
# pl.plot(ref1)
#pl.show()
#pl.subplot(1,2,1)
#pl.plot(abs(ref1_f))
#pl.plot(abs(comp_f))
#pl.subplot(1,2,2)
#pl.plot(np.angle(ref1_f))
#pl.plot(np.angle(comp_f))
f_div = fs/(N)
nyquist_region = int(f_current/(.5*fs))
f_mod = f_current % fs
if (nyquist_region % 2) == 1:
idx = N/2 - int((f_mod/f_div))
else:
idx = int((f_mod/f_div))
#print idx
print(abs(comp_f[idx])/abs(ref1_f[idx]))
print((np.angle(comp_f[idx])-np.angle(ref1_f[idx])))
# pl.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.