 What to do when your SPICE models run out of gas: Part 5 - Embedded.com

# What to do when your SPICE models run out of gas: Part 5

Peeling is one of the truly awesome procedures in the repertoire ofsignal analysis and modeling tools available to the engineer. It isused by oil companies looking for oil deep in the ground. It has shownup in movies, with famous actors searching for bones of Jurassic-eradinosaurs. And now it is right here in your hands.

Peeling is a procedure that generates a physical profile of someembedded object by deconvolving reflections in the time domain. Whatthat means in English, what that means for you, is that you can quicklygenerate high-quality SPICE models directly from theelectricalmeasurements of connectors and similar devices.

How Peeling Works
Envision a voltage source with known output resistance (Rs) thatconnects to a second resistor (Rl) to ground. It is easy to calculatethe value of the second resistor by measuring the open circuit voltage(Vo) and the connected voltage (Vc) using the voltage-divider equation: The schematic version of the system is displayed in Figure 4.1below. This circuit is one of the most fundamental inelectronics, andthe sad truth is that if you connect the two resistors with atransmission line, the voltage divider equation only works when a lotof qualifiers surround its application.  Figure4.1. A Voltage Divider

If the load is a transmission line and the source voltage is atransient, the transient impedance of the line can be measured by thesame procedure. However, if the source impedance is a transmissionline, this procedure no longer works very well. When something drivesor is driven by a transmission line, you calculate with the reflectioncoefficient for transients, rather than the voltage-divider equation,as follows: This equation holds for all cases, whether Z1 and Z2 are discreteelements, the characteristic impedances of transmission lines, or anymixture of the two. Through use of this equation, we can extract theimpedance profile of a device by analysis of the reflected signal seenat the input to that device.

Starting with some known source impedance and given a losslessnon-uniform transmission line, it is possible to extract the impedancevalues of each segment of the line. The procedure is called the peelingalgorithm. Figure 4.2 below depicts a coaxial transmission line with variation in characteristicimpedances through the length of the line.  Figure4.2. A Co-ax of Varying Impedance

When a transient is applied to such a line, each discontinuity inimpedance produces a reflection. Some energy goes forward, somereflects back. Following the path of the forward energy, the forwardwave encounters the next discontinuity and, again, some reflects andsome goes forward. Now, following the reflected energy, it encountersthat original discontinuity and, again, some reflects and someproceeds.

A typical device has numerous small and sometimes big variations inimpedance, each causing another point of reflection. The result canbecome very complicated. A reflection diagram, sometimes called alattice diagram, can be used to help envision the situation, asillustrated in Figure 4.3 below.  Figure4.3.Signal Reflections

Here's the way peeling works. You start at the first segment, andfrom the known voltage and source impedance, calculate the reflectioncoefficient at the beginning time (t0) of the first segment.

Knowing this coefficient, you can calculate the impedance of thatfirst segment and, more importantly, the forward and reverse energy inthat first segment. Knowing these waves, you can calculate thereflection coefficient at the right side of this second segment.

That information allows you to calculate the forward and reverseenergy in the next segment, and so on. This series of calculations iscalled the peeling algorithm. Though it sounds easy, it can be a littletricky to implement.

Implementing Peeling
For an applied transient starting at the left, the information neededto calculate the reflection coefficient at the right side of anysegment is the forward and reverse energy in this segment.

The procedure starts with only the voltage at the driving port, thestep response. The calculation presumes that the device under test(DUT) was in a non-excited steady state to begin with, and so it has noinitial reflected wave in any segment of the device.

If this is true, the first change of voltage from the initialvoltage must include two pieces of information: first, the magnitude ofthe change is the sum of the initial applied wave and the reflectedwave; second, the time it took for this to return equals twice theelectrical length of the first segment. Using this information, youcalculate reflection coefficient at the right side of the firstsegment.

To continue the solution, calculate the forward and reverse wave ineach segment, one segment at time. Two vectors, Wr and Wl, representingthe energy traveling to the right and to the left in whatever is thecurrent segment, F must be defined. The number of elements in each ofthese vectors is equal to the number of time samples in the impulseresponse.

