What is this about?

Solving electronic circuits analytically can be a pain, if not impossible with complex constellations. So therefore, simulation is the only real alternative. For such purposes, programs like Spice and its many derivates have been developed. The original Berkley Spice has the disadvantage that the code is quite old and does not compile cleanly on new Linux machines. Commercial Spice variants are powerful and easy to use, but usually just run on Windows. Sure, there is Wine - but are there OpenSource alternatives? Yes. gEDA and GnuCap are two programs, one which is for design of electronic ciruits (doing EDA, Electronic Design Automation, as the name suggests) and the other which is for simulating circuits.

Why are you writing this?

Basically, when I started working with those programs and got a little into them - I instantly fell in love. They are powerful and pretty, easy to use (sure, there's a learning curve - but there's a learning curve virtually with anything) and very customizable. I was capable of writing pretty cool scripts within a few days with the help of one other GnuCap tutorial and the really nice and competent people on the GnuCap and gEDA mailing lists.

Alright, what are we trying to do?

Basically, we'll first start out with some very simple transient analysis and then try do do some filter analysis and frequency response.

What software do I need?

Well, you need gEDA from here, GnuCap from here and gwave from here. They're all free software and wonderful tools. It will not hurt if you read the documentation, especially of GnuCap (which has a bit of a learning curve) - it is available as a wonderfully verbose PDF-file. As I was looking for assistance the mailing lists of the projects were also of great help - although the questions I asked there will also be answered in this little tutorial :-)

Getting started...

We start out with a very simple simulation: Have a square wave of 50Hz with 50% duty cycle. This means the wave is 0.01 seconds on and 0.01 seconds off and the total period is 0.02 seconds. This signal we will push trough a RC-lowpass-filter and after that through a simple diode - see what happens. We expect the rising and falling edges of the waveform to be deformed because of the RC and also some voltage drop across the diode which we should be able to measure. So lets start drawing:

Simple circuit with R, C and D in gEDA

We set the "value" attribute of R1 to 1k (which imples 1 kOhm), "value" of C1 to 1u (which implies 1 µF) and the "value" attribute of D1 to "1N4004". Then we set the "netname" attributes of the input and output ot "Vin" and "Vout" respectively. And we're already done with the drawing. Now we can export the whole thing with using the gnetlist tool (which is part of the gEDA package):

joeserver joe [~/gnucap-example-jb/01_simple]: gnetlist -g spice -o mycircuit.net mycircuit.sch
gEDA/gnetlist version 1.5.0.20080706
gEDA/gnetlist comes with ABSOLUTELY NO WARRANTY; see COPYING for more details.
This is free software, and you are welcome to redistribute it under certain
conditions; please see the COPYING file for more details.

Remember to check that your schematic has no errors using the drc2 backend.
You can do it running 'gnetlist -g drc2 your_schematic.sch -o drc_output.txt'
and seeing the contents of the file drc_output.txt.

Loading schematic [/home/joe/gnucap-example-jb/01_simple/mycircuit.sch]
joeserver joe [~/gnucap-example-jb/01_simple]:

We can check what we got there:

joeserver joe [~/gnucap-example-jb/01_simple]: cat mycircuit.net
* Spice netlister for gnetlist
R1 Vin 1 1k
C1 0 1 1u
D1 1 Vout 1N4004
.END
joeserver joe [~/gnucap-example-jb/01_simple]:

What this means is, line by line:

  1. Just a comment.
  2. Create a resistor. Call it R1. Connect it to the nets "Vin" and "1". Let it have a value of 1kOhm.
  3. Create a capacitor. Call it C1. Connect it to the net "0" (which is always impled GND) and "1". Let it have a value of 1µF.
  4. Create a diode. Call it D1. Connect it to net "1" and "Vout". Let it have the model "1N4004".
  5. That's it.

Okay, this is neat - now let's try to simulate something. Start up gnucap. I'll always highlight the prompt (which is usually not the case) to improve readability.

joeserver joe [~/gnucap-example-jb/01_simple]: gnucap
Gnucap 2008.03.24 RCS 26.76
The Gnu Circuit Analysis Package
Never trust any version less than 1.0
Copyright 1982-2007, Albert Davis
Gnucap comes with ABSOLUTELY NO WARRANTY
This is free software, and you are welcome
to redistribute it under the terms of 
the GNU General Public License, version 3 or later.
See the file "COPYING" for details.
gnucap> 

It will await you with a prompt. So we'll load the netlist in there:

gnucap> get mycircuit.net 
* Spice netlister for gnetlist 

So it has loaded the file. Let's check that it's really there:

