%pylab inline

from __future__ import division
from deltasigma import *
import warnings
from scipy.signal import ss2zpk
warnings.filterwarnings('ignore')

order = 3
osr = 32
nlev = 2
f0 = 0.
Hinf = 1.5
tdac = [0, 1]
form = 'FB'
M = nlev - 1
#generic options
plotsize = (15, 8) #inches

print '\t\t\t %dth-Order Continuous-Time Lowpass Example\n' % order
print 'Doing NTF synthesis...',
ntf0 = synthesizeNTF(order, osr, 2, Hinf, f0) # 2: Optimized zero placement
print 'done.'
print "\nSynthesized a %d-order NTF, with roots:\n" % order
print " Zeros:\t\t\t Poles:"
for z, p in zip(ntf0[0], ntf0[1]):
    print "(%f, %fj)\t(%f, %fj)" % (np.real(z), np.imag(z), np.real(p), np.imag(p))
print ""

figure(figsize=plotsize)
plotPZ(ntf0, showlist=True)
changeFig(10, 1.5, 7)

figure(figsize=plotsize)
DocumentNTF(ntf0, osr, f0)

print 'Doing time-domain simulations...'
figure(figsize=plotsize)
PlotExampleSpectrum(ntf0, M, osr, f0)
title('Example Spectrum');

if nlev == 2:
    snr_pred, amp_pred, k0, k1, se = predictSNR(ntf0, osr)
snr, amp = simulateSNR(ntf0, osr, None, f0, nlev)

figure(figsize=plotsize)
if nlev == 2:
    plot(amp_pred, snr_pred, '-', label='Predicted')
    hold(True)
plot(amp, snr,'o-.g', label='simulated')
#figureMagic(np.array([- 100,0]).reshape(1,-1),10,2,np.array([0,100]).reshape(1,-1),10,2,np.array([7,3]).reshape(1,-1),'Discrete-Time Simulation')
xlabel('Input Level (dBFS)')
ylabel('SQNR (dB)')
peak_snr, peak_amp = peakSNR(snr, amp)
msg = 'peak SQNR = %4.1fdB  \n@ amp = %4.1fdB  ' % (peak_snr, peak_amp)
text(peak_amp-10, peak_snr, msg, horizontalalignment='right', verticalalignment='center');
msg = 'OSR = %d ' % osr
text(-2, 5, msg, horizontalalignment='right');
hold(False)
figureMagic([-80, 0], 10, None, [0, 80], 10, None, [12, 6], 'Time-Domain Simulations')
legend(loc=2);

print 'Mapping  to continuous-time... ',
ABCDc, tdac2 = realizeNTF_ct(ntf0, form, tdac)
Ac, Bc, Cc, Dc = partitionABCD(ABCDc)
sys_c = []
for i in range(Bc.shape[1]):
    sys_c.append(ss2zpk(Ac, Bc, Cc, Dc, input=i))
print 'done.'
print 'ABCD matrix:'
print ABCDc
print "DAC timing (tdac2):"
print tdac2

figure(figsize=plotsize)
n_imp = 10
y = -impL1(ntf0, n_imp)
lollipop(np.arange(n_imp + 1), y)
hold(True)
grid(True)
dt = 1./16
tppulse = np.vstack((np.zeros((1, 2)), tdac2[1:, :])).tolist()
yy = -pulse(sys_c, tppulse, dt, n_imp).squeeze()
t = np.linspace(0, n_imp + dt, 10/dt + 1)
plot(t, yy, 'g', label='continuous-time')
legend()
title('Loop filter pulse/impulse responses (negated)');

sys_d, Gp = mapCtoD(ABCDc, tdac2)
ABCD = np.vstack((
                  np.hstack((sys_d[0], sys_d[1])),
                  np.hstack((sys_d[2], sys_d[3]))
                ))
ntf, G = calculateTF(ABCD)
ntf = cancelPZ(ntf)
figure(figsize=plotsize)
DocumentNTF(ntf0, osr, f0)
subplot(121)
plotPZ(ntf, 'c', 10)
hold(True)
plotPZ(ntf0, 'k')
hold(False)

L0 = sys_c[0]
f = linspace(0, 0.5)
G = evalTFP(L0, ntf, f)
subplot(122)
hold(True)
plot(f, dbv(G), 'm')
hold(False)
title('NTF and STF');

print 'Re-evaluating the SNR... ',
snrR, ampR = simulateSNR(ABCD, osr, None, f0, nlev)
print 'done'

figure(figsize=plotsize)
plot(ampR, snrR, 'o-')
hold(True)
plot(amp, snr, '-')
peak_snrR, peak_ampR = peakSNR(snrR, ampR)
msg = 'Peak SNR %.1fdB at amp = %-.1fdB' % (peak_snrR, peak_ampR)
text(peak_ampR - 10, peak_snrR, msg, horizontalalignment='right', verticalalignment='center');
figureMagic([-100, 0], 10, 1, [0, 100], 10, 1, None,'SQNR vs. Input Amplitude')
xlabel('Input Amplitude (dBFS)')
ylabel('SNR (dB)')
title('Continuous-Time Implementation');

print 'Doing dynamic range scaling... ',
ABCDs, umax, S = scaleABCD(ABCD, nlev, f0, 1, None, None, 10000.0)
S = S[:order, :order]
Sinv = inv(S)
Acs = np.dot(np.dot(S, Ac), Sinv)
Bcs = np.dot(S, Bc)
Ccs = np.dot(Cc, Sinv)
ABCDcs = np.vstack((np.hstack((Acs, Bcs)),
                    np.hstack((Ccs, Dc))
                  ))
sys_cs = (Acs, Bcs, Ccs, Dc)
print 'Done.'
print "During scaling, umax was found to be %g"% umax
print "Scaled ABCD matrix:"
print ABCDcs

adc = {
    'order':order,
    'osr':osr,
    'M':M,
    'f0':f0,
    'ntf':ntf,
    'ABCD':ABCD,
    'umax':umax,
    'peak_snr':peak_snr,
    'form':form,
    'ABCDc':ABCDc,
    'tdac':tdac,
    'tdac2':tdac2,
    'L0':L0,
    'sys_c':sys_c,
    'ABCDcs':ABCDcs,
    'sys_cs':sys_cs
}
for k in sort(adc.keys()):
    print "%s:" % k,
    if str(adc[k]).count('\n'):
        print ""
    print adc[k]

#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
%load_ext version_information
%reload_ext version_information

%version_information numpy, scipy, matplotlib, deltasigma