#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('pylab', 'inline')


# In[2]:


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


# Example continuous-time lowpass modulator - dsexample #3
# ========================================================
# 
# Demonstration of a low-pass, 3rd-order continuous-time delta sigma modulator, in a very similar fashion as as done in the **MATLAB Delta Sigma Toolbox**, albeit employing the Python port of the MATLAB package, **`deltasigma`**.
# 
# 
#  * The **Noise Transfer Function** (NTF) is synthesized for a **3th-order**, **low-pass** modulator.
#  * Time domain simulations are then performed to **extract the SNR**.
#  * The continuous-time NTF is realized and mapped to an **equivalent discrete-time structure**.
#  * The SNQR is again assessed.
#  * The modulator coefficients are then **scaled** and the final results are presented.

# ##General modulator parameters

# In[3]:


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


# ## Synthesis

# In[4]:


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 ""


# In[5]:


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


# In[6]:


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


# ##Time-domain NTF simulation

# In[7]:


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


# ##Signal to Noise Ratio

# If we have a two level quantizer, we run the prediction routine `predictSNR`, which employs the describing function method of Ardalan and Paulos.
# 
# Otherwise, we skip the prediction step and we move directly to time simulation of the NTF, to evaluate the SNR.

# In[8]:


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


# In[9]:


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);


# ## Realize the NTF and map it to an equivalent discrete-time structure
# 
# Performed in two steps:
# 
#  * We call `realizeNTF_ct` to realize an NTF with a continuous-time loop filter, obtaining ABCD matrix and DAC timings. 
# 
#  * We map the synthesized NTF to a discrete-time equivalent. The criterion for equivalence is set in `mapCtoD`: the sampled pulse response of the CT system must be identical to the impulse response of the DT system.

# In[10]:


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


# In[11]:


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)');


# In[12]:


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');


# In[13]:


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


# In[14]:


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='bottom');
figureMagic([-80, 0], 10, 1, [0, 80], 10, 1, None,'SQNR vs. Input Amplitude')
xlabel('Input Amplitude (dBFS)')
ylabel('SNR (dB)')
title('Continuous-Time Implementation');


# ## Dynamic range scaling

# In[15]:


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


# ##Conclusions: synthesis results
# 
# The ADC parameters were found to be:

# In[16]:


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]


# ### System version information

# In[17]:


#%install_ext http://raw.github.com/jrjohansson/version_information/master/version_information.py
get_ipython().run_line_magic('load_ext', 'version_information')
get_ipython().run_line_magic('reload_ext', 'version_information')

get_ipython().run_line_magic('version_information', 'numpy, scipy, matplotlib, deltasigma')