#' % Playing with diode modulators...
#' % Wolfgang Scherr
#' % 13th March 2021

#' All rights reserved. For educational purposes only.

#+ echo=False
#-------------------------------------------------------------------
import math
import numpy as np
import matplotlib.pyplot as plt
from lcapy import Circuit
#-------------------------------------------------------------------


#+ echo=False
#-------------------------------------------------------------------
# see: https://kluedo.ub.uni-kl.de/frontdoor/deliver/index/docId/4293/file/exact_fft_measurements.pdf
# or:  https://stackoverflow.com/questions/15382076/plotting-power-spectrum-in-python/20171984
def fft_new(t,sig,xl,yl,title='',freq=-1):
    N = len(t)                                  # no. of samples

    two_peak = np.abs(np.fft.fft(sig)) / N      # sample scaling
    #two_rms = two_peak/np.sqrt(2)
    #two_pwr = two_rms**2
    #two_dbm = 10*np.log10(two_pwr/(50*0.001))
    ts = t[2]-t[1]                              # sampling timestep
    freqs = np.fft.fftfreq(sig.size, ts)
    idx = np.argsort(freqs)
    size = int(max(idx)/2)

    one_peak = two_peak[0:size]                 # dual sided -> single sided
    one_peak[1:size] = one_peak[1:size]*2
    one_rms = one_peak/np.sqrt(2)               # Vrms based on sine signals
    one_pwr = one_rms**2                        # Vrms -> Power
    one_dbm = 10*np.log10(one_pwr/(50*0.001))   # logarithmic conv.
    one_f = freqs[0:size]

    if title!='':
        fig, ax = plt.subplots()                    # plot FFT
        ax.plot(one_f, one_dbm)
        ax.set(xlabel='freq [Hz]', ylabel='pwr [dBm]', title=title)
        ax.grid()
        if xl!=[]:
            plt.xlim(xl)
        if yl!=[]:
            plt.ylim(yl) 
        plt.show()

    if freq==-1:
        return max(one_dbm)
    else:
        fmax = 0.5 / ts         # 0...fs/2 is 0...size ->  f  is   ?    
        fidx = int(freq*size/fmax+0.5)
        return one_dbm[fidx]
#-------------------------------------------------------------------


#+ echo=False
#-------------------------------------------------------------------
def ac_plot(f,re1,im1,l1,re2,im2,l2,title):
    size = len(f)

    y1 = []
    y2 = []
    
    for x in range(0, size):
        y1.append(20*math.log10(math.sqrt(re1[x]**2 + im1[x]**2)+1e-6))
        y2.append(20*math.log10(math.sqrt(re2[x]**2 + im2[x]**2)+1e-6))

    fig, ax = plt.subplots()
    ax.plot(f, y1, label=l1)
    ax.plot(f, y2, label=l1)

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

    plt.show()
#-------------------------------------------------------------------


#+ echo=False
#-------------------------------------------------------------------
def xy_plot(sigx,sigy,xmin,xmax,xlabel,ylabel,title):
    size = len(sigx)

    x = []
    y = []
    for i in range(0, size):
        if (sigx[i]>=xmin) and (sigx[i]<=xmax):
            x.append(sigx[i])
            y.append(sigy[i])

    fig, ax = plt.subplots()
    ax.plot(x, y, label='XY plot')
    ax.set(xlabel=xlabel, ylabel=ylabel, title=title)
    ax.grid()
    #ax.legend()
    plt.show()
#-------------------------------------------------------------------


#+ echo=False
#-------------------------------------------------------------------
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()
#-------------------------------------------------------------------


#' #Introduction

#' One day I stumbled over this web link 
#' [about modulator circuitries](https://rickettslab.org/bits2waves/design/mixer-discrete/mixer-discrete-theory/)
#' (external link).

#' Instead of playing games, why not once playing with circuit simulators?

#' So I thought it is fun to try this examples in SystemC-AMS ELN using the
#' piece-wise linear sources already available in the actual proof of concept library.

#' First, I build some basic blocks, afterwards I use them for the modulators.
#' Have fun trying them yourself!


#' #Elements for the circuits with basic tests

#' ###Low pass and high pass filter models

#' SystemC-AMS source code:

#' LC lowpass: [H](../elements/lc_lowpass.h) - 
#' [CPP](../elements/lc_lowpass.cpp) - 
#' [(TB) CPP](../elements/lc_lowpass_tb.cpp)

#+ echo=False
#-------------------------------------------------------------------
# check out: http://lcapy.elec.canterbury.ac.nz/schematics.html#introduction
cct = Circuit("""
 L1 inp inner; right=2
 C1 inner _gnd; down
 W _gnd __gnd; down=0.1, rground
 L2 inner outp; right=2""")
cct.draw()
#-------------------------------------------------------------------

#' </BR>
#' LC highpass: [H](../elements/lc_highpass.h) - 
#' [CPP](../elements/lc_highpass.cpp) - 
#' [(TB) CPP](../elements/lc_highpass_tb.cpp)