That is to say, element Wr(i) and Wl(i) are the forward and reverseenergy waves in the current segment. And i is the number of timesamples since the stimulus first arrived at this segment.

Initialize the system by setting each element except the first, inWr, to zero. Set the voltage of the first element of Wr to 1 – theapplied step is normalized to one. Initialize all elements of Wl tozero. Now, from the second element to the last, Wr, the value isinitialized according to this next equation. Here Sr(i) is the measured step response at time slot i at the portof the device under test. It is the data set that came out of your TDRmachine or that you generated from applying the IFFT to your frequencydomain measurements.

As each layer is peeled off, the first thing needed is thereflection coefficient at the right side of the segment. Note that inthis one instance, r(n) – the reflection coefficient at  the rightof segment n – equals the voltage of Wl(2). If that equivalence looks alittle like smoke and mirrors, more detail should help clarify thisstep.

Initially the DUT is idle, dead, with no forward or reflected waves.Then at some time (t0), a unit step is applied. A forward wave of unitamplitude is injected into the initial segment of the transmissionline. To easily envision this action, presume that the initial segmentis a section of a 50-ohm transmission line. Presume a reflectioncoefficient of “0.25. This is shown in Figure4.4 below.  Figure4.4. The First Reflection

The voltage measured at the first node at time t1, the sum of theforward and reflected waves, is 0.75 volt. The voltage applied is stillone volt. So there must be a negative quarter volt being reflected. Thedifference between the voltage at time zero and the voltage at time oneis the reflection coefficient. Look at it from the perspective of thecurrent.

The driver puts out the current required to make one volt at itsoutput. In reflecting 25 percent of it, the load sees 0.75 volt. Thedifference in voltage between time zero and time one shows what theforward and reverse waves must be in that time interval. Repeat theprocess for each time interval, using the complete set of voltagesmaking up the step response. The situation at the first time intervalsis depicted in Figure 4.5 below.  Figure4.5. Initializing the Vectors

So, the measured impulse response is given, and we need the forwardand reverse energy in the first segment of the transmission line. Tofind these vectors, take the differences of the applied voltagesbetween the time steps.

The voltage at t2 cannot yet be predicted, nor can the reflectioncoefficient be calculated from t1 or t2 data because a reflection fromfurther to the right may have added to this voltage.

Only the first step can guarantee that there is no impact fromfurther right, so only the first two elements can be used to calculatethe reflection coefficient. All the remaining forward and reflectedelements are calculated by taking the differences, and the system isinitialized, ready for peeling.

Peeling requires two relationships. First, you need the sum of theenergy waves approaching the right side reflecting boundary, whichequals the sum of the energy waves moving away from that boundary.Second, you need to know the reflection coefficient for waves from theleft. Name the energy components involved at the reflecting position asshown in Figure 4.6 below.  Figure4.6.Naming the Components

The reason that we assert that the sum of waves coming in must equalthe sum of waves going out is that the segment is presumed lossless.Though this presumption is never strictly true in any real system, itis close enough for use.

In a real system, the sum of all the energy reflected out of anyjunction will always be less than the sum of the energy inserted. Itwill be less by the amount of energy lost to resistance and toradiation. You must weigh the impact of that presumption to determinewhen peeling is appropriate for the particular system you aredesigning.

The conservation of energy gives one of the two required equations;the definition of the reflection coefficient gives another. Incombination, they yield all the information needed to derive all thewave elements in the next segment to the right from those in thecurrent segment.

The reflection coefficient for a wave approaching from the rightdoes not equal the reflection coefficient for a wave approaching fromthe left. When you examine the equation, you see that though the twodiffer, they differ only in sign. Trading places between Z1 and Z2 onlychanges the sign. The reflection coefficient from the left is the negative of thereflection coefficient from the right. Due to conservation of energy,what doesn't get reflected goes through. The transmission coefficientis one minus the reflection coefficient.

For me it works best to use the convention that, in all cases wherethe reflection coefficient is used, the value is as seen by a wavecoming from the left. So, with knowledge of the forward and reverseenergy left of the boundary and knowledge of the reflection coefficientat the boundary, you can calculate the forward and reverse energy atthe right side if the boundary.

