#' % Continuous-Time Sigma Delta Converters
#' % Wolfgang Scherr
#' % 14th February 2021

#' (c) Carinthia University of Applied Sciences<br />
#' All rights reserved. For educational purposes only.

#' SystemC-AMS models created with [Coside](https://www.coseda-tech.com/coside-overview).

#+ echo=False
#-------------------------------------------------------------------
import shutil

import os

import math

import subprocess

import schemdraw as sd
from schemdraw import dsp

import numpy as np
from numpy.fft import rfft

import matplotlib
import matplotlib.pyplot as plt

import sympy as sym
from sympy import symbols

def bode_plot(H,start,stop,title):
    # Frequency vector
    freq = np.logspace(start,stop) # 10^n Hz
    s_num = 2j * np.pi * freq
    f = sym.lambdify(s, H, modules='numpy')
    Hs_num = np.abs(f(s_num))

    # Plotting the Bode diagram
    fig, ax = plt.subplots()
    ax.loglog(freq, Hs_num )
    ax.set(xlabel='freq (Hz)', ylabel='magnitude |H|', title=title)
    ax.grid()
    plt.show()


def sdadc_svg(filt='$\int$',svg_file='figures/sdadc_ideal.svg'):
    d = sd.Drawing()

    # forward path X -> Y
    d.add(dsp.Arrow(l=d.unit/2, label='X'))
    smx = d.add(dsp.SumSigma)
    d.add(dsp.Arrow(xy=smx.E, l=d.unit/2))
    intg = d.add(dsp.Box(w=2, h=2, label=filt, anchor='W'))
    d.add(dsp.Arrow('right', l=d.unit/2, xy=intg.E))
    smn = d.add(dsp.SumSigma)
    d.add(dsp.Line(l=d.unit/3))
    fb = d.add(dsp.Dot)
    d.add(dsp.Arrow(l=d.unit/3, label='Y'))

    # feedback path Y -> sum X
    d.add(dsp.Line('down', xy=fb.center, l=d.unit/2))
    d.add(dsp.Line('left', tox=smx.S))
    d.add(dsp.Line('up', toy=smx.S))
    d.add(dsp.Arrowhead(botlabel='-'))

    # noise input N
    d.add(dsp.Arrow('up', xy=smn.N, l=d.unit/2, reverse=True, label='N'))

    #d.draw()
    d.save(svg_file)


# see e.g.: https://stackoverflow.com/questions/27306474/plotting-fft-from-a-wav-file-using-python
def fft_plot(t,sig,xlabel,ylabel,title,log=0):
    Fs = len(t) / max(t)                   # Sampling frequency
    nfft = 2**np.floor(np.log2(len(t)))    # Use next highest power of 2 greater
                                           # than or equal to length(x) to calculate FFT.
    fftx = rfft(sig,int(nfft))             # Take FFT, padding with zeros so that length(fftx) is equal to nfft
    NumUniquePts = int(np.floor((nfft+1)/2.0)) # Calculate the number of unique points
    fftx = fftx[0:NumUniquePts]                # FFT is symmetric, throw away second half
    mx = abs(fftx) / float(nfft)               # Take the magnitude of FFT of x and 
                                               # scale the FFT so that it is not a function of the length of x
    # multiply by two (see technical document for details)
    # odd nfft excludes Nyquist point
    if nfft % 2 > 0:                                   # odd number of points FFT
            mx[1:len(mx)] = mx[1:len(mx)] * 2
    else:                                              # even number of points FFT
            mx[1:len(mx) -1] = mx[1:len(mx) - 1] * 2
    f = np.arange(0, NumUniquePts, 1.0) * (Fs / nfft)  # create frequency array

    fig, ax = plt.subplots()
    ax.loglog(f, mx)

    ax.set(xlabel=xlabel, ylabel=ylabel, title=title)
    ax.grid()

    plt.show()


def ac_plot(f,re,im,title):
    size = len(f)

    y = []
    
    for x in range(0, size):
        y.append(math.sqrt(re[x]**2 + im[x]**2))

    fig, ax = plt.subplots()
    ax.loglog(f, y)

    ax.set(xlabel='freq [Hz]', ylabel='|magnitude|', title=title)
    ax.grid()

    plt.show()


def tran_plot(t,sig,tmin,tmax,ylabel,title):
    size = len(t)

    x = []
    y = []
    
    for i in range(0, size):
        if (t[i]>=tmin) and (t[i]<=tmax):
            x.append(t[i])
            y.append(sig[i])

    fig, ax = plt.subplots()
    ax.plot(x, y)

    ax.set(xlabel='time [s]', ylabel=ylabel, title=title)
    ax.grid()

    plt.show()


