From 539c2f74818087f12bce85532bde0396a1bd4cab Mon Sep 17 00:00:00 2001 From: Luke Date: Fri, 20 Jul 2018 19:07:01 -0700 Subject: [PATCH] Major refactor to ease duplicate computations and plotting --- FreqClass.py | 62 ++++++++++++++++ TankGlobals.py | 188 +++++++++++++++++++++++++++++++++++++---------- tankComputers.py | 52 +++++++++++++ tankPlot.py | 144 +++++++++++++++--------------------- 4 files changed, 326 insertions(+), 120 deletions(-) create mode 100644 FreqClass.py create mode 100644 tankComputers.py diff --git a/FreqClass.py b/FreqClass.py new file mode 100644 index 0000000..6a78c00 --- /dev/null +++ b/FreqClass.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python3 +import numpy as np + +class FreqClass: + def __init__(self, steps, f0, bw): + self.f0 = f0 + self._bw = bw + self._steps = steps; + self._update_delta() + + def _update_delta(self): + self._delta = self._bw/self.f0*np.linspace(-1/2,1/2,self._steps) + + def __repr__(self): + return self.__str__() + def __str__(self): + return "%gGHz, %gGHz BW sweep [%d points]" % \ + (self.f0, self._bw, self._steps) + + @property + def hz_range(self): + return (np.min(self.hz), np.max(self.hz)) + + @property + def delta(self): + return self._delta + @property + def bw(self): + return self._bw + @bw.setter + def bw(self, bw): + self._bw = bw + self._update_delta() + + @property + def steps(self): + return self._steps + @steps.setter + def steps(self, steps): + self._steps = steps + self._update_delta() + + @property + def hz(self): + return self.f0*(1+self._delta) + @property + def f(self): + return self.hz + @property + def rad(self): + return 2*np.pi*self.f0*(1+self._delta) + @property + def w(self): + return self.rad + @property + def jw(self): + return 2j*np.pi*self.f0*(1+self._delta) + @property + def delta(self): + return self._delta + + diff --git a/TankGlobals.py b/TankGlobals.py index 6129836..a73a4d2 100644 --- a/TankGlobals.py +++ b/TankGlobals.py @@ -1,50 +1,164 @@ #!/usr/bin/env python3 import numpy as np +import sys ################################################################################ +# BEWARE, FOR BEYOND THIS POINT THERE BE DRAGONS! THIS IS ONLY FOR EASE OF +# GENERATING ACADEMIC PUBLICATIONS AND FIGURES, NEVER DO THIS SHIT! +################################################################################ + +def g1_map_default(system): + # compute correction factor for g1 that will produce common gain at f0 + g1_swp = system.g1 * np.sin(np.pi/2-system.phase_swp) / system.alpha_swp + return g1_swp + # Operating Enviornment ##### -f0 = 28 -bw0 = 8 # assumed tuning range (GHz) -bw_plt = 4 # Plotting range (GHz) -fbw = bw0/f0 # fractional bandwidth +class ampSystem: + """define global (hardware descriptive) variables for use in our system.""" + def __init__(self, quiet=False): + self.f0 = 28 # design frequency (GHz) + self.bw0 = 8 # assumed extreme tuning range (GHz) + self.bw_plt = 4 # Plotting range (GHz) -frequency_sweep_steps = 101 -gamma_sweep_steps = 8 + # Configuration Of Hardware + ##### + self.q1_L = 25 + self.q1_C = 8 + self.l1 = 140e-3 # nH + self.gm1 = 25e-3 # S -gamma = 1 - np.power(f0 / (f0 + bw0/2),2) -gamma_limit_ratio = 0.99 # how close gamma can get to theoretical extreme -phase_limit_requested = (1-1/gamma_sweep_steps)*np.pi/2 + self._gamma_steps=8 + self._gamma_cap_ratio = 0.997 + self.alpha_min=1 + if not quiet: + ## Report System Descrption + print(' L1 = %.3fpH, C1 = %.3ffF' % (1e3*self.l1, 1e6*self.c1)) + print(' Rp = %.3f Ohm' % (1/self.g1)) + print(' Q = %.1f' % (self.Q1)) + self._gamma_warn = False + + self._g1_map_function = g1_map_default + + @property + def w0(self): + return self.f0*2*np.pi + @property + def fbw(self): # fractional bandwidth + return self.bw0/self.f0 + @property + def phase_max(self): + return np.pi/2 * (1 - 1/self.gamma_len) -# Configuration Of Hardware -##### -q1_L = 20 -q1_C = 7 -l1 = 180e-3 # nH -gm1 = 25e-3 # S + # Compute system + ##### + @property + def c1(self): + return 1/(self.w0*self.w0*self.l1) + @property + def g1(self): + g1_L = 1 / (self.q1_L*self.w0*self.l1) + g1_C = self.w0 * self.c1 / self.q1_C + return g1_L + g1_C + @property + def Q1(self): + return np.sqrt(self.c1/self.l1)/self.g1 -# Compute frequency sweep -##### -w0 = f0*2*np.pi -fbw_plt = bw_plt/f0 -delta = np.linspace(-fbw_plt/2,fbw_plt/2,frequency_sweep_steps) -w = w0*(1+delta) -f = f0*(1+delta) -jw = 1j*w + @property + def gamma_len(self): + return self._gamma_steps + + @property + def gamma(self): + gamma = 1 - np.power(self.f0 / (self.f0 + self.bw0/2),2) + phase_limit_requested = (1-1/self.gamma_len)*np.pi/2 + + # Verify gamma is valid + ##### + gamma_max = 1/(self.alpha_min*self.Q1) + if gamma > (self._gamma_cap_ratio * gamma_max): + if not self._gamma_warn: + self._gamma_warn = True + print("==> WARN: Gamma to large, reset to %.1f%% (was %.1f%%) <==" % \ + (100*self._gamma_cap_ratio * gamma_max, 100*gamma)) + gamma = self._gamma_cap_ratio * gamma_max + return gamma + + @property + def alpha_swp(self): + range_partial = np.ceil(self.gamma_len/2) + lhs = np.linspace(np.sqrt(self.alpha_min),1, range_partial) + rhs = np.flip(lhs,0) + swp = np.concatenate((lhs,rhs[1:])) if np.mod(self.gamma_len,2) == 1 \ + else np.concatenate((lhs,rhs)) + return np.power(swp,2) + + @property + def gamma_swp(self): + return np.cos(np.pi/2-self.phase_swp) / self.Q1 / self.alpha_swp + @property + def phase_swp(self): + #def phaseSweepGenerate(g1, gamma, c, l, phase_extreme, phase_steps): + # Linear PHASE gamma spacing + # First compute the most extreme phase given the extreme gamma + # if gamma is tuned to the limit, and we want to match the gain performance, + # then this is the required tuned g1 value. + gamma = self.gamma + g1_limit = np.sqrt(np.power(self.g1,2) - np.power(gamma,2)*self.c1/self.l1) + # This implies a Q in that particular setting + Q_limit = self.Q1*self.g1/g1_limit + # given this !, I compute the delta phase at that point. + phase_limit = np.pi/2 - np.arctan(1/(Q_limit*gamma)) + + phase_swp = np.linspace(-1,1,self.gamma_len) * self.phase_max + + if phase_limit < self.phase_max: + print( "==> ERROR: Phase Beyond bounds. Some states will be ignored") + print( " %.3f requested\n" + " %.3f hardware limit" % \ + (180/np.pi*self.phase_max, 180/np.pi*abs(phase_limit))) + print( " To increase tuning range, gamma must rise or native Q must rise") + phase_swp = np.where(phase_swp > phase_limit, phase_swp, np.NaN) + + # This gives us our equal phase spacing points + return phase_swp + + @property + def c1_swp(self): + return self.c1 * (1 + self.gamma_swp) + + def set_g1_swp(self, g1_swp_function): + self._g1_map_function = g1_swp_function + + @property + def g1_swp(self): + return self._g1_map_function(self) + + def compute_block(self, f_dat): + g1_swp = self.g1_swp + c1_swp = self.c1_swp + y_tank = np.zeros((self.gamma_len,f_dat.steps), dtype=complex) + tf = np.zeros((self.gamma_len,f_dat.steps), dtype=complex) + for itune,gamma_tune in enumerate(self.gamma_swp): + c1_tune = c1_swp[itune] + g1_tune = g1_swp[itune] + y_tank[itune,:] = g1_tune + f_dat.jw*c1_tune + 1/(f_dat.jw * self.l1) + tf[itune,:] = self.__class__.tf_compute(f_dat.delta, gamma_tune, g1_tune, self.gm1, self.l1, self.c1) + + tf = tf.T + return (y_tank, tf) + + def compute_ref(self, f_dat): + y_tank = self.g1 + f_dat.jw*self.c1 + 1/(f_dat.jw * self.l1) + tf = self.__class__.tf_compute(f_dat.delta, 0, self.g1, self.gm1, self.l1, self.c1) + return (y_tank, tf) + + @classmethod + def tf_compute(cls, delta, gamma, gx, gm, l, c): + Q = np.sqrt(c/l)/gx + return gm / gx \ + * 1j*(1+delta) \ + / (1j*(1+delta) + Q*(1-np.power(1+delta,2)*(1+gamma))) -################## -# Compute system -##### -c1 = 1/(w0*w0*l1) -g1_L = 1 / (q1_L*w0*l1) -g1_C = w0 * c1 / q1_C -g1 = g1_L + g1_C -# Verify gamma is valid -##### -gamma_max = g1 * np.sqrt(l1/c1) -if gamma > (gamma_limit_ratio * gamma_max): - print("==> WARN: Gamma to large, reset to %.3f (was %.3f) <==" % \ - (gamma_limit_ratio * gamma_max, gamma)) - gamma = gamma_limit_ratio * gamma_max diff --git a/tankComputers.py b/tankComputers.py new file mode 100644 index 0000000..d178362 --- /dev/null +++ b/tankComputers.py @@ -0,0 +1,52 @@ +#!/usr/bin/env python3 +import numpy as np + +################################################################################ +# Define my helper functions. +def dB20(volt_tf): + """Describe signal gain of a transfer function in dB (i.e. 20log(x))""" + return 20*np.log10(np.abs(volt_tf)) +def ang(volt_tf): + """Describe phase of a transfer function in degrees. Not unwrapped.""" + return 180/np.pi*np.angle(volt_tf) +def ang_unwrap(volt_tf): + """Describe phase of a transfer function in degrees. With unwrapping.""" + return 180/np.pi*np.unwrap(np.angle(volt_tf)) +def dB10(pwr_tf): + """Describe power gain of a transfer function in dB (i.e. 10log(x))""" + return 10*np.log10(np.abs(pwr_tf)) + +def dB2Vlt(dB20_value): + return np.power(10,dB20_value/20) + +def wrap_rads(angles): + return np.mod(angles+np.pi,2*np.pi)-np.pi +def atand(x): + return 180/np.pi*np.arctan(x) + +def rms_v_bw(err_sig, bandwidth_scale=1): + """compute the rms vs bandwidth assuming a fixed center frequency""" + # First compute the error power + err_pwr = np.power(np.abs(err_sig),2) + steps = len(err_pwr) + isodd = True if steps%2 != 0 else False + + # We want to generate the midpoint to the left, and midpoint to the right + # as two distinct sets. + pt_rhs_start = int(np.floor(steps/2)) + pt_lhs_stop = int(np.ceil(steps/2)) + + folded = err_pwr[pt_rhs_start:] + np.flip(err_pwr[:pt_lhs_stop],0) + # Now, we MIGHT have double counted the mid point + # if the length is odd, correct for that + if isodd: folded[0]*=0.5 + + # Now we need an array that describes the number of points used to get here. + # this one turns out to be pretty easy. + frac_step = np.arange(int(not isodd),steps,2)/(steps-1) + ind = 2*np.arange(0,frac_step.shape[0])+1+int(not isodd) + + # Now actually compute the RMS values. First do the running sum + rms = np.sqrt(np.cumsum(folded,0) / (ind*np.ones((folded.shape[1],1))).T ) + return (frac_step*bandwidth_scale, rms) + diff --git a/tankPlot.py b/tankPlot.py index 2d6ce45..964936f 100644 --- a/tankPlot.py +++ b/tankPlot.py @@ -10,25 +10,6 @@ sys.path.append("./pySmithPlot") import smithplot from smithplot import SmithAxes -################################################################################ -# Define my helper functions. -def dB20(volt_tf): - """Describe signal gain of a transfer function in dB (i.e. 20log(x))""" - return 20*np.log10(np.abs(volt_tf)) -def ang(volt_tf): - """Describe phase of a transfer function in degrees. Not unwrapped.""" - return 180/np.pi*np.angle(volt_tf) -def ang_unwrap(volt_tf): - """Describe phase of a transfer function in degrees. With unwrapping.""" - return 180/np.pi*np.unwrap(np.angle(volt_tf)) -def dB10(pwr_tf): - """Describe power gain of a transfer function in dB (i.e. 10log(x))""" - return 10*np.log10(np.abs(pwr_tf)) - -def atan(x): - return 180/np.pi*np.arctan(x) - - ################################################################################ # Override the defaults for this script rcParams['figure.figsize'] = [10,7] @@ -36,76 +17,57 @@ default_window_position=['+20+80', '+120+80'] ################################################################################ # Operating Enviornment (i.e. circuit parameters) -from TankGlobals import * +import TankGlobals +from FreqClass import FreqClass +from tankComputers import * + +S=TankGlobals.ampSystem() +f=FreqClass(501, S.f0, S.bw_plt) ################################################################################ -# Now generate the sweep of resonance tuning (gamma, and capacitance) - -# Linear based gamma spacing -#gamma_swp = np.linspace(-gamma,gamma,gamma_sweep_steps) - -# Linear PHASE gamma spacing -# First compute the most extreme phase given the extreme gamma -g1_limit = np.sqrt( g1*g1 - (gamma*gamma) * c1/l1 ) -K_limit = np.sqrt(c1/l1)*1/g1_limit -phase_limit = np.mod(np.pi/2 - np.arctan( -1/K_limit * 1/gamma ),np.pi) - np.pi - -if abs(phase_limit) < phase_limit_requested: - print("==> WARN: Phase Beyond bounds, leaving at limits. <==") - print("==> %.3f requested, but hardware limit is %.3f <==" % \ - (180/np.pi*phase_limit_requested, 180/np.pi*abs(phase_limit))) - sys.exit(-1) -else: - phase_limit = phase_limit_requested - - -# This gives us our equal phase spacing points -phase_swp = np.linspace(-1,1,gamma_sweep_steps) * phase_limit -# Then use this to compute the gamma steps to produce arbitrary phase given -# our perfect gain constraint. -gamma_swp = np.sign(phase_swp)/np.sqrt(np.power(np.tan(np.pi/2 - phase_swp),2)+1) * g1 / np.sqrt(c1/l1) +# We want a smooth transition out to alpha. So For now assume a squares +# weighting out to the maximum alpha at the edges. +gain_variation = -8*0 # dB +S.alpha_min = dB2Vlt(gain_variation) # compute correction factor for g1 that will produce common gain at f0 -g1_swp = np.sqrt( g1*g1 - (gamma_swp*gamma_swp) * c1/l1 ) +# this is defined as the class default +g1_swp = S.g1_swp # and compute how much of a negative gm this requres, and it's relative # proportion to the gm of the assumed main amplifier gm. -g1_boost = (g1_swp - g1) -g1_ratio = -g1_boost / gm1 +g1_boost = (g1_swp - S.g1) +g1_ratio = -g1_boost / S.gm1 -c1_swp = c1 * (1 + gamma_swp) - -## Report System Descrption -print(' L1 = %.3fpH, C1 = %.3ffF' % (1e3*l1, 1e6*c1)) -print(' Rp = %.3f Ohm' % (1/g1)) print(' Max G1 boost %.2fmS (%.1f%% of gm1)' % \ (1e3*np.max(np.abs(g1_boost)), 100*np.max(g1_ratio))) -y_tank = np.zeros((len(gamma_swp),len(f)), dtype=complex) -tf = np.zeros((len(gamma_swp),len(f)), dtype=complex) -for itune,gamma_tune in enumerate(gamma_swp): - c1_tune = c1_swp[itune] - g1_tune = g1_swp[itune] - K = np.sqrt(c1/l1)/g1_tune - y_tank_tmp = g1_tune + jw*c1_tune + 1/(jw * l1) - y_tank[itune,:] = y_tank_tmp - tf_tmp = gm1 / g1_tune * \ - 1j*(1+delta) / \ - ( 1j*(1+delta) + K*(1 - (1+gamma_tune)*np.power(1+delta,2)) ) - tf[itune,:] = tf_tmp - -tf = tf.T +################################################################################ +# Generate a reference implementation +(y_tank, tf) = S.compute_block(f) +(_, tf_ref) = S.compute_ref(f) # double to describe with perfect inversion stage tf = np.column_stack((tf,-tf)) -ref_index = int(gamma_swp.shape[0]/2) -tf_r = tf / (tf[:,ref_index]*np.ones((tf.shape[1],1))).T -y_tank = y_tank.T +# compute the relative transfer function thus giving us flat phase, and +# flat (ideally) gain response if our system perfectly matches the reference +tf_r = tf / (tf_ref*np.ones((tf.shape[1],1))).T + +# We will also do a direct angle comparison +tf_r_ang_ideal = wrap_rads(np.concatenate((-S.phase_swp, -np.pi - S.phase_swp))) +tf_r_ang = np.angle(tf_r) +tf_r_ang_rms = np.sqrt(np.mean(np.power(tf_r_ang-tf_r_ang_ideal,2),0)) + +y_tank = y_tank.T +################################################################################ +# Compute RMS phase error relative to ideal reference across plotting bandwidth +(bw_ang, rms_ang_swp)=rms_v_bw(tf_r_ang-tf_r_ang_ideal, S.bw_plt) +(bw_mag, rms_gain_swp)=rms_v_bw(tf_r, S.bw_plt) -print(ang(tf[f==28,:])) ################################################################################ h1 = pp.figure() h2 = pp.figure(figsize=(5,7)) +h3 = pp.figure(figsize=(5,7)) mgr = pp.get_current_fig_manager() ################################################################################ ax1 = h1.add_subplot(2,2,1, projection='smith') @@ -115,14 +77,19 @@ ax4 = h1.add_subplot(2,2,4) ax1.plot(y_tank, datatype=SmithAxes.Y_PARAMETER, marker="None") ax2.plot(np.angle(tf), dB20(tf)) -ax3.plot(f,dB20(tf)) -ax4.plot(f,ang_unwrap(tf)) +ax3.plot(f.hz,dB20(tf)) +ax4.plot(f.hz,ang_unwrap(tf)) ################################################################################ -ax8 = h2.add_subplot(2,1,1) -ax9 = h2.add_subplot(2,1,2) -ax8.plot(f,dB20(tf_r)) -ax9.plot(f,ang_unwrap(tf_r.T).T) +ax6 = h2.add_subplot(2,1,1) +ax7 = h2.add_subplot(2,1,2) +ax6.plot(f.hz,dB20(tf_r)) +ax7.plot(f.hz,ang_unwrap(tf_r.T).T) + +ax8 = h3.add_subplot(2,1,1) +ax9 = h3.add_subplot(2,1,2) +ax8.plot(bw_mag,dB20(rms_gain_swp)) +ax9.plot(bw_ang,rms_ang_swp*180/np.pi) ax1.set_title('Tank Impedance') ax2.set_title('Transfer Function') @@ -131,20 +98,31 @@ ax3.set_title('TF Gain') ax3.set_ylabel('Gain (dB)') ax4.set_title('TF Phase') ax4.set_ylabel('Phase (deg)') -ax8.set_title('TF Relative Gain') -ax8.set_ylabel('Relative Gain (dB)') -ax9.set_title('TF Relative Phase') -ax9.set_ylabel('Relative Phase (deg)') -for ax_T in [ax3, ax4, ax8, ax9]: +ax6.set_title('TF Relative Gain') +ax6.set_ylabel('Relative Gain (dB)') +ax7.set_title('TF Relative Phase') +ax7.set_ylabel('Relative Phase (deg)') +for ax_T in [ax3, ax4, ax6, ax7]: ax_T.grid() ax_T.set_xlabel('Freq (GHz)') - ax_T.set_xlim(( np.min(f), np.max(f) )) + ax_T.set_xlim(f.hz_range) + +ax8.set_title('RMS Gain Error') +ax8.set_ylabel('RMS Gain Error (dB)') +ax9.set_title('RMS Phase Error') +ax9.set_ylabel('RMS Phase Error (deg)') +for ax_T in [ax8, ax9]: + ax_T.grid() + ax_T.set_xlim((0,S.bw_plt)) + ax_T.set_xlabel('Bandwidth (GHz)') ################################################################################ h1.tight_layout() h2.tight_layout() +h3.tight_layout() mgr.window.geometry(default_window_position[0]) h1.show() mgr.window.geometry(default_window_position[1]) h2.show() +h3.show()