In [1]:

    

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import scipy.signal
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()


    

In [2]:

    

obs_rate = 100

def load_file(filename):
    phase = np.fromfile(filename, dtype = 'float32').astype('double')
    start = np.datetime64(filename.strip('.f32').split('_')[-1])
    timestamps = start + np.timedelta64(1, 'ns') * np.arange(phase.shape[0]) * 1e9 / obs_rate
    return xr.Dataset({'phase' : ('time', np.unwrap(phase))}, coords = {'time' : timestamps})


    

In [3]:

    

data = load_file('phase_bpsk_2019-11-20T20:40:13.986230.f32').sel(time = slice('2019-11-20T20:41', '2019-11-23'))


    

In [4]:

    

freq_corr = -1.5
freq_aliasing = 1
phase_corrected = ((data['phase'].diff('time')+freq_aliasing) % (2*np.pi) -freq_aliasing+freq_corr).cumsum('time')
data['phase'][1:] = phase_corrected
data['phase'][0] = 0


    

In [5]:

    

plt.figure(figsize = (12,6), facecolor = 'w')
(data['phase'][::100]/(2*np.pi)).differentiate('time', datetime_unit = 's').plot()
plt.title('BPSK beacon frequency offset')
plt.xlabel('UTC time')
plt.ylabel('Frequency (Hz)');


    

In [6]:

    

(data['phase'].coords['time'][-1] - phase_corrected.coords['time'][0]).astype('float')*1e-9


    

    
Out[6]:





<xarray.DataArray 'time' ()>
array(148986.92)

In [7]:

    

f_obs = 10489.8e6

def adev(series, skip, freq = 10e9, overlapping = False):
    x = series.values/(2*np.pi*freq)
    
    tau = skip / obs_rate
    if overlapping:
        y = x[:-2*skip] - 2*x[skip:-skip] + x[2*skip:]
    else:
        z = x[:x.size//skip*skip].reshape((-1,skip))[:,0]
        y = z[:-2] - 2*z[1:-1] + z[2:]
    return np.sqrt(0.5/tau**2*np.average(y**2))

def get_skips(n):
    if n <= 0:
        return np.array([], dtype = 'int')
    a = int(np.log10(n))
    step = max(10**(a-2), 1)
    return np.concatenate((get_skips(10**a - 1) , np.arange(10**a, n+1, step)))

def compute_adev(data, overlapping = False):
    skips = get_skips(data.coords['time'].size//2)
    taus = skips / obs_rate
    adevs = [adev(data['phase'], skip, f_obs, overlapping) for skip in skips]
    return xr.Dataset({'adev' : ('tau', adevs)}, coords = {'tau' : taus})


    

In [8]:

    

#adevs = compute_adev(data, overlapping = True)
#adevs.to_netcdf('adevs_qo100_vectron.nc')


    

In [9]:

    

adevs = xr.open_dataset('adevs_qo100_vectron.nc')


    

In [10]:

    

def plot_adev(a, label):
    plt.loglog(a.coords['tau'], a['adev'], label = f'{label}')


    

In [11]:

    

plt.figure(figsize = (12,6), facecolor = 'w')
plot_adev(adevs, 'QO-100 BPSK -- Vectron MD-011 (10.5GHz)')
plt.xlabel('$\\tau$ (s)')
plt.ylabel('$\\sigma(\\tau)$')
plt.legend()
plt.grid(which = 'both')
plt.title('Allan deviation');