sdr_averaging_plot.py

#
# Assume white phase noise sampled at 1 MS/s
# visualize what happens when we decimate 10x six times to 1 S/s
#
import allantools as at
import numpy
import matplotlib.pyplot as plt
import scipy

# simulated data
f_sample =1e6
b0 = pow(2e-15,2) # desired PSD
num_pts=pow(2,25)
print "seconds", num_pts/f_sample
x6 = at.noise.white(num_points=num_pts, b0=b0,fs=f_sample)

def ddc10(x):
    # n=400, FIR order, default is 20x downsampling-factor
    return scipy.signal.decimate(x-numpy.mean(x), 10, ftype='fir')
    
x5 = ddc10(x6)

x4 = ddc10(x5)
x3 = ddc10(x4)
x2 = ddc10(x3)
x1 = ddc10(x2)
x0 = ddc10(x1)

def db_gain(x1,x2):
    return 10*numpy.log10( numpy.var(x1)/numpy.var(x2) )

print len(x6), numpy.std(x6)
print len(x5), numpy.std(x5), db_gain(x5,x6)
print len(x4), numpy.std(x4), db_gain(x4,x6)
print len(x3), numpy.std(x3), db_gain(x3,x6)
print len(x2), numpy.std(x2), db_gain(x2,x6)
print len(x1), numpy.std(x1), db_gain(x1,x6)
print len(x0), numpy.std(x0), db_gain(x0,x6)

#print len(x0)

# ADEVs
(tau6,dev6,err,n) = at.oadev(x6, rate=1e6)
(tau5,dev5,err,n) = at.oadev(x5, rate=1e5)
(tau4,dev4,err,n) = at.oadev(x4, rate=1e4)
(tau3,dev3,err,n) = at.oadev(x3, rate=1e3)
(tau2,dev2,err,n) = at.oadev(x2, rate=1e2)
(tau1,dev1,err,n) = at.oadev(x1, rate=1e1)
(tau0,dev0,err,n) = at.oadev(x0, rate=1e0)
print "ADEV done"

# PSDs
(f6, psd6) = at.noise.scipy_psd(x6, f_sample=1e6, nr_segments=4)
(f5, psd5) = at.noise.scipy_psd(x5, f_sample=1e5, nr_segments=4)
(f4, psd4) = at.noise.scipy_psd(x4, f_sample=1e4, nr_segments=4)
(f3, psd3) = at.noise.scipy_psd(x3, f_sample=1e3, nr_segments=4)
(f2, psd2) = at.noise.scipy_psd(x2, f_sample=1e2, nr_segments=4)
print "PSD done"

plt.figure()
plt.loglog(f6,psd6)
plt.loglog(f5,psd5)
plt.loglog(f4,psd4)
plt.loglog(f3,psd3)
plt.loglog(f2,psd2)


plt.loglog([1, 1e5], [b0, b0],'--')

plt.grid()

plt.figure()
plt.loglog(tau6,dev6,'o',label='1e6 S/s')
plt.loglog(tau4,dev4,'o',label='1e4 S/s')
plt.loglog(tau2,dev2,'o',label='1e2 S/s')
plt.loglog(tau0,dev0,'o',label='1e0 S/s')
plt.loglog(tau6, [2e-12/xx for xx in tau6],'--',label='2e-12/tau')
plt.loglog(tau6, [2e-13/xx for xx in tau6],'--',label='2e-13/tau')
plt.loglog(tau6, [2e-14/xx for xx in tau6],'--',label='2e-14/tau')
plt.loglog(tau6, [2e-15/xx for xx in tau6],'--',label='2e-15/tau')


#plt.loglog(tau2,dev2,'o',label='10 S/s, 0 IF')
#plt.loglog(tau1, [22e-13/xx for xx in tau1],'--',label='22e-13/tau')
#plt.loglog(tau3,dev3,'o',label='100 S/s, 0 IF')
#plt.loglog(tau1, [60e-13/xx for xx in tau1],'--',label='60e-13/tau')
#plt.loglog(tau4,dev4,'o',label='1000 S/s')
#plt.title('Ettus N210, 10MHz into ch1/ch2')
#plt.loglog(tau1, [7e-14/xx for xx in tau1],'--',label='7e-14/tau')
plt.grid()
plt.legend()
plt.figure()
#plt.plot(x1-numpy.mean(x1),label='1S/s')

plt.plot(x4-numpy.mean(x4),label='1000S/s')
plt.plot(x3-numpy.mean(x3),label='100S/s')
plt.plot(x2-numpy.mean(x2),label='10S/s')

plt.legend()

plt.figure()
plt.hist( x4-numpy.mean(x4),histtype='step',normed=True,label='1000S/s')
plt.hist( x3-numpy.mean(x3),histtype='step',normed=True,label='100S/s')
plt.hist( x2-numpy.mean(x2),histtype='step',normed=True,label='10S/s')
plt.hist( x1-numpy.mean(x1),histtype='step',normed=True,label='1S/s')
plt.show()