# Laplace variable "s":
s = symbols('s')
#-------------------------------------------------------------------


#' # Data Converter Basics

#' This document does not elaborate on basics of signal processing and
#' sampling or history of A/D converters. If you are interested in that,
#' there is plenty of information on the web, e.g from
#' [Analog Devices here](https://www.analog.com/media/en/training-seminars/tutorials/MT-022.pdf).

#' # Ideal continuous-time sigma-delta converter

#' For an ideal c.t. sigma-delta converter, an integrator is used.
#' An equivalent SystemC-AMS (IEEE 1666.1) model using LSF (linear signal flow) looks like this:

#+ echo=False, results = 'hidden'
#-------------------------------------------------------------------
shutil.copyfile('../sigma_delta/doc/lsf_simple_schematic.svg','figures/lsf_simple_schematic.svg');
#-------------------------------------------------------------------
#' <img  src="figures/lsf_simple_schematic.svg" />

#' This leads to this signal transfer characteristics:

#' $Hs(s) = \frac{Y(s)}{X(s)} = \frac{\frac{1}{s}}{1 + \frac{1}{s}} =
#'                             \frac{1}{s + 1}$

#' And to this noise transfer characteristics:

#' $Hn(s) = \frac{Y(s)}{N(s)} = \frac{1}{1 + \frac{1}{s}}$

#' Which can be nicely checked by frequency-domain simulations of the shown model:

#+ echo=False
#-------------------------------------------------------------------
#change to simulation directory
cpath=os.getcwd();
os.chdir('../DEBUG/sigma_delta');
#create command line with parameter of the SystemC simulation
cmd = 'lsf_simple_simple_tb.exe --x_ac=1.0 --n_ac=0.0 ';
#start systemc simulation
print ('Running : ' + cmd);
subprocess.call(cmd,stdout=subprocess.PIPE,stderr=subprocess.PIPE);
#return to script directory
os.chdir(cpath);
dat = []
dat = np.loadtxt('../DEBUG/sigma_delta/lsf_simple_ac.dat', comments='%', unpack=True)
ac_plot(dat[0],dat[5],dat[6],'AC: 1.O. Sigma-Delta signal path H=Y/X')
#-------------------------------------------------------------------

#+ echo=False
#-------------------------------------------------------------------
#change to simulation directory
cpath=os.getcwd();
os.chdir('../DEBUG/sigma_delta');
#create command line with parameter of the SystemC simulation
cmd = 'lsf_simple_simple_tb.exe --x_ac=0.0 --n_ac=1.0 ';
#start systemc simulation
print ('Running : ' + cmd);
subprocess.call(cmd,stdout=subprocess.PIPE,stderr=subprocess.PIPE);
#return to script directory
os.chdir(cpath);
dat = []
dat = np.loadtxt('../DEBUG/sigma_delta/lsf_simple_ac.dat', comments='%', unpack=True)
ac_plot(dat[0],dat[5],dat[6],'AC: 1.O. Sigma-Delta noise path H=Y/N')
#-------------------------------------------------------------------

#' This circuit has some problems. First of all, there is no ideal integrator
#' in the analog world. In principle, they any implementation is more an
#' active low-pass filter with some low corner frequency.

#' Secondly, the ideal integrator does not really provide a proper bandwidth
#' to be used. For example: if we want a SD-ADC with 256MHz sampling rate
#' to convert a 1MHz signal, we need more bandwidth for the signal path.

#+ echo=False

#' # A bit more realistic sigma-delta converter

#' Now we use a low-pass filter instead of an ideal integrator.

#+ echo=False, results = 'hidden'
#-------------------------------------------------------------------
shutil.copyfile('../sigma_delta/doc/lsf_realistic_schematic.svg','figures/lsf_realistic_schematic.svg');
#-------------------------------------------------------------------
#' <img  src="figures/lsf_realistic_schematic.svg" />

#' The first order lowpass filter function is:
#' $Hlp(s) = \frac{1}{1 + \tau s}$.

#' We set the corner frequency $f_c$ to about 2MHz  ($\tau = \frac{1}{2\pi f_c}$).
#' The transfer function of this low-pass filter alone will look as follows:

#+ echo=False
#-------------------------------------------------------------------

# signal transfer function
tau = 1.0/(2.0*math.pi*2e6)
print('Tau = ' + str(tau))
Hlp = 1 / (1+tau*s)
print('Hlp(s) = ' + str(Hlp))

bode_plot(Hlp,3,9,'Bode diagram of 1.O. LP filter')
#-------------------------------------------------------------------

#' In the closed loop setup of the SD-ADC we will now get for the signal transfer function:

