Linear least squares phase estimation of a sine wave

This snippet extracts the phase out of a noisy sinusoidal signal. The algorithm is a clever linear least squares (LLS) scheme, courtesy of the IEEE waveform measurement standard for characterizing ADCs and DACs. It works in spite of the fact that the sine wave is non-linear -- which just underscores the point that the linear bit of LLS is linearity in the parameters, not the fit function. It works quite a bit faster than the Levenberg-Marquardt nonlinear least squares routine.

 

	# Least squares fitting
	# IEEE algorithm version for known frequency
	from __future__ import division
	import numpy as np
	from numpy import pi
	from numpy import fft
	from numpy import linalg as lg
	import pylab as plt
	import matplotlib.gridspec as gridspec
	## 3-parameter sine fit, with linear least squares
	dt = 1e-2
	T = np.random.uniform(3,15) # record length
	A = 1.0 # sine amplitude
	noise_amplitude = np.random.uniform(0,20)
	f = np.random.uniform(50,80) # sine frequency
	N = int(T/dt) + 1
	t = np.arange(0,T,dt)
	phase = np.random.uniform(0,2*pi)
	noise = np.random.normal(0,noise_amplitude,N) # noise vector
	y = A * np.sin(2*pi*f*t + phase) + noise
	y = np.matrix(y).T
	D = np.matrix( [np.cos(2*pi*f*t), np.sin(2*pi*f*t), np.ones(N)] ).T
	x_opt = lg.inv(D.T * D) * D.T * y # parameter vector
	a,b,c = x_opt # = A * np.sin(phase), A * np.cos(phase), np.average(y)
	residuals = y - D * x_opt
	S = lg.norm(residuals)**2
	x_cov = np.asscalar( S/(N - len(x_opt)) ) * lg.inv(D.T * D) # covariance matrix
	phi = np.arctan2(a,b)%pi # phase estimate
	delta_a, delta_b = np.sqrt(x_cov[0,0]), np.sqrt(x_cov[1,1])
	delta_phi = np.cos(phi)**2 * np.abs(a/b) * np.sqrt( (delta_a/a)**2 + (delta_b/b)**2 )
	def corr(x,y): return np.dot(x,y)/np.sqrt( np.dot(x,x) * np.dot(y,y) )
	def correlation_length(x):
	l = 0
	while abs(corr( x,np.roll(x,l) )) > 0.5: l = l+1
	return l
	N_eff = N/correlation_length(noise)
	delta_phi_estimated = np.sqrt(2) * np.std(residuals)/np.sqrt( (a**2 + b**2) * N_eff )
	phase_error = phi%pi - phase%pi
	print "SNR = ", round( lg.norm(x_opt[:-1])/np.std(residuals) ,3)
	print
	print "phase error = ",round(phase_error,3), "("+`round(delta_phi,3)`+")", "(fit)"
	print "predicted uncertainty = ", "("+`round(delta_phi_estimated,3)`+")"
	residuals, y, y_reconstructed = np.array(residuals), np.array(y), np.array( D * x_opt )
	## Plotting
	fig = plt.figure()
	gs = gridspec.GridSpec(3,6)
	ax1 = plt.subplot(gs[0:2,:-1])
	ax2 = plt.subplot(gs[2:,:-1])
	ax3 = plt.subplot(gs[2:,-1])
	ax3.xaxis.set_visible(False)
	ax3.yaxis.set_visible(False)
	ax1.plot(t,y,'b',lw=2)
	ax1.plot(t,y_reconstructed,'g',lw=2)
	ax1.set_ylabel("y")
	ax1.grid(axis='y')
	ax2.plot(t,residuals, 'g')
	ax2.set_ylabel("residuals")
	ax2.set_xlabel("Time, t [s]")
	n, bins, patches = ax3.hist(residuals, 50, normed=True, \
	orientation='horizontal', facecolor='green', \
	histtype='stepfilled')
	plt.show()

view raw sine_fitting - known frequency.py hosted with ❤ by GitHub

Exit fluff. Enter IEP !

After dinking around with getting iPython notebooks to work without trouble on my linux laptop, I've finally hit the limit of my patience for the sluggishness and foppish foofery of iPython notebooks. I tried to revert to my old gedit + Dreampie combination, but having been spoiled by executing code inline, I looked around for other alternatives. And found an amazing option: IEP 

The features I like include