gnucap> list
R1 ( vin 1 )  1.K
C1 ( 0 1 )  1.u
D1 ( 1 vout )  1N4004 NA( 1.)

Alright, so we now define that the 1N4004 model should just be an ordinary diode (no special model so far):

gnucap> model 1N4004 D

On some versions of GnuCap this might not work. If this is the case with your version, try:

gnucap> build
> .model 1N4004 D

And hit the return key to return to the standard GnuCap prompt. Thanks go to Brian Fuller who reported this to me on the geda-user mailing list.

If you now re-check the circuit listing, one line was added:

gnucap> list
[...]
.model 1N4004 d ( tnom=NA( 27.) is=NA( 10.f) rs=NA( 0.) n=NA( 1.) tt=NA( 0.) cjo=NA( NA) pb=NA( NA) mj=NA( 0.5) egap=NA( 1.11) xti=NA( 3.) fc=NA( 0.5))

These values describe the diode parameters, GnuCap has done some guessing for us, as we did not specify what diode exactly we wanted. Now let's tell GnuCap that the voltage source we connected is as we specified before. As the voltage source is part of the circuit description we loaded from the netfile, we will have to use the "build" command to insert it:

gnucap> build
> Vcc (Vin 0) pulse(iv=0, pv=5, period=0.02, width=0.01)
>

After that, we hit Ctrl-D or the return key to exit build mode. Now the circuit should be complete:

gnucap> list
R1 ( vin 1 )  1.K
C1 ( 0 1 )  1.u
D1 ( 1 vout )  1N4004 NA( 1.)
.model 1N4004 d ( tnom=NA( 27.) is=NA( 10.f) rs=NA( 0.) n=NA( 1.) tt=NA( 0.) cjo=NA( NA) pb=NA( NA) mj=NA( 0.5) egap=NA( 1.11) xti=NA( 3.) fc=NA( 0.5))
Vcc ( vin 0 ) pulse iv= 0. pv= 5. delay=NA( 0.) rise=NA( 0.) fall=NA( 0.) width= 0.01 period= 0.02

So we should tell GnuCap what exactly we want to analyze. Let's just analyze everything:

gnucap> print tran v(nodes)

