Major refactor to ease duplicate computations and plotting

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

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
+class ampSystem:
-bw0		= 8 # assumed tuning range (GHz)
+	"""define global (hardware descriptive) variables for use in our system."""
-bw_plt	= 4 # Plotting range (GHz)
+	def __init__(self, quiet=False):
-fbw		= bw0/f0 # fractional bandwidth
+		self.f0		= 28 # design frequency (GHz)
-
+		self.bw0	= 8 # assumed extreme tuning range (GHz)
-frequency_sweep_steps = 101
+		self.bw_plt	= 4 # Plotting range (GHz)
-gamma_sweep_steps = 8
-
-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
 
 		# Configuration Of Hardware
 		#####
-q1_L	= 20
+		self.q1_L	= 25
-q1_C	= 7
+		self.q1_C	= 8
-l1		= 180e-3 # nH
+		self.l1		= 140e-3 # nH
-gm1		= 25e-3 # S
+		self.gm1	= 25e-3 # S
 
-# Compute frequency sweep
+		self._gamma_steps=8
-#####
+		self._gamma_cap_ratio = 0.997
-w0		= f0*2*np.pi
+		self.alpha_min=1
-fbw_plt	= bw_plt/f0
+		if not quiet:
-delta	= np.linspace(-fbw_plt/2,fbw_plt/2,frequency_sweep_steps)
+			## Report System Descrption
-w		= w0*(1+delta)
+			print('  L1 = %.3fpH, C1 = %.3ffF' % (1e3*self.l1, 1e6*self.c1))
-f		= f0*(1+delta)
+			print('    Rp = %.3f Ohm' % (1/self.g1))
-jw		= 1j*w
+			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)
 
-##################
 	# Compute system 
 	#####
-c1		= 1/(w0*w0*l1)
+	@property
-g1_L	= 1 / (q1_L*w0*l1)
+	def c1(self):
-g1_C	= w0 * c1 / q1_C
+		return 1/(self.w0*self.w0*self.l1)
-g1		= g1_L + g1_C
+	@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
+
+	@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 = g1 * np.sqrt(l1/c1)
+		gamma_max = 1/(self.alpha_min*self.Q1)
-if gamma > (gamma_limit_ratio * gamma_max):
+		if gamma > (self._gamma_cap_ratio * gamma_max):
-	print("==> WARN: Gamma to large, reset to %.3f (was %.3f) <==" % \
+			if not self._gamma_warn:
-		(gamma_limit_ratio * gamma_max, gamma))
+				self._gamma_warn = True
-	gamma = gamma_limit_ratio * gamma_max
+				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)))
+
+
 

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

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)
+# We want a smooth transition out to alpha. So For now assume a squares
-
+# weighting out to the maximum alpha at the edges.
-# Linear based gamma spacing
+gain_variation = -8*0	# dB
-#gamma_swp = np.linspace(-gamma,gamma,gamma_sweep_steps)
+S.alpha_min = dB2Vlt(gain_variation)
-
-# 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)
 
 # 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_boost = (g1_swp - S.g1)
-g1_ratio = -g1_boost / gm1
+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)
+# Generate a reference implementation
-for itune,gamma_tune in enumerate(gamma_swp):
+(y_tank, tf) = S.compute_block(f)
-	c1_tune = c1_swp[itune]
+(_, tf_ref) = S.compute_ref(f)
-	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
 # double to describe with perfect inversion stage
 tf = np.column_stack((tf,-tf))
 
-ref_index = int(gamma_swp.shape[0]/2)
+# compute the relative transfer function thus giving us flat phase, and
-tf_r = tf / (tf[:,ref_index]*np.ones((tf.shape[1],1))).T
+# flat (ideally) gain response if our system perfectly matches the reference
-y_tank = y_tank.T
+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))
+ax3.plot(f.hz,dB20(tf))
-ax4.plot(f,ang_unwrap(tf))
+ax4.plot(f.hz,ang_unwrap(tf))
 
 ################################################################################
-ax8 = h2.add_subplot(2,1,1)
+ax6 = h2.add_subplot(2,1,1)
-ax9 = h2.add_subplot(2,1,2)
+ax7 = h2.add_subplot(2,1,2)
-ax8.plot(f,dB20(tf_r))
+ax6.plot(f.hz,dB20(tf_r))
-ax9.plot(f,ang_unwrap(tf_r.T).T)
+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')
+ax6.set_title('TF Relative Gain')
-ax8.set_ylabel('Relative Gain (dB)')
+ax6.set_ylabel('Relative Gain (dB)')
-ax9.set_title('TF Relative Phase')
+ax7.set_title('TF Relative Phase')
-ax9.set_ylabel('Relative Phase (deg)')
+ax7.set_ylabel('Relative Phase (deg)')
-for ax_T in [ax3, ax4, ax8, ax9]:
+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()