These are the Wr and Wl elements in the next segment. And it isagain true that the top, or the first, reflection from the right forthis segment involves no reflected wave from further out and so can beused to calculate the reflection coefficient at the right of thissegment. That reflection reconstitutes from this segment all theinformation needed to calculate the next segment, and so on. That iswhy we call it peeling.

The result of peeling can be a list of reflection coefficients or alist of impedances, one for each time increment. All the rest of theinformation about forward and reverse waves can now be discarded. It isknown that the system starts with a 50-ohm driver source resistance, sothe next impedance is calculated from this impedance and the reflectioncoefficient, the next from that, and so on. That step reconstitutes the impedance profile of a losslesstransmission line equivalent to the original data. Here is a simpleexample. This SPICE code is for a simple lossless transmission linewith varying impedance through its length.

You can use this model in two ways. You can send a voltage step intoit and record and peel the reflected wave, or you can do a frequencysweep of it, use the IFFT to convert that to time domain, and implementthe peeling. The following examples show both. Given this SPICEcircuit: Example #1: Analysis in the TimeDomain
The first example using this code results in the graph of Figure 4.7, below. For this case,the circuit was analyzed in the time domain. The step response wascalculated by SPICE. An extremely fast edge-rate on the step results innice square corners.  Figure4.7. SPICE Step Response

The next step would be to peel that reflection profile to recoverthe original impedance profile. But that would be too easy. We'll do itlater. First, let's look at the same model in the frequency domain.

The SPICE model was frequency swept and S11 was calculated. Theobjective now is to take this calculated S11 and attempt toreconstitute the original waveform shown in Figure 4.7 above.

From that reconstituted wave, we'll apply peeling to regenerate thetransmission lines in the SPICE model. In examining this loop fromSPICE model to S11 to step-response to SPICE model, numerous issues areuncovered and insights become available.

Example #2: Creating S Parameters
At this time, only S11 is used. As an example, a lossless, non-uniform,single-ended transmission line is modeled in SPICE. A losslesstransmission line is used because the accuracy of the method decreasesas loss is introduced.

For now we go lossless and minimize these complications. The exampleis a non-uniform transmission line because the characteristics of auniform transmission line can be calculated by hand in moments, so areuninteresting. Also, note that it is the impedance that is non-uniform.The transmission medium is presumed linear. Single-ended circuitry, asopposed to differential, is used for simplicity.

As the above SPICE code shows, the DUT consists of a section of50-ohm transmission line followed by a section of 30 ohms, a section of60 ohms, and a final section of 45 ohms. S11 will be generated totallyby hand; that is, by use of only fundamental SPICE operations to makeit clear how you generate S parameters. In application, you would usewhatever facilities your version of SPICE has to simplify this step.

To generate S11, you need to drive the DUT with an AC source thathas a series resistance of 50 ohms. The far end, or the receive end, ofthe DUT is terminated with 50 ohms to ground. Analysis is in the formof an AC sweep of the DUT. The data you need to acquire is the voltageand current at the input node of the DUT.

An important point to note is that, since S parameters are complexnumbers, the complex voltage and current must be used. Most versions ofSPICE, possibly all versions, default to outputting voltage and currentmagnitudes. Magnitudes are not adequate for the calculations. Print outthe real and imaginary parts of the voltage at the input of the DUT,and the real and imaginary parts of the current through the AC source.

Looking forward, the S11 is used to calculate the impulse responsethrough the inverse fast Fourier transform (IFFT). That procedurerequires many values evenly spaced in frequency—a harmonic sequence.That spacing translates to a requirement that the AC sweep be a linearsweep. SPICE can also do others. “Many” translates to a number in therange of 2,000 to 5,000. Thus, the analysis line in the SPICE filemight look something like:

.ac lin 5000 1 2e10

You should sweep to go to quite high frequencies because therise-time achieved is determined by the highest frequency in the sweep.It is desirable to go to quite low frequencies if the DUT responds toDC.

In this example, the frequency sweep starts at one, which is a goodstarting point because it is near zero. SPICE cannot sweep from zero.Yet, a true harmonic sequence would start at zero. To get a reasonableapproximation to a harmonic sequence, start at a frequency that is lowor small compared to the frequency steps in the sequence.