#' $Hs(s) = \frac{Y(s)}{X(s)} = \frac{\frac{1}{1+\tau s}}{1 + \frac{1}{1+\tau s}} =
#'                             \frac{1}{\tau s + 2}$

#' Looking at the noise transfer function we have now:

#' $Hn(s) = \frac{Y(s)}{N(s)} = \frac{1}{1 + \frac{1}{1+\tau s}}$

#' Which we again check by a frequency domain simulation:

#+ echo=False
#-------------------------------------------------------------------
#change to simulation directory
cpath=os.getcwd();
os.chdir('../DEBUG/sigma_delta');
#create command line with parameter of the SystemC simulation
cmd = 'lsf_realistic_simple_tb.exe --x_ac=1.0 --n_ac=0.0 ';
#start systemc simulation
print ('Running : ' + cmd);
subprocess.call(cmd,stdout=subprocess.PIPE,stderr=subprocess.PIPE);
#return to script directory
os.chdir(cpath);
dat = []
dat = np.loadtxt('../DEBUG/sigma_delta/lsf_realistic_ac.dat', comments='%', unpack=True)
ac_plot(dat[0],dat[5],dat[6],'AC: 1.O. Sigma-Delta signal path (LP-Filter) H=Y/X')
#-------------------------------------------------------------------

#+ echo=False
#-------------------------------------------------------------------
#change to simulation directory
cpath=os.getcwd();
os.chdir('../DEBUG/sigma_delta');
#create command line with parameter of the SystemC simulation
cmd = 'lsf_realistic_simple_tb.exe --x_ac=0.0 --n_ac=1.0 ';
#start systemc simulation
print ('Running : ' + cmd);
subprocess.call(cmd,stdout=subprocess.PIPE,stderr=subprocess.PIPE);
#return to script directory
os.chdir(cpath);
dat = []
dat = np.loadtxt('../DEBUG/sigma_delta/lsf_realistic_ac.dat', comments='%', unpack=True)
ac_plot(dat[0],dat[5],dat[6],'AC: 1.O. Sigma-Delta noise path (LP-Filter) H=Y/N')
#-------------------------------------------------------------------


#' This does not really like good noise shaping. So we need to do a bit more,
#' which is introducing an active LP filter with some gain, looking at the filter alone we get now:

#+ echo=False
#-------------------------------------------------------------------

# signal transfer function
tau = 1.0/(2.0*math.pi*2e6)
print('Tau = ' + str(tau))
Hlp = 100 / (1+tau*s)
print('Hlp(s) = ' + str(Hlp))

bode_plot(Hlp,3,9,'Bode diagram of 1.O. LP filter with gain=100')
#-------------------------------------------------------------------

#' What is now the effect in the closed-loop simulation of the converter?
#' Let's check the frequency-domain simulation again:

#+ echo=False
#-------------------------------------------------------------------
#change to simulation directory
cpath=os.getcwd();
os.chdir('../DEBUG/sigma_delta');
#create command line with parameter of the SystemC simulation
cmd = 'lsf_realistic_simple_tb.exe --x_ac=1.0 --n_ac=0.0 --gain=100.0';
#start systemc simulation
print ('Running : ' + cmd);
subprocess.call(cmd,stdout=subprocess.PIPE,stderr=subprocess.PIPE);
#return to script directory
os.chdir(cpath);
dat = []
dat = np.loadtxt('../DEBUG/sigma_delta/lsf_realistic_ac.dat', comments='%', unpack=True)
ac_plot(dat[0],dat[5],dat[6],'AC: 1.O. Sigma-Delta signal path (act. LP-Filter) H=Y/X')
#-------------------------------------------------------------------

#+ echo=False
#-------------------------------------------------------------------
#change to simulation directory
cpath=os.getcwd();
os.chdir('../DEBUG/sigma_delta');
#create command line with parameter of the SystemC simulation
cmd = 'lsf_realistic_simple_tb.exe --x_ac=0.0 --n_ac=1.0 --gain=100.0';
#start systemc simulation
print ('Running : ' + cmd);
subprocess.call(cmd,stdout=subprocess.PIPE,stderr=subprocess.PIPE);
#return to script directory
os.chdir(cpath);
dat = []
dat = np.loadtxt('../DEBUG/sigma_delta/lsf_realistic_ac.dat', comments='%', unpack=True)
ac_plot(dat[0],dat[5],dat[6],'AC: 1.O. Sigma-Delta noise path (act. LP-Filter) H=Y/N')
#-------------------------------------------------------------------

#' Not so bad, right? We get a better signal bandwidth and shift the noise upwards,
#' so really nice noise-shaping, even for a very simple 1.O. ct. SD-ADC...

#' # Ideal implementation model

#' But how does this look more on electrical level? Let's do a basic
#' circuit model for our SD-ADC using SystemC-AMS (IEEE 1666.1) ELN (electrical linear network).