#+ echo=False
#-------------------------------------------------------------------
cct = Circuit("""
 C1 inp inner; right=2
 L1 inner _gnd; down
 W _gnd __gnd; down=0.1, rground
 C2 inner outp; right=2""")
cct.draw()
#-------------------------------------------------------------------

#' </BR>
#' AC simulation result:

#+ echo=False
#-------------------------------------------------------------------
dat1 = []
dat1 = np.loadtxt('../lp_ac.dat', comments='%', unpack=True)
dat2 = []
dat2 = np.loadtxt('../hp_ac.dat', comments='%', unpack=True)
ac_plot(dat1[0],dat1[3],dat1[4],'lp',dat2[3],dat2[4],'hp','AC: LC filter')
#-------------------------------------------------------------------

#' ###Diode model

#' SystemC-AMS source code:

#' Diode: [H](../elements/diode.h) - 
#' [CPP](../elements/diode.cpp) - 
#' [(TB) CPP](../elements/diode_tb.cpp)

#+ echo=False
#-------------------------------------------------------------------
cct = Circuit("""
   W anode _P; right=2
   E1 _P _N anode cathode IPWL; down=1, l=I_P_W_L
   W _N cathode; right=2
 """)
cct.draw()
#-------------------------------------------------------------------

#' </BR>
#' DC plot generated from a transient simulation using a voltage ramp:

#+ echo=False
#-------------------------------------------------------------------
dat3 = []
dat3 = np.loadtxt('../diode_tran.tab', comments='%', unpack=True)

xy_plot(dat3[1],-dat3[2],-21,2,'diode voltage [V]','diode current [A]','diode behaviour')
#-------------------------------------------------------------------

#' #Modulator circuits

#' We can now use the elements previously defined:

#' ###Single diode modulator

#' SystemC-AMS source code:

#' Single diode mixer: [H](../circuits/singlediode_mixer.h) - 
#' [CPP](../circuits/singlediode_mixer.cpp) - 
#' [(TB) CPP](../circuits/singlediode_mixer_tb.cpp)

#+ echo=False
#-------------------------------------------------------------------
cct = Circuit("""
   D1 LO_i inner; right=2
   W inner _innerlp; right=1
   TR1 _innerlp _if A; right=2, l=lowpass
   W _if IF_i; right=1
   W inner _innerhp; down=2
   W _innerhp _innerhp2; right=1
   TR2 _innerhp2 _rf A; right=2, l=highpass
   W _rf RF_o; right=1
   """)
cct.draw()
#-------------------------------------------------------------------

#' </BR>
#' Transient plot of the modulated RF signal:

#+ echo=False
#-------------------------------------------------------------------
dat4 = []
dat4 = np.loadtxt('../singlemix_tran.tab', comments='%', unpack=True)
t = dat4[0]
s = dat4[3]
size = len(t)
tran_plot(t,s,6e-3,6.4e-3,'RF_o','RF voltage over time')
#-------------------------------------------------------------------

#' </BR>
#' Some spectrum plots of this circuitry and detailed measurements:

#+ echo=False
#-------------------------------------------------------------------
dat4 = []
dat4 = np.loadtxt('../singlemix_tran.tab', comments='%', unpack=True)

dbm=fft_new(t[int(size/2):size]-t[int(size/2)],s[int(size/2):size],[0,5e6],[-120,00],'RF overall power spectrum')

lores_dbm=fft_new(t[int(size/2):size]-t[int(size/2)],s[int(size/2):size],[0.98e6,1.02e6],[-50,-10],'RF detail power spectrum (residual LO)',1e6)
ifres_dbm=fft_new(t[int(size/2):size]-t[int(size/2)],s[int(size/2):size],[0,10e3],[-150,-50],'RF detail power spectrum (residual IF)',5e3)
usb_dbm=fft_new(t[int(size/2):size]-t[int(size/2)],s[int(size/2):size],[],[],'',1.005e6)
lsb_dbm=fft_new(t[int(size/2):size]-t[int(size/2)],s[int(size/2):size],[],[],'',0.995e6)

s = dat4[5]
lo_dbm=fft_new(t[int(size/2):size]-t[int(size/2)],s[int(size/2):size],[0.99e6,1.01e6],[0,20],'LO input power',1e6)

s = dat4[6]
if_dbm=fft_new(t[int(size/2):size]-t[int(size/2)],s[int(size/2):size],[0,10e3],[-20,0],'IF input power',5e3)
#-------------------------------------------------------------------

#' Extracting some relevant data:

#+ echo=False
#-------------------------------------------------------------------
print('residual LO power: '+str(lores_dbm)+' dBm')
print('residual IF power: '+str(ifres_dbm)+' dBm')
print('USB power: '+str(usb_dbm)+' dBm')
print('LSB power: '+str(lsb_dbm)+' dBm')
print('LO input power '+str(lo_dbm)+' dBm')
print('IF input power '+str(if_dbm)+' dBm')
#-------------------------------------------------------------------


#' More to follow...