One more important point to which many might say “well of course”and others will spend a week scratching their heads is why SPICE errorsout on such a simple program, as shown. To avoid this error, insert aSPICE line such as:

.option limpts=15000

Many versions of SPICE default to presuming a programming mistake ifyou attempt to print out more than a few hundred data points. You needto print out thousands. The limpts option enables you to print outthose thousands.

The numeric value you use is unimportant as long as it is bigenough. The only reason this option defaults to a small number is thatback in the early seventies when SPICE originated, printers were slowand printer time was expensive.

Generation of the Impulse and StepResponse
The next procedure is to generate the impulse response and the stepresponse from the S11 data. You can do this step with various toolssuch as Matlab or Octave. If you simply take the S11 values and runthem through IFFT, the result is a set of complex numbers. What youneed is a set of real numbers. The impulse and step responses are real,not complex.

When you take the FFT of a real function, the result is a set ofcomplex numbers from zero to some maximum frequency, and then a mirrorimage of that data set extending again as far. The S11 parametersconsist of only half of this set. The required set needs to be formedfrom the S11 data.

This is done by mirroring each data point across the maximumfrequency value. In implementing this step, note that the mirroredvalue is not precisely equal to the original value. Rather, themirrored value has the same real part and the negative of the imaginarypart. It is the complex conjugate. In Matlab code, this might be:

Se=S11;
Se(n+1:2*n)=S11(n:-1:1)';

When this is done, applying the IFFT to the Se vector can result ina real, rather than complex, impulse response. The step response is theintegral of the impulse response. Of course, to get a step responsethat looks like the step response of a TDR, it is necessary to add anintegration constant of one to each value.

However, when that has been done, except for possibly the rise time,you can expect to see the exact same step response as would be seen ifa transient analysis had been done in SPICE. If the same step responseis not seen, check the source amplitude used in generating the ACsweep.The first 100 points of the impulse response calculated accordingto the above procedures looks like this:  Figure4.8. Peeled Impulse Response

It appears in this graph that the values from 20 to 50 and from 60to 100 are all zero, but they are not. This difference will show upwhen the impulse response is integrated to yield the step response. Atthat point, it shows:  Figure4.9. Peeled Step Response

The details to the twentieth time increment of this graph areexpanded for comparability to the SPICE simulation shown above. Figure4.10 shows just the first 20 samples of Figure 4.9 above.  Figure4.10. Expanding the First 20 Samples

The SPICE time domain simulation graph is repeated in Figure 4.11, below for convenientcomparison.  Figure4.11. The Original for Comparison

Comparing the S11 derived plot to the original SPICE time-domainplot of the network, it is clear that they are similar, but therise-times are much slower in the plot we derived from the S11 data.

You might be tempted to say they are identical except for the risetimes. That wouldn't be precisely true. Careful examination shows thatregions that were totally flat in the original simulation now have avery slight slope to them.

The graph corresponds quite well to what SPICE showed for the stepresponse, but the rise time is too slow because the frequency sweepused to generate the S11 values was only swept out to 5 gigahertz. Thesame procedure with the maximum frequency extended to 50 gigahertzyields Figure 4.12 below .  Figure4.12. A Failed Attempt

Now the edges are sharper, but the flat parts are badly misshapen.What was done in SPICE was to change:

From .ac lin 2000 1 5e9
To .ac lin 2000 1 5e10

The impact of this change was to increase both the highest frequencyand the spacing between measurement points by a factor of ten.Increasing the highest frequency yielded sharper transitions, asdesired.

Increasing the spacing between data points had the undesirableimpact seen in all the low frequency components of this graph. Yes, wecould, and should, change from 2,000 data points to 20,000 data points.But what would you learn if we did it right?

As a compromise between adequate rise time and massive quantities ofdata, the next attempt uses the SPICE analysis line:

.ac lin 5000 1 2e10

This compromise, sweeping to 20 gigahertz and acquiring 5,000 datapoints, yields an improvement and serves to show how the dataparameters interact for this problem.  Figure4.13. A Better Compromise