#' For that we define a simple 1.O. operational amplifier. The advantage of using
#' SystemC-AMS for the model is, that it runs way faster than
#' using SPICE, despite I can still use a SPICE-like entry many are familar with.
#' Basically, it is nearly as fast as the LSF model we used before (we just have a bit larger network to solve)...

#+ echo=False, results = 'hidden'
#-------------------------------------------------------------------
shutil.copyfile('../basic_blocks/doc/op_amp_model_schematic.svg','figures/op_amp_model_schematic.svg');
#-------------------------------------------------------------------
#' <img  src="figures/op_amp_model_schematic.svg" />

#' Its frequency response simulated looks like this:

#+ echo=False
#-------------------------------------------------------------------
dat = np.loadtxt('../DEBUG/basic_blocks/i_opa_top_ac.dat', comments='%', unpack=True)
ac_plot(dat[0],dat[3],dat[4],'AC: Opamp 1.O, fc=2.0e6, gain=100.0')
#-------------------------------------------------------------------

#' The operational amplifier is included in the toplevel SD-ADC model. The output
#' of the SD-ADC is connected to two 1.O. CIC with a frequency response corresponding
#' to R=64, but we don't do the actual decimation to show the residual bins which
#' would be folded back into the signal band.

#' For decimation, I used a simple averaging filter. More about decimation can be found on
#' [one of my other modelling pages here](../filter/average_filter.html).

#+ echo=False, results = 'hidden'
#-------------------------------------------------------------------
shutil.copyfile('../sigma_delta/doc/sdadc_1o_schematic.svg','figures/sdadc_1o_schematic.svg');
#-------------------------------------------------------------------
#' <img  src="figures/sdadc_1o_schematic.svg" />
#' <!-- img  src="figures/sdadc_1o_schematic.svg" style="margin:0px -50% 0px 0px; clip: rect(0px, 0px, 0px, -50%);" -->

#' Let's look at the result of a transient simulation of this model. For the
#' FFT, I simulated 5ms, which took just a few seconds on my notebook. Here
#' I plot just a few periods of the input signal. But be aware, the sample
#' rate of the ADC is 256MHZ, so a lots of points to calculate and post-process.

#+ echo=False
#-------------------------------------------------------------------
dat = np.loadtxt('../DEBUG/sigma_delta/i_sdadc_1o_tran.dat', comments='%', unpack=True)
tran_plot(dat[0],dat[1],1e-6,3e-6,'vin','TRAN: SD-ADC input voltage')
tran_plot(dat[0],dat[2],1e-6,3e-6,'vlp','TRAN: SD-ADC integrator voltage')
tran_plot(dat[0],dat[3],1e-6,3e-6,'dout','TRAN: SD-ADC digital bitstream')
tran_plot(dat[0],dat[4],1e-6,3e-6,'dout_dec','TRAN: output 1.O. CIC filtered')
tran_plot(dat[0],dat[5],1e-6,3e-6,'dout_dec2','TRAN: output 2.O. CIC filtered')
#-------------------------------------------------------------------

#' The digital output stream can be better analysed using an FFT transform.
#' This is before any decimation:

#+ echo=False
#-------------------------------------------------------------------
fft_plot(dat[0],dat[3],'output freqency [Hz]','dout','TRAN: SD-ADC 1.O, 256MHz, bitstream')
#-------------------------------------------------------------------

#' This is the FFT transform after the 1.O CIC filter:

#+ echo=False
#-------------------------------------------------------------------
fft_plot(dat[0],dat[4],'output freqency [Hz]','dout_dec','TRAN: 1.O. CIC, filtered')
#-------------------------------------------------------------------

#' And here the FFT transform after the second 1.O CIC filter:

#+ echo=False
#-------------------------------------------------------------------
fft_plot(dat[0],dat[5],'output freqency [Hz]','dout_dec','TRAN: 2.O. CIC, filtered')
#-------------------------------------------------------------------

#' Finally, the FFT transform after both CIC filter and decimation by 64
#' (which means the new output sampling frequency is 4MHz):

#+ echo=False
#-------------------------------------------------------------------
dat = np.loadtxt('../DEBUG/sigma_delta/sdadc_1o_dec.dat', comments='%', unpack=True)
fft_plot(dat[0],dat[1],'output freqency [Hz]','dout_dec','TRAN: 2.O. CIC, R=64, decimated')
#-------------------------------------------------------------------

#' Please note that I didn't use any window function for the FFT, this is of course
#' not optimal. But this report is for the purpose of demonstration, not more. You
#' can find plenty of serious and much more detailed information on post processing of
#' such converters in books, in papers and on the web.