From 54f1c18e07800f6a821c8a0e3834a2554f3b1b29 Mon Sep 17 00:00:00 2001 From: Luke Date: Wed, 1 Aug 2018 13:05:02 -0700 Subject: [PATCH] Added RMS error plotting to measured results. Theoretical is still a WiP. --- LPRDefaultPlotting.py | 2 + comparisonPlots.py | 88 +++++++++++++++++++++-- parsePy.py | 161 ++++++++++++++++++++++++++++++++++-------- runParse.sh | 2 +- tankComputers.py | 24 +++++-- 5 files changed, 235 insertions(+), 42 deletions(-) diff --git a/LPRDefaultPlotting.py b/LPRDefaultPlotting.py index 09b3fd2..071b661 100644 --- a/LPRDefaultPlotting.py +++ b/LPRDefaultPlotting.py @@ -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 = [ diff --git a/comparisonPlots.py b/comparisonPlots.py index 416c1fa..307abb9 100755 --- a/comparisonPlots.py +++ b/comparisonPlots.py @@ -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] diff --git a/parsePy.py b/parsePy.py index 3e704a3..24f08e8 100755 --- a/parsePy.py +++ b/parsePy.py @@ -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() diff --git a/runParse.sh b/runParse.sh index 6a34ff1..301c582 100755 --- a/runParse.sh +++ b/runParse.sh @@ -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/ diff --git a/tankComputers.py b/tankComputers.py index 631191b..190a249 100644 --- a/tankComputers.py +++ b/tankComputers.py @@ -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