• syntax highlighting and the simple, uncluttered interface
• organized file structure for projects
• easy integration with matplotlib (and not having to switch back and forth with %pylab tags to make code run on different machines)
• being able to run selected code snippets with just Alt + Enter !

One word for it. Liberating! 

Arduino gets its Due

After a few years away from Arduinos, I started thinking about them a few days ago as a possible solution for a laser feedback controller project. In the past, I have been unimpressed with the Arduino's ADC, and with its outputs: either PWMs that need to be post-filtered, or klugey DAC chips. And so I was delighted to see the specs of the new(ish) Arduino Due: 12 bit analog ADCs sampling (as fast as ~50 kHz), and real 12 bit DACs right on the board. Finally, something with enough juice to do useful things in the lab!

To celebrate, here's a simple arbitrary waveform generator, taking advantage of the Due Timer library.

	// Arduino Due simple synthesizer
	// arbitrary waveforms on DAC0
	// using Timer library to run off internal clock
	#include <DueTimer.h>
	int led = 13;
	int sensorPin = A0;
	int triggerPin = 4;
	int outputValue = 0;
	int s = 0;
	long samplePeriod = 5; // us. Needs to be >= 4 us
	int *waveform;
	int sine20[] = {2047, 2680, 3250, 3703, 3994, 4095, 3994, 3703, 3250, 2680, 2047, 1414, 844, 391, 100, 0, 100, 391, 844, 1414, -1};
	int cosine100[] = {4095, 4090, 4078, 4058, 4030, 3994, 3951, 3900, 3841, 3776, 3703, 3625, 3540, 3449, 3352, 3250, 3144, 3033, 2919, 2801, 2680, 2556, 2431, 2304, 2176, 2047, 1918, 1790, 1663, 1538, 1414, 1293, 1175, 1061, 950, 844, 742, 645, 554, 469, 391, 318, 253, 194, 143, 100, 64, 36, 16, 4, 0, 4, 16, 36, 64, 100, 143, 194, 253, 318, 391, 469, 554, 645, 742, 844, 950, 1061, 1175, 1293, 1414, 1538, 1663, 1790, 1918, 2047, 2176, 2304, 2431, 2556, 2680, 2801, 2919, 3033, 3144, 3250, 3352, 3449, 3540, 3625, 3703, 3776, 3841, 3900, 3951, 3994, 4030, 4058, 4078, 4090,-1};
	int sine100[] = {2047, 2176, 2304, 2431, 2556, 2680, 2801, 2919, 3033, 3144, 3250, 3352, 3449, 3540, 3625, 3703, 3776, 3841, 3900, 3951, 3994, 4030, 4058, 4078, 4090, 4095, 4090, 4078, 4058, 4030, 3994, 3951, 3900, 3841, 3776, 3703, 3625, 3540, 3449, 3352, 3250, 3144, 3033, 2919, 2801, 2680, 2556, 2431, 2304, 2176, 2047, 1918, 1790, 1663, 1538, 1414, 1293, 1175, 1061, 950, 844, 742, 645, 554, 469, 391, 318, 253, 194, 143, 100, 64, 36, 16, 4, 0, 4, 16, 36, 64, 100, 143, 194, 253, 318, 391, 469, 554, 645, 742, 844, 950, 1061, 1175, 1293, 1414, 1538, 1663, 1790, 1918,-1};
	int square[] = {0,4095,-1};
	int sawtooth512[] = {0, 8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 88, 96, 104, 112, 120, 128, 136, 144, 152, 160, 168, 176, 184, 192, 200, 208, 216, 224, 232, 240, 248, 256, 264, 272, 280, 288, 296, 304, 312, 320, 328, 336, 344, 352, 360, 368, 376, 384, 392, 400, 408, 416, 424, 432, 440, 448, 456, 464, 472, 480, 488, 496, 504, 512, 520, 528, 536, 544, 552, 560, 568, 576, 584, 592, 600, 608, 616, 624, 632, 640, 648, 656, 664, 672, 680, 688, 696, 704, 712, 720, 728, 736, 744, 752, 760, 768, 776, 784, 792, 800, 808, 816, 824, 832, 840, 848, 856, 864, 872, 880, 888, 896, 904, 912, 920, 928, 936, 944, 952, 960, 968, 976, 984, 992, 1000, 1008, 1016, 1024, 1032, 1040, 1048, 1056, 1064, 1072, 1080, 1088, 1096, 1104, 1112, 1120, 1128, 1136, 1144, 1152, 1160, 1168, 1176, 1184, 1192, 1200, 1208, 1216, 1224, 1232, 1240, 1248, 1256, 1264, 1272, 1280, 1288, 1296, 1304, 1312, 1320, 1328, 1336, 1344, 1352, 1360, 1368, 1376, 1384, 1392, 1400, 1408, 1416, 1424, 1432, 1440, 1448, 1456, 1464, 1472, 1480, 1488, 1496, 1504, 1512, 1520, 1528, 1536, 1544, 1552, 1560, 1568, 1576, 1584, 1592, 1600, 1608, 1616, 1624, 1632, 1640, 1648, 1656, 1664, 1672, 1680, 1688, 1696, 1704, 1712, 1720, 1728, 1736, 1744, 1752, 1760, 1768, 1776, 1784, 1792, 1800, 1808, 1816, 1824, 1832, 1840, 1848, 1856, 1864, 1872, 1880, 1888, 1896, 1904, 1912, 1920, 1928, 1936, 1944, 1952, 1960, 1968, 1976, 1984, 1992, 2000, 2008, 2016, 2024, 2032, 2040, 2048, 2056, 2064, 2072, 2080, 2088, 2096, 2104, 2112, 2120, 2128, 2136, 2144, 2152, 2160, 2168, 2176, 2184, 2192, 2200, 2208, 2216, 2224, 2232, 2240, 2248, 2256, 2264, 2272, 2280, 2288, 2296, 2304, 2312, 2320, 2328, 2336, 2344, 2352, 2360, 2368, 2376, 2384, 2392, 2400, 2408, 2416, 2424, 2432, 2440, 2448, 2456, 2464, 2472, 2480, 2488, 2496, 2504, 2512, 2520, 2528, 2536, 2544, 2552, 2560, 2568, 2576, 2584, 2592, 2600, 2608, 2616, 2624, 2632, 2640, 2648, 2656, 2664, 2672, 2680, 2688, 2696, 2704, 2712, 2720, 2728, 2736, 2744, 2752, 2760, 2768, 2776, 2784, 2792, 2800, 2808, 2816, 2824, 2832, 2840, 2848, 2856, 2864, 2872, 2880, 2888, 2896, 2904, 2912, 2920, 2928, 2936, 2944, 2952, 2960, 2968, 2976, 2984, 2992, 3000, 3008, 3016, 3024, 3032, 3040, 3048, 3056, 3064, 3072, 3080, 3088, 3096, 3104, 3112, 3120, 3128, 3136, 3144, 3152, 3160, 3168, 3176, 3184, 3192, 3200, 3208, 3216, 3224, 3232, 3240, 3248, 3256, 3264, 3272, 3280, 3288, 3296, 3304, 3312, 3320, 3328, 3336, 3344, 3352, 3360, 3368, 3376, 3384, 3392, 3400, 3408, 3416, 3424, 3432, 3440, 3448, 3456, 3464, 3472, 3480, 3488, 3496, 3504, 3512, 3520, 3528, 3536, 3544, 3552, 3560, 3568, 3576, 3584, 3592, 3600, 3608, 3616, 3624, 3632, 3640, 3648, 3656, 3664, 3672, 3680, 3688, 3696, 3704, 3712, 3720, 3728, 3736, 3744, 3752, 3760, 3768, 3776, 3784, 3792, 3800, 3808, 3816, 3824, 3832, 3840, 3848, 3856, 3864, 3872, 3880, 3888, 3896, 3904, 3912, 3920, 3928, 3936, 3944, 3952, 3960, 3968, 3976, 3984, 3992, 4000, 4008, 4016, 4024, 4032, 4040, 4048, 4056, 4064, 4072, 4080, 40880, 8, 16, 24, 32, 40, 48, 56, 64, 72, 80, 88, 96, 104, 112, 120, 128, 136, 144, 152, 160, 168, 176, 184, 192, 200, 208, 216, 224, 232, 240, 248, 256, 264, 272, 280, 288, 296, 304, 312, 320, 328, 336, 344, 352, 360, 368, 376, 384, 392, 400, 408, 416, 424, 432, 440, 448, 456, 464, 472, 480, 488, 496, 504, 512, 520, 528, 536, 544, 552, 560, 568, 576, 584, 592, 600, 608, 616, 624, 632, 640, 648, 656, 664, 672, 680, 688, 696, 704, 712, 720, 728, 736, 744, 752, 760, 768, 776, 784, 792, 800, 808, 816, 824, 832, 840, 848, 856, 864, 872, 880, 888, 896, 904, 912, 920, 928, 936, 944, 952, 960, 968, 976, 984, 992, 1000, 1008, 1016, 1024, 1032, 1040, 1048, 1056, 1064, 1072, 1080, 1088, 1096, 1104, 1112, 1120, 1128, 1136, 1144, 1152, 1160, 1168, 1176, 1184, 1192, 1200, 1208, 1216, 1224, 1232, 1240, 1248, 1256, 1264, 1272, 1280, 1288, 1296, 1304, 1312, 1320, 1328, 1336, 1344, 1352, 1360, 1368, 1376, 1384, 1392, 1400, 1408, 1416, 1424, 1432, 1440, 1448, 1456, 1464, 1472, 1480, 1488, 1496, 1504, 1512, 1520, 1528, 1536, 1544, 1552, 1560, 1568, 1576, 1584, 1592, 1600, 1608, 1616, 1624, 1632, 1640, 1648, 1656, 1664, 1672, 1680, 1688, 1696, 1704, 1712, 1720, 1728, 1736, 1744, 1752, 1760, 1768, 1776, 1784, 1792, 1800, 1808, 1816, 1824, 1832, 1840, 1848, 1856, 1864, 1872, 1880, 1888, 1896, 1904, 1912, 1920, 1928, 1936, 1944, 1952, 1960, 1968, 1976, 1984, 1992, 2000, 2008, 2016, 2024, 2032, 2040, 2048, 2056, 2064, 2072, 2080, 2088, 2096, 2104, 2112, 2120, 2128, 2136, 2144, 2152, 2160, 2168, 2176, 2184, 2192, 2200, 2208, 2216, 2224, 2232, 2240, 2248, 2256, 2264, 2272, 2280, 2288, 2296, 2304, 2312, 2320, 2328, 2336, 2344, 2352, 2360, 2368, 2376, 2384, 2392, 2400, 2408, 2416, 2424, 2432, 2440, 2448, 2456, 2464, 2472, 2480, 2488, 2496, 2504, 2512, 2520, 2528, 2536, 2544, 2552, 2560, 2568, 2576, 2584, 2592, 2600, 2608, 2616, 2624, 2632, 2640, 2648, 2656, 2664, 2672, 2680, 2688, 2696, 2704, 2712, 2720, 2728, 2736, 2744, 2752, 2760, 2768, 2776, 2784, 2792, 2800, 2808, 2816, 2824, 2832, 2840, 2848, 2856, 2864, 2872, 2880, 2888, 2896, 2904, 2912, 2920, 2928, 2936, 2944, 2952, 2960, 2968, 2976, 2984, 2992, 3000, 3008, 3016, 3024, 3032, 3040, 3048, 3056, 3064, 3072, 3080, 3088, 3096, 3104, 3112, 3120, 3128, 3136, 3144, 3152, 3160, 3168, 3176, 3184, 3192, 3200, 3208, 3216, 3224, 3232, 3240, 3248, 3256, 3264, 3272, 3280, 3288, 3296, 3304, 3312, 3320, 3328, 3336, 3344, 3352, 3360, 3368, 3376, 3384, 3392, 3400, 3408, 3416, 3424, 3432, 3440, 3448, 3456, 3464, 3472, 3480, 3488, 3496, 3504, 3512, 3520, 3528, 3536, 3544, 3552, 3560, 3568, 3576, 3584, 3592, 3600, 3608, 3616, 3624, 3632, 3640, 3648, 3656, 3664, 3672, 3680, 3688, 3696, 3704, 3712, 3720, 3728, 3736, 3744, 3752, 3760, 3768, 3776, 3784, 3792, 3800, 3808, 3816, 3824, 3832, 3840, 3848, 3856, 3864, 3872, 3880, 3888, 3896, 3904, 3912, 3920, 3928, 3936, 3944, 3952, 3960, 3968, 3976, 3984, 3992, 4000, 4008, 4016, 4024, 4032, 4040, 4048, 4056, 4064, 4072, 4080, 4088,-1};
	int random64[] = {1964, 1366, 646, 1394, 1480, 1632, 1456, 1356, 1257, 1953, 862, 384, 2021, 1974, 45, 37, 1844, 1222, 1435, 77, 1061, 1128, 376, 160, 1198, 365, 1315, 658, 1837, 1800, 656, 239, 691, 823, 446, 355, 1870, 840, 354, 1979, 1830, 1727, 1029, 1642, 1453, 344, 2, 1343, 471, 1227, 803, 1273, 311, 1191, 1919, 1513, 1365, 536, 338, 467, 1175, 237, 1301, 1375, -1};
	void setup() {
	// Pick a waveform
	waveform = sine20;
	analogWriteResolution(12);
	Timer0.attachInterrupt(wave).start(samplePeriod);
	// setting time intervals is more reliable than frequency, because of integer math
	// Timer1.attachInterrupt(trigger).start(20*samplePeriod);
	}
	void wave(){
	outputValue = *(waveform + s);
	if (outputValue==-1){ // wrap around and repeat the waveform
	s = 0;
	outputValue = *(waveform + s);
	}
	analogWrite(DAC0, outputValue); // write the selected waveform on DAC1
	s++;
	}
	void trigger(){
	digitalWrite(triggerPin, HIGH);
	delayMicroseconds(50);
	digitalWrite(triggerPin, LOW);
	}
	void loop(){}

