Pole-Zero example ================= Giuseppe Venturini, Thu May 7, 2015 In this short example we will simulate a simple `RLC circuit `__ with the `ahkab `__ simulator. In particular, we consider a series resonant RLC circuit. If you need to refresh your knowledge on 2nd filters, you may `take a look at this page `__. The plan: what we'll do ----------------------- 0. A brief analysis of the circuit ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ This should be done with pen and paper, we'll just mention the results. The circuit is pretty simple, feel free to skip if you find it boring. 1. How to describe the circuit with ahkab ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ We'll show this: - from Python, - and briefly with a netlist deck. 2. Pole-zero analysis ~~~~~~~~~~~~~~~~~~~~~ - We will extract poles and zeros. - We'll use them to build input-output transfer function, which we evaluate. 3. AC analysis ~~~~~~~~~~~~~~ - We will run an AC analysis to evaluate numerically the transfer function. 4. Symbolic analysis ~~~~~~~~~~~~~~~~~~~~ - We'll finally run a symbolic analysis as well. - Once we have the results, we'll substitute for the real circuit values and verify both AC and PZ analysis. 5. Conclusions ~~~~~~~~~~~~~~ We will check that the three PZ, AC and Symbolic analysis match! The circuit ----------- The circuit we simulate is a very simple one: .. figure:: ../images/pz_example/rlc_series.svg :alt: RLC series circuit 0. Theory --------- Once one proves that the current flowing in the only circuit branch in the Laplace domain is given by: .. math:: I(s) = \frac{1}{L}\cdot\frac{s}{s^2 + 2\alpha\cdot s + \omega_0^2} Where: - :math:`s` is the Laplace variable, :math:`s = \sigma + j \omega`: - :math:`j` is the imaginary unit, - :math:`\omega` is the angular frequency (units rad/s). - :math:`\alpha` is known as the *Neper frequency* and it is given by :math:`R/(2L)`, - :math:`\omega_0` is the \*(undamped) resonance frequency, equal to :math:`(\sqrt{LC})^{-1}`. It's easy to show that the pass-band transfer function we consider in our circuit, :math:`V_{OUT}/V_{IN}`, has the expression: .. math:: H(s) = \frac{V_{OUT}}{V_{IN}}(s) = k_0 \cdot\frac{s}{s^2 + 2\alpha\cdot s + \omega_0^2} Where the coeffiecient :math:`k_0` has value :math:`k_0 = R/L`. Solving for poles and zeros, we get: - One zero: - :math:`z_0`, located in the origin. - Two poles, :math:`p_0` and :math:`p_1`: - :math:`p_{0,1} = - \alpha \pm \sqrt{\alpha^2 - \omega_0^2}` 1. Describe the circuit with ahkab ---------------------------------- Let's call ``ahkab`` and describe the circuit above. First we need to import ``ahkab``: .. code-block:: python # libraries we need import ahkab .. raw:: html :: Populating the interactive namespace from numpy and matplotlib .. raw:: html .. code-block:: python print "We're using ahkab %s" % ahkab.__version__ .. raw:: html :: We're using ahkab 0.16 Then we create a new circuit object titled 'RLC bandpass', which we name ``bpf`` from Band-Pass Filter: .. code-block:: python bpf = ahkab.Circuit('RLC bandpass') A circuit is made of, internally, components and nodes. For now, our ``bpf`` circuit is empty and really of not much use. We wish to define our nodes, our components, specifying their connection to the appropriate nodes and inform the circuit instance about the what we did. It sounds complicated, but it is actually very simple, also thanks to the convenience functions ``add_*()`` in the ``Circuit`` instances (`circuit documentation `__). We now add the inductor ``L1``, the capacitor ``C1``, the resistor ``R1`` and the input source ``V1``: .. code-block:: python bpf = ahkab.Circuit('RLC bandpass') bpf.add_inductor('L1', 'in', 'n1', 1e-6) bpf.add_capacitor('C1', 'n1', 'out', 2.2e-12) bpf.add_resistor('R1', 'out', bpf.gnd, 13) # we also give V1 an AC value since we wish to run an AC simulation # in the following bpf.add_vsource('V1', 'in', bpf.gnd, dc_value=1, ac_value=1) Notice that: - the nodes to which they get connected (``'in'``, ``'n1'``, ``'out'``...) are nothing but strings. If you prefer handles, you can call the ``create_node()`` method of the circuit instance ``bpf`` (`create\_node documentation `__). - Using the convenience methods ``add_*``, the nodes are not explicitly added to the circuit, but they are in fact automatically taken care of behind the hood. Now we have successfully defined our circuit object ``bpf``. Let's see what's in there and generate a netlist: .. code-block:: python print(bpf) .. raw:: html :: * RLC bandpass L1 in n1 1e-06 C1 n1 out 2.2e-12 R1 out 0 13 V1 in 0 type=vdc value=1 vac=1 The above text defines the same circuit in netlist form. It has the advantage that it's a very concise piece of text and that the syntax resembles (not perfectly yet) that of simulators such as `SPICE `__. If you prefer to run ``ahkab`` from the command line, be sure to check the `Netlist syntax doc page `__ and to add the simulation statements, which are missing above. 2. PZ analysis -------------- The analysis is set up easily by calling ``ahkab.new_pz()``. Its signature is: .. code-block:: python ahkab.new_pz(input_source=None, output_port=None, shift=0.0, MNA=None, outfile=None, x0=u'op', verbose=0) And you can find `the documentation for ahkab.new\_pz here `__. We will set: - Input source and output port, to enable the extraction of the zeros. - the input source is ``V1``, - the output port is defined between the output node ``out`` and ground node (``bpf.gnd``). - We need no linearisation, since the circuit is linear. Therefore we set ``x0`` to ``None``. - I inserted a non-zero shift in the initial calculation frequency below. You may want to fiddle a bit with this value, the algorithm internally tries to kick the working frequency away from the exact location of the zeros, since we expect a zero in the origin, we help the simulation find the zero quickly by shifting away the initial working point. .. raw:: html .. code-block:: python pza = ahkab.new_pz('V1', ('out', bpf.gnd), x0=None, shift=1e3) r = ahkab.run(bpf, pza)['pz'] The results are in the ``pz_solution`` object ``r``. It has an interface that works like a dictionary. Eg. you can do: .. code-block:: python r.keys() .. raw:: html .. code-block:: python [u'p0', u'p1', u'z0'] Check out the `documentation on pz\_solution for more `__. Let's see what we got: .. code-block:: python :emphasize-lines: 2,3 print('Singularities:') for x, _ in r: print "* %s = %+g %+gj Hz" % (x, np.real(r[x]), np.imag(r[x])) .. raw:: html :: Singularities: * p0 = -1.03451e+06 -1.07297e+08j Hz * p1 = -1.03451e+06 +1.07297e+08j Hz * z0 = -1.44751e-13 +0j Hz **Note that the results are frequencies expressed in Hz** (and *not* angular frequencies in rad/s). Graphically, we can see better where the singularities are located: .. code-block:: python figure(figsize=figsize) # plot o's for zeros and x's for poles for x, v in r: plot(np.real(v), np.imag(v), 'bo'*(x[0]=='z')+'rx'*(x[0]=='p')) # set axis limits and print some thin axes xm = 1e6 xlim(-xm*10., xm*10.) plot(xlim(), [0,0], 'k', alpha=.5, lw=.5) plot([0,0], ylim(), 'k', alpha=.5, lw=.5) # plot the distance from the origin of p0 and p1 plot([np.real(r['p0']), 0], [np.imag(r['p0']), 0], 'k--', alpha=.5) plot([np.real(r['p1']), 0], [np.imag(r['p1']), 0], 'k--', alpha=.5) # print the distance between p0 and p1 plot([np.real(r['p1']), np.real(r['p0'])], [np.imag(r['p1']), np.imag(r['p0'])], 'k-', alpha=.5, lw=.5) # label the singularities text(np.real(r['p1']), np.imag(r['p1'])*1.1, '$p_1$', ha='center', fontsize=20) text(.4e6, .4e7, '$z_0$', ha='center', fontsize=20) text(np.real(r['p0']), np.imag(r['p0'])*1.2, '$p_0$', ha='center', va='bottom', fontsize=20) xlabel('Real [Hz]'); ylabel('Imag [Hz]'); title('Singularities'); .. figure:: ../images/pz_example/singularities_plot.png :alt: Plot of the singularities returned by the PZ analysis As expected, we got two complex conjugate poles and a zero in the origin. **The resonance frequency** Let's check that indeed the (undamped) resonance frequency :math:`f_0` has the expected value from the theory. It should be: .. math:: f_0 = \frac{1}{2\pi\sqrt{LC}} Since we have little damping, :math:`f_0` is very close to the damped resonant frequency in our circuit, given by the absolute value of the imaginary part of either :math:`p_0` or :math:`p_1`. In fact, the damped resonant frequency :math:`f_d` is given by: .. math:: f_d = \frac{1}{2\pi}\sqrt{\alpha^2 -w_0^2} Since this is an example and we have Python at our fingertips, we'll compensate for the frequency pulling due to the damping anyway. That way, the example is analytically correct. :: C = 2.2e-12 L = 1e-6 f0 = 1./(2*np.pi*np.sqrt(L*C)) print 'Resonance frequency from analytic calculations: %g Hz' %f0 .. raw:: html :: Resonance frequency from analytic calculations: 1.07302e+08 Hz .. raw:: html :: alpha = (-r['p0']-r['p1'])/2 a1 = np.real(abs(r['p0'] - r['p1']))/2 f0 = np.sqrt(a1**2 - alpha**2) f0 = np.real_if_close(f0) print 'Resonance frequency from PZ analysis: %g Hz' %f0 .. raw:: html :: Resonance frequency from PZ analysis: 1.07292e+08 Hz That's alright. 3. AC analysis -------------- Let's perform an AC analysis: :: aca = ahkab.new_ac(start=1e8, stop=5e9, points=5e2, x0=None) rac = ahkab.run(bpf, aca)['ac'] Next, we use sympy to assemble the transfer functions from the singularities we got from the PZ analysis. :: import sympy sympy.init_printing() .. raw:: html :: from sympy.abc import w from sympy import I p0, p1, z0 = sympy.symbols('p0, p1, z0') k = 13/1e-6 # constant term, can be calculated to be R/L H = 13/1e-6*(I*w + z0*6.28)/(I*w +p0*6.28)/(I*w + p1*6.28) Hl = sympy.lambdify(w, H.subs({p0:r['p0'], z0:abs(r['z0']), p1:r['p1']})) We need a function to evaluate the absolute value of a transfer function in decibels. Here it is: :: def dB20(x): return 20*np.log10(x) Next we can plot :math:`|H(\omega)|` in dB and inspect the results visually. :: figure(figsize=figsize) semilogx(rac.get_x()/2/np.pi, dB20(abs(rac['vout'])), label='TF from AC analysis') semilogx(rac.get_x()/2/np.pi, dB20(abs(Hl(rac.get_x()))), 'o', ms=4, label='TF from PZ analysis') legend(); xlabel('Frequency [Hz]'); ylabel('|H(w)| [dB]'); xlim(4e7, 3e8); ylim(-50, 1); .. figure:: ../images/pz_example/plot_pz_ac.png :alt: Transfer function plot of AC and PZ simulation data 4. Symbolic analysis -------------------- Next, we setup and run a symbolic analysis. We set the input source to be ``'V1'``, in this way, ``ahkab`` will calculate all transfer functions, together with low-frequency gain, poles and zeros, with respect to *every* variable in the circuit. It is done very similarly to the previous cases: :: symba = ahkab.new_symbolic(source='V1') rs, tfs = ahkab.run(bpf, symba)['symbolic'] Notice how to the ``'symbolic'`` key corresponds a tuple of two objects: the symbolic results and the TF object that was derived from it. Let's inspect their contents: :: print(rs) .. raw:: html :: Symbolic simulation results for 'RLC bandpass' (netlist None). Run on 2015-05-07 04:24:42. I[L1] = C1*V1*s/(C1*L1*s**2 + C1*R1*s + 1.0) I[V1] = -C1*V1*s/(C1*L1*s**2 + C1*R1*s + 1.0) VIN = V1 VN1 = V1*(C1*R1*s + 1.0)/(C1*L1*s**2 + C1*R1*s + 1.0) VOUT = C1*R1*V1*s/(C1*L1*s**2 + C1*R1*s + 1.0) .. raw:: html :: print tfs .. raw:: html :: Symbolic transfer function results for 'RLC bandpass' (netlist None). Run on 2015-05-07 04:24:42. I[L1]/V1: gain: C1*s/(C1*L1*s**2 + C1*R1*s + 1.0) gain0: 0 poles: 0.5*(-C1*R1 + sqrt(C1*(C1*R1**2 - 4.0*L1)))/(C1*L1) -0.5*(C1*R1 + sqrt(C1*(C1*R1**2 - 4.0*L1)))/(C1*L1) zeros: 0 I[V1]/V1: gain: -C1*s/(C1*L1*s**2 + C1*R1*s + 1.0) gain0: 0 poles: 0.5*(-C1*R1 + sqrt(C1*(C1*R1**2 - 4.0*L1)))/(C1*L1) -0.5*(C1*R1 + sqrt(C1*(C1*R1**2 - 4.0*L1)))/(C1*L1) zeros: 0 VIN/V1: gain: 1 VN1/V1: gain: (C1*R1*s + 1.0)/(C1*L1*s**2 + C1*R1*s + 1.0) gain0: 1.00000000000000 poles: 0.5*(-C1*R1 + sqrt(C1*(C1*R1**2 - 4.0*L1)))/(C1*L1) -0.5*(C1*R1 + sqrt(C1*(C1*R1**2 - 4.0*L1)))/(C1*L1) zeros: -1/(C1*R1) VOUT/V1: gain: C1*R1*s/(C1*L1*s**2 + C1*R1*s + 1.0) gain0: 0 poles: 0.5*(-C1*R1 + sqrt(C1*(C1*R1**2 - 4.0*L1)))/(C1*L1) -0.5*(C1*R1 + sqrt(C1*(C1*R1**2 - 4.0*L1)))/(C1*L1) zeros: 0 In particular, to our transfer function corresponds: :: tfs['VOUT/V1'] .. raw:: html :: {u'gain': C1*R1*s/(C1*L1*s**2 + C1*R1*s + 1.0), u'gain0': 0, u'poles': [0.5*(-C1*R1 + sqrt(C1*(C1*R1**2 - 4.0*L1)))/(C1*L1), -0.5*(C1*R1 + sqrt(C1*(C1*R1**2 - 4.0*L1)))/(C1*L1)], u'zeros': [0]} It's easy to show the above entries are a different formulation that corresponds to the theoretical results we introduced at the beginning of this example. We'll do it graphically. First of all, let's isolate out TF: :: Hs = tfs['VOUT/V1']['gain'] Hs .. math:: \frac{C_{1} R_{1} s}{C_{1} L_{1} s^{2} + C_{1} R_{1} s + 1.0} We wish to substitute the correct circuit values to ``R1``, ``L1`` and ``C1`` to be able to evaluate numerically the results. In order to do so, the ``symbolic_solution`` class in the ``results`` module has a method named ``as_symbols`` that takes a string of space-separed symbol names and returns the ``sympy`` symbols associated with them (`symbolic\_solution.as\_symbols documentation `__). :: s, C1, R1, L1 = rs.as_symbols('s C1 R1 L1') HS = sympy.lambdify(w, Hs.subs({s:I*w, C1:2.2e-12, R1:13., L1:1e-6})) Did we get the same results, let's sat within a 1dB accuracy? :: np.allclose(dB20(abs(HS(rac.get_x()))), dB20(abs(Hl(rac.get_x()))), atol=1) .. raw:: html :: True Good. 5. Conclusions -------------- Let's take a look at PZ, AC and symbolic results together: :: figure(figsize=figsize); title('Series RLC passband: TFs compared') semilogx(rac.get_x()/2/np.pi, dB20(abs(rac['vout'])), label='TF from AC analysis') semilogx(rac.get_x()/2/np.pi, dB20(abs(Hl(rac.get_x()))), 'o', ms=4, label='TF from PZ analysis') semilogx(rac.get_x()/2/np.pi, dB20(abs(HS(rac.get_x()))), '-', lw=10, alpha=.2, label='TF from symbolic analysis') vlines(1.07297e+08, *gca().get_ylim(), alpha=.4) text(7e8/2/np.pi, -45, '$f_d = 107.297\\, \\mathrm{MHz}$', fontsize=20) legend(); xlabel('Frequency [Hz]'); ylabel('|H(w)| [dB]'); xlim(4e7, 3e8); ylim(-50, 1); .. figure:: ../images/pz_example/plot_complete.png :alt: Plot of the transfer function from the three datasets I hope this example helped show how to use `ahkab `__ and in particular how to perform PZ, AC and symbolic analysis. If it also cleared up some doubts, great! Please remember this is an experimental simulator and you may find bug... it's getting better but we're not really ready for prime time yet: please report any and all bugs you may encounter on `the issue tracker `__. This document was written with Jupiter running with a Python kernel (project formerly named IPython). You can find it here: `Jupyter/IPython `__ and you may `access the whole notebook `__, which will allow you to download and modify this example. Have fun!