This result looks fairly good, perhaps good enough for an example.The next step is to run the peeling algorithm and see the results.Certainly this diagram is not “as good as it can get”, but it doesexemplify the impact of trade-offs.

I have no doubt that clever engineers will make further improvementson such trade-offs by interpolation or extrapolation from data sets.Beyond just showing the procedures, these diagrams show some of theissues you encounter when you put these procedures into practice.

Figure 4.14 below shows theresult of peeling the data set of this example. Two things should beclear at this time. First, the result of this peeling is not a perfectreproduction of the original transmission line. Second, any desireddegree of perfection can be achieved by adding more data points.

Recalling that the original model was 50, 30, 60, and 45 ohms, it isclear that the procedure is converging toward that result. It is alsoclear that with this fast a rise-time, a lot of data is required to getgood results.  Figure4.14. The Impedance Approximation

It isn't necessary to start with frequency domain data to dopeeling. Frequency domain data is used in this example only to make itclear that you can do peeling whether you start in the frequency domainor the time domain. You need to understand what characteristics make agood dataset for this type of application. Certainly, if you havetime-domain data in the first place, use it.

You also have undoubtedly noticed that Figures 4.8 through 4.14 , with asingle exception, use sample numbers rather than time in the horizontalaxis. For the points that needed to be made in that section, time wasirrelevant. One point is worth mentioning.

The timing of events you observe in a reflected signal includes boththe time to the feature and the time from the feature. Thus, the widthof features in reflections appear to be twice their actual width. Thistime doubling has various implications. For example, if you observe areflection that has a 50-picosecond rise time, the actual point ofreflection saw a 25-picosecond rise time.

This example and Figure 4.14 above show an impedance profile derived from S11 data. There are otherpossibilities. Often your interest is in what signal makes it throughto a second port rather than what gets reflected back at this port.

A typical case is a two-port device in which the signal observed atport two, resulting from the signal applied at port one, is ofinterest. The signal that gets through is quantified as S21 data.Generating an equivalent S11 data set and applying the peelingalgorithm to that S11 set can model this.

In a lossless system, energy is conserved so S11 can be fullydetermined by S21. When the peeling algorithm was described, a key step was enabled bythe presumption of a lossless system. Here again the statement was madethat S11 is fully determined by S21 in a lossless system. Real systemsare never lossless.

Real systems have resistive losses, dielectric loss, perhapsradiative loss, and maybe others. In a real system, loss causes adifference between the total energy injected into a system and the sumof energy that is transmitted and reflected.

As the portion of total energy lost in the device increases, thedeviation from that predicted by this equation increases. The peelingalgorithm is useful, even in lossy systems, but you must be aware thataccuracy decreases as loss increases. Yet, in a lossless or effectivelylossless system, since S21 is totally determined by S11, if yougenerate a model that gets S11 right, it must as a consequence also getS21 right.

So, even in a lossy system where S11 cannot be derived from S21, thepeeling algorithm can be useful. In practical applications, peelingwould not be used on the bulk transmission lines across the board, butrather on connectors and similar devices. These are usually selected tobe minimal loss in the first place, so peeling works well.

Dealing with Loss
A reasonable question to ask at this time is how severe is the loss ofaccuracy that results from the presumption of losslessness. To presenta graphical answer to that question, we modeled a transmission linewith frequency-dependent loss, using the HSPICE W model.

The W model does a fairly good job of correlating to the lossesmeasured on an FR4 circuit board. Physical and electrical conditionsother than loss were similar to the lossless system that has beenanalyzed here, but lengths were made larger. The result of peeling thismodel is displayed here.

The most obvious observation is that transitions are crisp at thebeginning but become increasingly rounded with distance, and thelargest step shows some numerical instability.

The experienced signal integrity engineer would look at this curveand recognize that it represents, at about 6 inches of FR4 trace pernanosecond, more than 6 inches of lossy trace. Most objects, such asconnectors, that you might need to model in this manner are likely tobe significantly smaller than this. Also, FR4 is probably one of themost lossy materials you will have to work with.  Figure4.15. Same Profile, Very Lossy Lines