Now we told GnuCap to output the voltage of all nodes (that's "Vin", "1" and "Vout") during transient simulation. Let's kick it off already - simulate 0.4 seconds with a 1ms timestep:

gnucap> tran 0 0.4 .001
#Time       v(1)       v(vin)     v(vout)   
 0.         5.         5.         4.8401    
 0.001      5.         5.         4.8401    
 0.002      5.         5.         4.8401    
 0.003      5.         5.         4.8401    
 0.004      5.         5.         4.8401    
 0.005      5.         5.         4.8401    
 0.006      5.         5.         4.8401    
 0.007      5.         5.         4.8401    
 0.008      5.         5.         4.8401    
 0.009      5.         5.         4.8401    
 0.01       3.3333     33.333n    3.1842    
 0.011      1.2        12.n       1.0787    
 0.012      0.43914    4.3914n    0.34669   
 0.013      0.1607     1.607n     0.098935  
 0.014      0.05881    588.09p    0.025814  
 0.015      0.021521   215.21p    0.0073179 
 0.016      0.0078757  78.757p    0.0023715 
 0.017      0.0028821  28.821p    827.01u   
 0.018      0.0010547  10.547p    297.21u   
 0.019      385.97u    3.86p      108.04u   
 0.02       179.68u    5.         50.191u   
[...]

Okay, nice as this is - Graphs are always nicer. So let's again restate the "print" command to only focus on what we want to know - that's just the voltages of "Vin" and "Vout":

gnucap> print tran v(vin) v(vout)
gnucap> tran 0 0.4 .001 >mycircuit.out

It is finished with the simulation almost instantly and has provided us with a "mycircuit.out" file - one we can open with gwave:

joeserver joe [~/gnucap-example-jb/01_simple]: gwave mycircuit.out

After pulling the "v(vin)" and "v(vout)" boxes into the black drawing area, we're rewarded with the pretty result:

The analysis of the circuit

The cursor has been placed at the 1 tau mark, where tau = R * C = 1ms. In the graph it can be seen that the voltage drop over the diode is about 0.16 Volts, therefore the the capacitor is about 3.0915V / 4.84V = 64% charged - we'd expect 63% (1 - exp(-1)) from an ideal capacitor, so the result really is accurate!

For more complex simulations or automation purposes, it is nice to not always have to edit the ".net" file, but have it autogenerated. For this purpose I made a Makefile which takes three files as input: The ".sch" schematic file, a ".hdr" header file (with circuit descriptions such as diode models or Vcc declarations, just as we added above) and a ".sim" simulation file (which describes what we want to simulate and where we want to put it). This makes it very easy to change the circuit and re-simulate, as you have to just hit "make" and it all works. So in future I will always use that style (and provide the necessary files like the Makefile, of course).

Now lets look at a voltage doubler, for example:

Voltage doubler circuit

In this circuit, R1 models the internal resistance of the voltage source Vin. An ideal voltage source could provide infinite current, thus the added resistor limits the current, in the shown schematic to 5mA (I = V / R = 5V / 1kOhm = 5mA). Also note that there is an added node called Vcc, which in the example provides a voltage source of +5V DC. This is achieved by adding to the ".hdr" file the following line:

Vcc (Vcc 0) pulse(iv=5, pv=5)

Also, should we ever need negative supply voltage, we can add that net, too:

Vcc (-Vcc 0) pulse(iv=-5, pv=-5)

This takes care of the supplies - whenever you will now connect something to "Vcc" or "-Vcc" in gEDA, it will be connected to the appropriate power supplies. Well, lets have a look at the simulation output with an internal resistance of 1kOhm:

Voltage doubler simulation with 1kOhm internal resistance

You can nicely see how the output voltage is charged up with every cycle of the pulsed input. The upper limit here is around 6.6 Volts, so its not an ideal voltage doubler - this is because of the high internal resistance we specified. Let's lower it to 0.1 Ohms and try out what the simulation shows:

Voltage doubler simulation with 0.1 Ohm internal resistance

Here we go - this looks much better. We're up to 9.3 Volts now after a few cycles.

Now for something a little bit more exciting - filters. We will first start out with the crappy, dirt-cheap, all-purpose operational amplifier LM324 (of which a model is freely available on the net, too). One warning beforehand: Be sure to use a recent version of GnuCap or you will receive a strange error message telling you there is "no subckt", because older versions of GnuCap will choke on the slash in the device name (LM324/NS). So let's look at the circuit we're going to simulate:

Active low-pass filter

First of all: We need to fetch the model for the LM324 off the net. It can be downloaded here. The package I provided already includes the file, as it may be copied and distributed (see copyright notice of National Semiconductor). However, if you plan to redistribute, keep in mind that that file is copyright protected and may not be sold in any way. If the link above is dead or in general you're looking for models, google for "xyz spice model" and try to find a ".mod" file that works.

We're doing AC analysis now, not transients. Our goal is to determine how much the signal is damped at which frequency. So therefore, we modify our Vin probe to be the AC source:

Vin (Vin 0) AC 1

Also, we have to change what we want to simulate. It is very important that before you run an AC analysis there is an operating point (OP) analysis or you will receive very strange messages about unconnected subcircuits and such. So here it just is:

.print op iter(0) v(Vout)
.op

.print ac vdb(Vout)
.ac 10 20k 10 >simplefilter.out

So we first perform the OP analysis (the manual says "This command also sets up the quiescent point for subsequent AC analysis. It is necessary to do this for nonlinear circuits. The last step in the sweep determines the quiescent point for the AC analysis.") and after that the actual analysis we want to perform. Note that we're not asking for v(Vout), but for vdb(Vout), which will yield the output voltage in Decibels relative to 1V AC. Therefore, we will get a logarithmic Y-axis (as dB is a logarithmic unit). So we're asking for an analysis from 10 Hz to 20 kHz in steps of 10 Hz. Let's see what we get:

Active low-pass filter simulation results

Well, this looks pretty good. Note that the cursor has highlighted the point at which the signal is damped approximately 3dB - this is at 4.50 kHz. If you do the math, the cutoff frequency of a active low-pass-filter should be f = 1 / (2 * Pi * R1 * C1) which in this case (R1 = 500 Ohm, C1 = 100nF) would be f = 3.18 kHz. Note, however, that we're not dealing with ideal parts, but with the crappy LM324 - the simulation is probably closer to the real world as the idealized operational amplifier is.

When your networks become more complicated and you want to try around part values in a more sophisticated manner, you'll love the flexibility of the GnuCap and gwave packages. I wrote a small Python script called "PartSweep", which will automatically vary parts values and generate a gwave-script to display those in a nice fashion. It's a bit clumsy to work with, but it's really neat once you figure out how it works. A small example will be immensely helpful. Let's say in the active low-pass filter we developed before we want to vary values of R1. Nothing is easier than:

joeserver joe [~/gnucap-example-jb/03_simplefilter]: ../PartSweep simplefilter.hdr simplefilter.net 10 10k 7 R1
Calculating R1 = 75.0...
Calculating R1 = 94.2...
Calculating R1 = 118...
Calculating R1 = 149...
Calculating R1 = 187...
Calculating R1 = 235...
Calculating R1 = 295...
Calculating R1 = 371...
Calculating R1 = 466...
Calculating R1 = 586...
Calculating R1 = 736...
Calculating R1 = 925...
Panel 0, Plot 0 (0): PartSweep_Output/ckt_tmp_0001.out 75.0
Panel 0, Plot 1 (1): PartSweep_Output/ckt_tmp_0002.out 94.2
Panel 0, Plot 2 (2): PartSweep_Output/ckt_tmp_0003.out 118
Panel 0, Plot 3 (3): PartSweep_Output/ckt_tmp_0004.out 149
Panel 0, Plot 4 (4): PartSweep_Output/ckt_tmp_0005.out 187
Panel 0, Plot 5 (5): PartSweep_Output/ckt_tmp_0006.out 235
Panel 1, Plot 0 (6): PartSweep_Output/ckt_tmp_0007.out 295
Panel 1, Plot 1 (7): PartSweep_Output/ckt_tmp_0008.out 371
Panel 1, Plot 2 (8): PartSweep_Output/ckt_tmp_0009.out 466
Panel 1, Plot 3 (9): PartSweep_Output/ckt_tmp_0010.out 586
Panel 1, Plot 4 (A): PartSweep_Output/ckt_tmp_0011.out 736
Panel 1, Plot 5 (B): PartSweep_Output/ckt_tmp_0012.out 925

It will first gather the initial value for R1, which in the example is 500 Ohms. Then it will logarithmically sweep through a area of values around that, ranging from 75 Ohms to 925 Ohms. It will run all simulations from 10 Hz to 10 kHz in 7 Hz steps, create the output and gwave file and display them immediately:

Part sweep displaying 12 graphs of frequency responses

Now you can just pick the graph that fits you best, note which one it is (ranging from "0" to "B") go back to the console and see which value was plugged in. You can iteratively develop filter circuits that way very easily.

You have also the ability to specify an upper and lower boundary of values that should be plugged in. Say you want for C1 12 graphs with values from 33pF to 330nF, in a linear fashion:

joeserver joe [~/gnucap-example-jb/03_simplefilter]: ../PartSweep simplefilter.hdr simplefilter.net 10 10k 7 C1 lin 12 33p 330n
Calculating C1 = 33.0 p...
Calculating C1 = 30.0 n...
Calculating C1 = 60.0 n...
Calculating C1 = 90.0 n...
Calculating C1 = 120 n...
Calculating C1 = 150 n...
Calculating C1 = 180 n...
Calculating C1 = 210 n...
Calculating C1 = 240 n...
Calculating C1 = 270 n...
Calculating C1 = 300 n...
Calculating C1 = 330 n...
Panel 0, Plot 0 (0): PartSweep_Output/ckt_tmp_0001.out 33.0 p
Panel 0, Plot 1 (1): PartSweep_Output/ckt_tmp_0002.out 30.0 n
Panel 0, Plot 2 (2): PartSweep_Output/ckt_tmp_0003.out 60.0 n
Panel 0, Plot 3 (3): PartSweep_Output/ckt_tmp_0004.out 90.0 n
Panel 0, Plot 4 (4): PartSweep_Output/ckt_tmp_0005.out 120 n
Panel 0, Plot 5 (5): PartSweep_Output/ckt_tmp_0006.out 150 n
Panel 1, Plot 0 (6): PartSweep_Output/ckt_tmp_0007.out 180 n
Panel 1, Plot 1 (7): PartSweep_Output/ckt_tmp_0008.out 210 n
Panel 1, Plot 2 (8): PartSweep_Output/ckt_tmp_0009.out 240 n
Panel 1, Plot 3 (9): PartSweep_Output/ckt_tmp_0010.out 270 n
Panel 1, Plot 4 (A): PartSweep_Output/ckt_tmp_0011.out 300 n
Panel 1, Plot 5 (B): PartSweep_Output/ckt_tmp_0012.out 330 n

And it will yield this result:

Part sweep displaying 12 graphs of frequency responses, now varying C1

Feel free to comment this application, improve it, play around with it or do whatever you wish - as long as it's under the GPL-3 or any later version, at your choice.

Download

The whole package of files can be downloaded here. Note that the PartSweep application is released unter the terms of the GNU General Public License version 3 (or any later version, at your choice). Note also that the LM324 model is copyrighted by National Semiconductor.