Added RMS error plotting to measured results. Theoretical is still a WiP.

This commit is contained in:
Luke 2018-08-01 13:05:02 -07:00
parent 852f4cad1d
commit 54f1c18e07
5 changed files with 235 additions and 42 deletions

View file

@ -13,6 +13,7 @@ from matplotlib import rcParams, pyplot as pp
from cycler import cycler
POLAR_YLIM_CONST=(-18,-6)
POLAR_YLIM_CONST_MEAS=(-22,-10)
POLAR_YLIM_CONST_ALT=(-32,-6)
fcFontList = FM.get_fontconfig_fonts()
@ -54,6 +55,7 @@ rcParams['mathtext.bf'] = 'serif:bold'
rcParams['mathtext.sf'] = 'serif'
rcParams['mathtext.tt'] = 'monospace'
rcParams['lines.linewidth'] = 1.0
#rcParams['axes.grid'] = True
# axes.prop_cycle
COLOR_CYCLE_LIST = [

View file

@ -66,12 +66,12 @@ from tankComputers import *
freq_pts = 501
S1=TankGlobals.ampSystem()
B=TankGlobals.bufferSystem()
S1.bw_plt=2
f=FreqClass(freq_pts, S1.f0, S1.bw_plt)
S1.q1_L = 15
S2 = copy.deepcopy(S1)
gain_variation = +4 # dB
gain_variation = +5 # dB
@ -105,8 +105,11 @@ tf_r_ang_ideal1 = wrap_rads(np.concatenate((-S1.phase_swp, -np.pi - S1.phase_swp
tf_r_ang_ideal2 = wrap_rads(np.concatenate((-S2.phase_swp, -np.pi - S2.phase_swp)))
tf_r_ang1 = np.angle(tf_r1)
tf_r_ang2 = np.angle(tf_r2)
tf_r_ang_rms1 = np.sqrt(np.mean(np.power(tf_r_ang1-tf_r_ang_ideal1,2),0))
tf_r_ang_rms2 = np.sqrt(np.mean(np.power(tf_r_ang2-tf_r_ang_ideal2,2),0))
#tf_r_ang_rms1 = np.sqrt(np.mean(np.power(tf_r_ang1-tf_r_ang_ideal1,2),0))
#tf_r_ang_rms2 = np.sqrt(np.mean(np.power(tf_r_ang2-tf_r_ang_ideal2,2),0))
tf_r_ang_rms1_f=delta_rms(tf_r_ang1, 2*np.pi/16)
tf_r_ang_rms2_f=delta_rms(tf_r_ang2, 2*np.pi/16)
################################################################################
# Compute RMS phase error relative to ideal reference across plotting bandwidth
@ -115,15 +118,13 @@ tf_r_ang_rms2 = np.sqrt(np.mean(np.power(tf_r_ang2-tf_r_ang_ideal2,2),0))
(bw_ang2, rms_ang_swp2)=rms_v_bw(tf_r_ang2-tf_r_ang_ideal2, S2.bw_plt)
(bw_mag2, rms_gain_swp2)=rms_v_bw(tf_r2, S2.bw_plt)
(y_buf, tf_buf) = B.compute_ref(f)
################################################################################
################################################################################
################################################################################
#mgr = pp.get_current_fig_manager()
################################################################################
if 3 in plot_list or 13 in plot_list:
if 3 in plot_list:
h3 = [pp.figure() for x in range(2)]
ax3a = h3[0].subplots(1,2)
ax3b = h3[1].subplots(1,2)
@ -195,3 +196,76 @@ if 3 in plot_list or 13 in plot_list:
else:
#mgr.window.geometry(default_window_position[0])
[hT.show() for hT in h3]
if 4 in plot_list:
h4 = [pp.figure() for x in range(2)]
ax4a = h4[0].subplots(1,2)
ax4b = h4[1].subplots(1,2)
ax4 = np.concatenate((ax4a, ax4b))
#ax4[0].plot(bw_mag1,dB20(rms_gain_swp1))
#ax4[1].plot(bw_mag2,dB20(rms_gain_swp2))
ax4[2].plot(f.hz,tf_r_ang_rms1_f*180/np.pi)
ax4[3].plot(f.hz,tf_r_ang_rms2_f*180/np.pi)
h4[0].suptitle('RMS Gain Error')
h4[1].suptitle('RMS Phase Error')
#ax4[0].set_title('RMS Gain Error')
ax4[0].set_ylabel('Gain Error (dB)')
#ax4[2].set_title('RMS Phase Error')
ax4[2].set_ylabel('Phase Error (deg)')
#ax4[1].set_title('RMS Gain Error w/GV')
ax4[1].set_ylabel('Gain Error (dB)')
#ax4[3].set_title('RMS Phase Error w/GV')
ax4[3].set_ylabel('Phase Error (deg)')
# Match Axes
limSetGain = [axT.get_ylim() for axT in ax4[:2]]
limSetPhase = [axT.get_ylim() for axT in ax4[2:]]
limSetGain = (np.min(limSetGain), np.max(limSetGain))
limSetPhase = (np.min(limSetPhase), np.max(limSetPhase))
for axT in ax4[:2]:
axT.set_ylim(limSetGain)
for axT in ax4[2:]:
axT.set_ylim(limSetPhase)
for axT in ax4[[1,3]]:
LPRDefaultPlotting.axAnnotateCorner(axT,
'%g dB gain variation' % (gain_variation), corner=2, ratio=0.04)
axT.yaxis.tick_right()
axT.yaxis.label_position='right'
axT.yaxis.labelpad = axT.yaxis.labelpad + axT.yaxis.label.get_size()
for axT in ax4[[0,2]]:
LPRDefaultPlotting.axAnnotateCorner(axT,
'%g dB gain variation' % (0), corner=2, ratio=0.04)
for axT in ax4:
axT.grid()
axT.set_xlim((np.min(f.hz),np.max(f.hz)))
axT.set_xlabel('Frequency (GHz)')
[hT.tight_layout() for hT in h4]
[hT.tight_layout() for hT in h4]
# Make XY mirror positions
for i in [0,2]:
p0 = ax4[i].get_position()
p1 = ax4[i+1].get_position()
p1.x1 = 1 - p0.x0
p1.x0 = 1 - p0.x1
ax4[i+1].set_position(p1)
for axT in ax4:
p=axT.get_position()
p.y1=0.88
axT.set_position(p)
if args.save:
h4[0].savefig('%s/%s.%s' % (figdir, 'dual_040-RMSGain', fig_ext))
h4[1].savefig('%s/%s.%s' % (figdir, 'dual_041-RMSPhase', fig_ext))
if HEADLESS:
pp.close()
else:
#mgr.window.geometry(default_window_position[0])
[hT.show() for hT in h4]

View file

@ -15,7 +15,7 @@ args_parser.add_argument('--polar','-p', action='store_true',
help='do polar plotting (wide bandwidth)')
args_parser.add_argument('--headless','-q', action='store_true',
help='Remain neadless even if we aren\'t saving files.')
args_parser.add_argument('-n', type=int, default=3,
args_parser.add_argument('-n', type=int, default=4,
help='plot testing number')
args = args_parser.parse_args()
@ -36,6 +36,7 @@ import skrf as rf
from scipy.io import loadmat
from collections import namedtuple
import LPRDefaultPlotting
from tankComputers import *
import re
import json
################################################################################
@ -70,8 +71,10 @@ class MeasurementConfig(namedtuple('config', ['r','c','inv','bias'])):
@property
def fn_str(self):
return "C%02d_R%1d_I%1d_B%0.4f" % (self.c, self.r, self.inv, self.bias)
Measurement = namedtuple('measurement', ['cfg','gain','phase','f','s21', 'slope'])
Measurement = namedtuple('measurement', ['cfg', 'pwr','gain','phase','f','s21', 'slope'])
plottingBandwidthMax = 2.01
plottingBandwidthFreq = 28+np.array([-1,1])*0.5*plottingBandwidthMax
slopeBandwidthMax = 1
slopeBandwidthFreq = 28+np.array([-1,1])*0.5*slopeBandwidthMax
@ -90,7 +93,7 @@ BDE_list.append(BDE(
np.array([ 4.06488853e-03, -5.11527396e-01, 2.53053550e+01]),
np.array([-1.62202706e-03, 6.94343608e-01, -1.80381551e+02]),
-60,
'S02bB_C+02dB_M0'
'S02bB_C+00dB_M0'
))
# 2018-05-16
BDE_list.append(BDE(
@ -98,7 +101,7 @@ BDE_list.append(BDE(
np.array([ 4.08875413e-03, -5.13017311e-01, 2.54047949e+01]),
np.array([-1.29541398e-03, 6.74431785e-01, -1.80127388e+02]),
-60,
'S02bB_C+02dB_M0'
'S02bB_C+00dB_M0'
))
# 2018-05-21
#PolyGain=np.array( [ 4.08875413e-03, -5.13017311e-01, 2.54047949e+01])
@ -108,7 +111,7 @@ BDE_list.append(BDE(
np.array([ 4.08875413e-03, -5.13017311e-01, 2.54047949e+01]),
np.array([-1.29541398e-03, 6.74431785e-01, -1.80127388e+02]),
-60,
'S02bB_C+02dB_M0'
'S02bB_C+00dB_M0'
))
# 2018-05-25
#PolyGain=np.array( [ 4.06488853e-03, -5.11527396e-01, 2.53053550e+01])
@ -118,7 +121,7 @@ BDE_list.append(BDE(
np.array([ 4.06488853e-03, -5.11527396e-01, 2.53053550e+01]),
np.array([-1.62202706e-03, 6.94343608e-01, -1.80381551e+02]),
-70,
'S02bB_C+06dB_M0'
'S02bB_C+00dB_M0'
))
source_directory='fromMat/%s_mat/' % SRC_DATA_NAME
@ -131,6 +134,29 @@ for BDEx in BDE_list:
FamStr=BDEx.mstr
break
with open(SRC_DATA_SUMMARY, 'r') as h_sumDat:
sumDat = json.load(h_sumDat)
def fetchSumDat_pwr(cfg):
global sumDat
mR = np.array(sumDat['r']) == cfg.r
mC = np.array(sumDat['c']) == cfg.c
mI = np.array(sumDat['inv']) == cfg.inv
mB = np.abs(np.array(sumDat['bias_dp_set'])-cfg.bias) < 0.0005
ind = np.squeeze(np.where(np.all((mR,mC,mI,mB),0)))
if ind.size == 0:
print("ERROR EVERYTHING IS BROKEN! AND i'M TIRED")
return -1
else:
return sumDat['ivdd'][ind]*sumDat['vdd'][ind]
def sumTuple_avgMinMax(data_list):
existing_data = []
for datum in data_list:
existing_data.extend([np.mean(datum), np.min(datum), np.max(datum)])
return tuple(existing_data)
combined_rms=np.array([])
for filename in os.listdir(source_directory):
filename=source_directory+filename
group_filename_string = filename.split('/')[-1][:-4]
@ -147,24 +173,31 @@ for filename in os.listdir(source_directory):
s2p_file = rf.Network(SRC_DATA_LOC + (FILE_PAT % pt.fn_str) )
freq = np.squeeze(s2p_file.f*1e-9)
inds_keep = np.where(np.all((freq >= plottingBandwidthFreq[0],
freq <= plottingBandwidthFreq[1]),0))
sdat_raw = np.squeeze(s2p_file.s21.s)
freq = freq[inds_keep]
sdat_raw = sdat_raw[inds_keep]
buffer_gain = np.polyval(PolyGain,freq)
buffer_phase = np.polyval(PolyPhase,freq)
buffer_phase = buffer_phase - np.mean(buffer_phase) + \
PhaseFixedRotationFactor*np.pi/180
buffer_sdat = np.power(10,buffer_gain/20)*np.exp(1j*buffer_phase)
sdat = np.squeeze(s2p_file.s21.s)/buffer_sdat
sdat = sdat_raw/buffer_sdat
slope_valid_inds = np.where(np.all((freq >= slopeBandwidthFreq[0],
freq <= slopeBandwidthFreq[1]),0))
sub_angles = np.unwrap(np.angle(sdat[slope_valid_inds]))*180/np.pi
sub_freq = freq[slope_valid_inds]-np.mean(freq[slope_valid_inds])
slope = np.polyfit(sub_freq,sub_angles-np.mean(sub_angles),1)[0]
index = np.squeeze(np.argwhere(freq==28))
collectedData.append(Measurement(pt,
dB20(sdat[index]),
ang_deg(sdat[index]),
freq, sdat, slope))
index_f0 = np.squeeze(np.argwhere(freq==28))
collectedData.append(Measurement(cfg=pt, pwr=fetchSumDat_pwr(pt),
gain=dB20(sdat[index_f0]),
phase=ang_deg(sdat[index_f0]),
f=freq, s21=sdat, slope=slope))
# Find the indicies close to 0 and 180 as my reference curves
phis = np.array([s.phase for s in collectedData])
@ -172,44 +205,105 @@ for filename in os.listdir(source_directory):
slope_list = np.array([s.slope for s in collectedData])
slope_avg = np.mean(slope_list[best_slopes])
h=pp.figure()
ref_index = np.argmin(np.abs(phis))
unwrapped_ref_phase = 180/np.pi*np.unwrap(ang(collectedData[ref_index].s21))
if args.polar:
h=pp.figure()
ax=h.add_subplot(1,1,1, projection='polar')
else:
h2=pp.figure()
h=pp.figure(figsize=(3.4*figScaleSize, 4.5*figScaleSize))
h2=pp.figure(figsize=(3.4*figScaleSize, 2.8*figScaleSize))
ax=h.subplots(2,1)
ax = np.append(ax, h2.subplots(1,1))
print("---------------------||------------------------------")
print(" _C R I _Bias_ || Gain Phase ")
print("---------------------||------------------------------")
ax = np.append(ax, ax[1].twinx())
ax = np.append(ax, ax[2].twinx())
summary_msg = \
"/---------------------\/----------------------------------------\\\n"\
"| _C R I _Bias_ || Gain Phase Power |\n"\
"|---------------------||----------------------------------------|\n"
all_sdat = np.column_stack([imeas.s21 for imeas in collectedData])
ang_rms = delta_rms(np.angle(all_sdat), 2*np.pi/16)*180/np.pi
for imeas in collectedData:
if args.polar:
#ax.plot(ang(imeas.s21)-buffer_phase, dB20(imeas.s21)-buffer_gain)
ax.plot(ang(imeas.s21), dB20(imeas.s21))
else:
#ax[0].plot(imeas.f, dB20(imeas.s21)-buffer_gain)
ax[0].plot(imeas.f, dB20(imeas.s21))
#unwrapped_phase = 180/np.pi*np.unwrap(ang(imeas.s21)-buffer_phase)
#ax[1].plot(imeas.f, unwrapped_phase)
unwrapped_phase = 180/np.pi*np.unwrap(ang(imeas.s21))
ax[1].plot(imeas.f, unwrapped_phase)
slope_relative = (imeas.f-28)*slope_avg
ax[2].plot(imeas.f, unwrapped_phase-slope_relative)
print(" %2d %d %d %.4f || %+7.1f dB %+9.2f deg" % \
#ax[1].plot(imeas.f, unwrapped_phase)
relative_phase_curve = unwrapped_phase-unwrapped_ref_phase
if np.any(relative_phase_curve < 0):
relative_phase_curve += 360
#relative_phase_curve -= 180-22.5/2
ax[1].plot(imeas.f, relative_phase_curve)
#slope_relative = (imeas.f-28)*slope_avg
#ax[2].plot(imeas.f, unwrapped_phase-slope_relative)
ax[2].plot(imeas.f, relative_phase_curve)
pwr_overage = int(2*(imeas.pwr*1e3 - 10))
pwr_string = (int(pwr_overage/2)*"=") + (np.mod(pwr_overage,2)*">")
summary_msg += "| %2d %d %d %.4f || "\
" %+7.1f dB %+9.2f deg %4.1f mW |%s\n" % \
(imeas.cfg.c, imeas.cfg.r, imeas.cfg.inv, imeas.cfg.bias, \
imeas.gain, imeas.phase))
print("---------------------||------------------------------")
imeas.gain, imeas.phase, imeas.pwr*1e3, pwr_string)
summary_msg += \
"\_____________________/\________________________________________/\n"
pwr_list=np.array([imeas.pwr*1e3 for imeas in collectedData])
gain_list=np.array([imeas.gain for imeas in collectedData])
summary_msg += \
"/ \\\n" \
"|===> Power: % 7.1f mW (% 7.1f mW - % 7.1f mW) |\n" \
"|===> Gain: %+7.1f dB (%+7.1f dB - %+7.1f dB) | \n" \
"|===> RMS: %6.1f deg (%6.1f deg - %6.1f deg) | \n" \
"\_______________________________________________________________/" % \
(sumTuple_avgMinMax([pwr_list, gain_list, ang_rms]))
if args.polar:
ax.set_ylim(LPRDefaultPlotting.POLAR_YLIM_CONST)
ax.set_ylim(LPRDefaultPlotting.POLAR_YLIM_CONST_MEAS)
if args.polar:
ax.set_title('Measured Performance')
else:
# Usually this also has crappy lower ylimits, so we fix that here.
# get ALL THE LOWER bounds
np.min([np.min(line.get_ydata()) for line in ax[2].get_lines()])
ax[0].set_title('Measured Performance')
ax[0].set_ylabel('Gain (dB)');
ax[1].set_ylabel('Phase (deg)');
ax[2].set_ylabel('Phase (deg)');
ax[0].set_ylabel('Gain (dB)')
ax[1].set_ylabel('Relative Phase (deg)')
ax[2].set_ylabel('Relative Phase (deg)')
ax[2].set_title('Relative Phase')
for i in range(3,5):
aT=ax[i]
aR=ax[i-2]
# make the ticks, and the y-axis line up in a tidy manner
# Recall that the ylimits should be 0-360 basically.
aT.set_ylabel('RMS Error (deg)')
aT.plot(imeas.f, ang_rms)
# The goal is to take the usual step size of 50,
# and then equate that with a 1-degree step in RMS Error
# and to then adjust the y-limit of the twin-axis to align
# the grid markers
if False:
yRscl=np.diff(aR.get_yticks()[-2:])
yTscl=np.diff(aT.get_yticks()[-2:])
# Now find the ratio of the ylimits margin verses their
# extreme tick marks.
yRmrks = aR.get_yticks()[[0,-1]]
yTmrks = aT.get_yticks()[[0,-1]]
tickTotal = max(len(aT.get_yticks()), len(aR.get_yticks()))
yRover = (aR.get_ylim()-yRmrks)/yRscl
yTover = (aT.get_ylim()-yTmrks)/yTscl
yRTover = np.stack((yRover,yTover))
yXover = np.array([np.min(yRTover[:,0]), np.max(yRTover[:,1])])
aR.set_ylim(yRscl*yXover + yRmrks)
aT.set_ylim(yTscl*yXover + yTmrks)
aT.set_ylim(aR.get_ylim()/np.array(50)+3)
aT.grid()
aT.get_lines()[0].set_linewidth(2.0)
aT.get_lines()[0].set_linestyle('-.')
aT.get_lines()[0].set_color('black')
for aT in ax:
aT.set_xlabel('Frequency (GHz)')
aT.grid()
@ -225,6 +319,11 @@ for filename in os.listdir(source_directory):
if not args.polar:
h2.tight_layout()
if args.save:
with open('%s/Summary-%s-%s.txt' % (figdir, FamStr,
group_filename_string), 'w') as summary_file:
summary_file.write(summary_msg)
summary_file.write("\n")
summary_file.close()
if args.polar:
h.savefig('%s/PolarGain-%s-%s.%s' % (figdir, FamStr,
group_filename_string, fig_ext))
@ -233,6 +332,8 @@ for filename in os.listdir(source_directory):
group_filename_string, fig_ext))
h2.savefig('%s/RelStdPlots-%s-%s.%s' % (figdir, FamStr,
group_filename_string, fig_ext))
else:
print(summary_msg)
if HEADLESS:
if not args.polar:
pp.close()

View file

@ -15,6 +15,6 @@ for n in $(seq 1 4); do
done
while [[ $(jobs -lr | wc -l) -gt 0 ]]; do sleep 0.1; done
SELECT_STRING="S02bB_C+03dB"
SELECT_STRING="S02bB_C+00dB"
rsync -aPv "figures-measured/"*"${SELECT_STRING}"* ../tex/figures-measured/

View file

@ -15,15 +15,15 @@ def ang_unwrap(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 setLimitsTicks(ax, data, steps):
targs = np.array([1, 2, 4, 5, 10, 20, 30, 50, 60, 100, 250, 1000])
lo = np.min(data)
@ -36,7 +36,7 @@ def setLimitsTicks(ax, data, steps):
marks = np.arange(0,steps+1)*step_size + lo
ax.set_ylim((lo,hi))
ax.set_yticks(marks)
def rms_v_bw(err_sig, bandwidth_scale=1):
"""compute the rms vs bandwidth assuming a fixed center frequency"""
# First compute the error power
@ -63,3 +63,19 @@ def rms_v_bw(err_sig, bandwidth_scale=1):
rms = np.sqrt(np.cumsum(folded,0) / (ind*np.ones((folded.shape[1],1))).T )
return (frac_step*bandwidth_scale, rms)
def delta_rms(signal, reference_delta, wrap_point=2*np.pi):
"""compute the rms difference between various states and a reference"""
# First compute the matrix difference including folding
signal_delta = np.column_stack((
signal[:,1:]-signal[:,:-1],
signal[:,0]-signal[:,-1]
))
signal_delta = np.where(signal_delta>wrap_point/2, \
signal_delta-wrap_point, signal_delta)
signal_delta = np.where(signal_delta<-wrap_point/2, \
signal_delta+wrap_point, signal_delta)
signal_error = np.abs(signal_delta)-reference_delta
signal_rms = np.sqrt(np.mean(np.power(signal_error,2),1))
return signal_rms