You can draw several conclusions. Foremost is that, though theprocedure is not perfect, it is going to be useful in manyapplications. A second conclusion is that, in any attempt to measureand model a device such as a connector, minimizing the loss between thetest equipment and the device being measured enhances accuracy. Thisconclusion is valid and important whether data acquired is intended forpeeling or not.

More Fourier TransformConsiderations
Many presumptions can be hidden in such developments. In doing theFourier transform and its inverse, where the data came from can make adifference. Consider the discrete Fourier transform itself. Evenlyspaced samples in one domain are translated to the complementarydomain.

As was seen above, that could turn out to be a lot of data points.Consider the frequency domain. In SPICE, a frequency sweep can be donein various ways. It is instructive to perform a wideband frequencysweep of a device, such as a resistive divider, that is simple and easyto examine. For a real simple example, try an analysis line such asthis:

.AC lin 10 40 100

This line asks SPICE for a linear sweep of 10 evenly spacedfrequencies starting at 40 hertz and ending at 100 hertz. To implementthis sweep, SPICE divides the difference (60) by nine, to achieve 10evenly spaced frequency points. The separation will be 6.6667 each.

Using that spacing, a Fourier series would also have had severalpoints at frequencies below 40. This simulation does not provide valuesfor several data points at the low end of the frequency spectrum.Therefore, as it stands, it is not suitable for use in the inverseFourier transform. It is not a harmonic sequence.

In a SPICE simulation, this seems easy to remedy: start the sweep ata lower frequency. If 40 is replaced by two, the range is large enough;there are almost enough data points, but it still is not right. What aFourier transform would have provided would have been the frequencysequence: 10, 20.

What this simulation provides is 12.9, 23.8. To actually get itright, the sweep has to start at zero. SPICE cannot do a frequencysweep that starts at zero.

In SPICE you could get it nearly right by starting at a really lowfrequency, but other cases are not so easy. Many network analyzers haveminimum frequencies that are fairly high. The same problem exists, butnot the same solution. The obvious solution to both situations is tomathematically resample the data set, interpolating the data to evenlyspaced intervals that fit the sequence: Another problem remains. Neither SPICE nor the network analyzer canprovide the frequency response at zero, but this problem need notbecome too serious an issue. The component at zero frequency has theimpact of setting the DC offset in the time domain.

Usually, the circuit conditions make it obvious what the right valueshould be. For example, if the through response of an AC coupled linkis being computed, the frequency response at zero must be zero.

In Figure 4.15 above , anoscillation in the response near a sharp transition was observed. Thisringing, called the can happen any time there is a large step in evenly spaced data points.

In fact, any time you have a step or a transition in a single timeor frequency increment, the Gibbs Phenomenon can occur. In fact, theringing will have a magnitude of about 18 percent of the step size,which is not reduced by taking more data points in the same range. But,ah the glories of math!

The Gibbs Phenomenon can be tricked largely out of existence. Ifmore data points are taken, the magnitude remains the same, but thedistance to which the ringing extends from the transition is reduced.

As long as the data is going to be resampled to place it on thecorrect frequency domain locations anyway, you might as well resampleto many data points and then use a filter to filter off the GibbsPhenomenon. If the number of data points is raised high enough, thefilter has arbitrarily small impact on the overall circuit response.

At this time, I'm not going to cover the subject of how many pointsare adequate to provide a specific level of error or what form offilter is most suitable. The Web has a lot of literature on the GibbsPhenomenon and ways of dealing with it.

For circuit design purposes, it is sufficient to know that you getgreater accuracy when you use slower calculations. The computation timeinvolved is not minutes but seconds, so the pain level is notsufficient to justify an extensive study of the trade-off.

Actually, the benefits are better than might be expected from theabove information alone. If the simulation or measured frequency datais going to be resampled anyway, it need not be evenly spaced in thefirst place.

If it doesn't have to be evenly spaced, the SPICE simulation can useone of the nonlinear sweeps to get much better frequency coverage atlower numbers of points. Summarizing all this, time is lost in thepost-processing end of the procedure, but time is saved in thesimulation end.

Next in Part 6: Mason's Rule andCrosstalk
To read Part 1, go to “Unmodelable featuresof high performance designs