fs = 1e6
nbias = (2**12)/2

fg = 20e3

def load_samples(path):
    for line in open(path):
        f = line.split()
    
        if f[0] != 'DS':
            continue
        if f[-1] != 'DE':
            continue
    
        x = array(map(float, f[1:-1]))

        x_adc1 = x[::2]        
        x_adc2 = x[1::2]
        
        return x_adc1, x_adc2
    
x_adc1, x_adc2 = load_samples("test_20khz.log")

N = len(x_adc1)
assert len(x_adc1) == len(x_adc2)

print N

n_adc1 = arange(1, N*2, 2)
n_adc2 = arange(0, N*2, 2)

plot(n_adc1, x_adc1, 'bo')
plot(n_adc2, x_adc2, 'ro')
x1, x2, y1, y2 = axis()
axis((0, 50, y1, y2));

x = empty(N*2)
x[1::2] = x_adc1
x[::2] = x_adc2

n = arange(0, N*2)

plot(n, x)
grid()
x1, x2, y1, y2 = axis()
axis((0, 100, y1, y2));

from scipy.optimize import leastsq

def x_model(c, n):
    return c[0] * sin(2. * pi * c[1] * (n + c[2])) + c[3]

def x_error(c, x, n):
    return x_model(c, n) - x

c0 = [ (max(x) - min(x))/2., fg/(2.*fs), -60, mean(x) ]

c = leastsq(x_error, c0, args=(x, n))[0]

c_adc2 = leastsq(x_error, c0, args=(x_adc2, n_adc2))[0]
c_adc1 = leastsq(x_error, c0, args=(x_adc1, n_adc1))[0]

xm = x_model(c, n)
xm_adc1 = x_model(c_adc1, n)
xm_adc2 = x_model(c_adc2, n)

plot(n, xm_adc1 - xm, 'r-')
plot(n_adc1, x_adc1 - xm[1::2], 'ro')
plot(n, xm_adc2 - xm, 'b-')
plot(n_adc2, x_adc2 - xm[::2], 'bo')
x1, x2, y1, y2 = axis()
axis((0, 100, y1, y2))
ylabel("difference from model");

dA_adc1, df_adc1, dph_adc1, doff_adc1 = c_adc1 - c

print "dA_adc1   = ", dA_adc1
print "df_adc1   = ", df_adc1
print "dph_adc1  = ", dph_adc1
print "doff_adc1 = ", doff_adc1

dA_adc2, df_adc2, dph_adc2, doff_adc2 = c_adc2 - c

print "dA_adc2   = ", dA_adc2
print "df_adc2   = ", df_adc2
print "dph_adc2  = ", dph_adc2
print "doff_adc2 = ", doff_adc2

print "A = ", c[0]

psd = fft.rfft(x)
plot(abs(psd[1:100]), 'o-');

x_adc1, x_adc2 = load_samples("test_500khz.log")

x = empty(N*2)
x[::2] = x_adc1
x[1::2] = x_adc2

plot(x[:50], 'o-');

psd = fft.rfft(x)
plot(abs(psd[1:]), 'o-');