view raw sine_test.ino hosted with ❤ by GitHub

Impedance of transmission lines

This code calculates the characteristic impedance of transmission lines of specified cross section, by calculating the electric potential through relaxation on a grid. You'll need to be able to run blitz from scipy.weave.

Here is an example of what it can do.

	# Transmission line characteristic impedance
	# Amar, 12 July 2012
	# version 3
	# All lengths in cm, E-fields in V/cm, time in ns
	import pylab as plt
	import numpy as np
	from numpy import pi
	from scipy.weave import blitz
	import cPickle, time
	eta0 = 376.730313461 # Ohms
	"""
	# Blitz test:
	a,b = 1,0
	print blitz("b = a+a")
	"""
	################### Computation grid ##########################
	size = 50.0 # mm
	numberOfGridBits = 10
	numberOfGridPoints = 2**numberOfGridBits
	dx = dy = 2*size/numberOfGridPoints
	x = np.arange(-size/2.,size/2.,dx)
	y = np.arange(-size/2.,size/2.,dy)
	X,Y = np.meshgrid(x,y)
	N = len(x)
	################## Electrodes #################################
	class CircularElectrode:
	def __init__(self,x0,y0,D,V0,innie=True):
	self.x0, self.y0, self.V0 = x0,y0,V0
	self.D = D
	if innie: self.mask = ( ((X-self.x0)**2 + (Y-self.y0)**2) <= self.D**2/4. )
	else: self.mask = ( ((X-self.x0)**2 + (Y-self.y0)**2) >= self.D**2/4. )
	#electrodes = [ CircularElectrode(0,0,.24425,V0,innie=True), CircularElectrode(0,0,.5625,0,innie=False) ] # GR900
	class RectangularElectrode:
	def __init__(self,x0,y0,w,h,V0,innie=True):
	self.x0, self.y0, self.V0 = x0,y0,V0
	self.w,self.h = w,h
	if innie: self.mask = np.logical_and( np.abs(X-self.x0) <= self.w/2., np.abs(Y-self.y0) <= self.h/2. )
	else: self.mask = np.logical_or( np.abs(X-self.x0) >= self.w/2., np.abs(Y-self.y0) >= self.h/2. )
	# Note this subtleness between logical_and vs. logical_or for innies vs. outies
	################# Calculation routine ########################
	expr = "V[step:-step, step:-step] = ((V[0:-2*step, step:-step] + V[2*step:, step:-step]) + (V[step:-step,0:-2*step] + V[step:-step, 2*step:]))/4.0"
	expr2 = "V = V_mask_inv*V + V_mask*V_electrode"
	def calculateVoltage(electrodes, V0=1.0):
	"""Calculation of voltage for an electrode configuration"""
	V_electrode = np.zeros(X.shape)
	for electrode in electrodes: V_electrode[electrode.mask] = electrode.V0
	V_mask = sum([electrode.mask for electrode in electrodes])
	V_mask_inv = np.logical_not(V_mask)
	V = V_mask_inv*np.average(V_electrode) # initialized
	print "Preconditioning ...",
	for i in xrange(4,numberOfGridBits):
	step = 2**numberOfGridBits/2**i
	for j in xrange(int(5e1)):
	blitz(expr, check_size=0)
	blitz(expr2,check_size=0)
	print "done"
	step = 1.0
	for i in xrange(int(4e3)):
	blitz(expr, check_size=0)
	blitz(expr2,check_size=0)
	V_unref = V.copy()
	for i in xrange(int(1e3)):
	blitz(expr, check_size=0)
	blitz(expr2,check_size=0)
	V_ref = V.copy()
	fractionalChange = np.sqrt( np.sum((V_ref - V_unref)**2)/np.sum(V_unref**2) )
	# Electric field
	ESquared = np.zeros(X.shape)
	ESquared[1:-1, 1:-1] = ((V[2:, 1:-1] - V[0:-2, 1:-1])/(2*dx))**2 + ((V[1:-1, 2:] - V[1:-1,0:-2])/(2*dy))**2
	# Characteristic impedance
	Z0 = eta0/( np.sum(ESquared)*dx*dy ) * V0**2
	return Z0, V, (fractionalChange, V_ref-V_unref) # characteristic impedance, Ohms
	##################### Test problem #############################
	V0 = 1.0
	outerConductor = CircularElectrode(0,0,46.06,0,innie=False)
	innerConductor = CircularElectrode(0,0,20.00,V0,innie=True)
	electrodes = [ outerConductor, innerConductor ]
	Z0,V,convergence = calculateVoltage(electrodes,V0)
	print "Z0 = ", round(Z0,2),
	print "converged to ", convergence[0]
	V0 = 2.0
	outerConductor = RectangularElectrode(0,0,49.6,49.6,0,innie=False)
	innerConductor = RectangularElectrode(0,0,19.8,19.8,V0,innie=True)
	electrodes = [ outerConductor, innerConductor ]
	Z0,V,convergence = calculateVoltage(electrodes,V0)
	print "Z0 = ", round(Z0,2),
	print "converged to ", convergence[0]
	######### Plots ##############################
	# Electric field meshgrid
	Ex,Ey = np.zeros(X.shape),np.zeros(X.shape)
	Ex[1:-1, 1:-1] = (V[2:, 1:-1] - V[0:-2, 1:-1])/(2*dx)
	Ey[1:-1, 1:-1] = (V[1:-1, 2:] - V[1:-1,0:-2])/(2*dy)
	E = np.sqrt( Ex**2 + Ey**2 )
	plt.figure()
	plt.title(r"Voltage (V)")
	plt.xlabel("x (mm)")
	plt.ylabel("y (mm)")
	C1 = plt.contourf(X,Y,V,30)
	plt.colorbar(C1)
	p.figure()
	p.title(r"Electric field, $\mathcal{E}_x$ (V/cm)")
	p.xlabel("x (mm)")
	p.ylabel("y (mm)")
	C2 = p.contourf(X,Y,E,30)
	p.colorbar(C2)
	##################### Batch problem setup ######################
	outerConductor = RectangularElectrode(0,0,30,50,0,innie=False)
	V0 = 1.0
	w = np.arange(5,25,2)
	h = np.arange(1,15,2)
	W,H = np.meshgrid(w,h)
	Z = np.zeros(W.shape)
	for i in xrange(len(w)):
	for j in xrange(len(h)):
	innerConductor = RectangularElectrode(0,-15,w[i],h[j],V0,innie=True)
	electrodes = [ outerConductor, innerConductor ]
	Z0,V,junk = calculateVoltage(electrodes)
	Z[j,i] = Z0
	print w[i],h[j],round(Z0)
	p.figure()
	p.title(r"Characteristic impedance ($\Omega$)")
	p.xlabel("Center conductor width, w (mm)")
	p.ylabel("Center conductor height, h (mm)")
	C = p.contourf(W,H,Z,30)
	p.colorbar(C)
	p.figure()
	p.title("Characteristic impedance ($\Omega$)\n"+"Inner conductor @ "+`(innerConductor.x0,innerConductor.y0)`)
	p.xlabel("Center conductor width, w (mm)")
	p.ylabel("Center conductor height, h (mm)")
	C = p.contour(W,H,Z,np.arange(0,100,5), linewidth=10, colors='b')
	p.clabel(C, inline=1, fontsize=12, linewidth=5)

view raw characteristicImpedanceTEM.py hosted with ❤ by GitHub