# Modeling with MatLab

Dr. Dobb's Journal November 1997: Modeling with MatLab

## A flexible tool for system solutions

Mark is a design engineer for Motorola and can be contacted at [email protected]

From PC boards and networking to cellular phones and pagers, communication systems are everywhere. Because designing and implementing communication systems is a complex process, modeling has become an important tool in the communications-design arena. In this article, I'll describe how I use the matrix laboratory (MatLab) from MathWorks to model communication systems. MatLab is an interactive system for mathematical computation and visualization that integrates numerical analysis, matrix computation, signal processing, and graphics in an interactive environment. MatLab has its own programming language, which supports control structures (FOR, WHILE, IF), string manipulation, ASCII and binary file I/O, debugging, and C code generation (if you have the MatLab compiler).

### Communication System Characteristics

All communication systems include a sender of data, a means or channel by which the data travels, and a receiver that listens to the data. Figure 1 illustrates the basic communication system I'll model in this article. The channel is a Category 3 10BaseT cable commonly used in LAN applications. The sender and receiver could be computers, network hubs, or routers. The modeling process is divided into three parts:

• Modeling the signal transmitted by the sender as it appears to the receiver. This signal is determined by the transmit characteristics of the sender and the cable attenuation.
• Modeling the receiver's ability to recover from signal loss due to the cable. Typically this is accomplished with an equalizer within the receiver. The equalizer adds flat DC gain and boosts the high-frequency parts of the received signal.
• Providing the model with timing-recovery capability to sample the received signal at the appropriate time to recover the data. Most often, some sort of phase-locked loop (PLL) is involved.

The ability to quickly put together system-level models lets you test different architectures for usability. MatLab's programming language allows for quick development and testing. And since it is extensible, you can create your own functions and/or purchase toolboxes that contain previously coded functions. (Available functions include Fourier transforms, matrix mathematics, complex numbers, butterworth filters, bilinear transform, and control system functions.)

### Examining the Signal

For simplicity, each signal will consist of three parts (similar to the method used in the IEEE 802.12 LAN protocol); see Figure 2. The first part consists of low-frequency tone signals used to adapt the automatic gain control block of the receiver. Because they are low frequency, tones do not lose much amplitude as they travel over the cable. They therefore give a good approximation for the peak voltage received. The second part of the signal is a constant high-frequency preamble used to lock the PLL in the clock-recovery circuit and adjust the high-frequency boost section. The third part of the packet is the data, which is delineated by a startdata/enddata pattern.

There are numerous ways to model the channel. Here, I'll make bench measurements of the bit response and use superposition to create an analog waveform like Figure 3. For example, take 100 meters of cable and measure the signal at the receive end while the sender signal changes from a high-impedance state to a logic 1 state, then back to the high-impedance state. A high-impedance state is the output state the sender remains in when not transmitting data. It is neither a logic 1 nor logic 0. Using the Superposition theorem, the measured bit response is then put through a shift-flip-add routine where a transmitted 0 is the measured bit response (see Figure 4) multiplied by -1, and a transmitted 1 is just the measured bit response. Listing One performs this manipulation. The analog waveform that results from the shift-flip-add algorithm models the transmitted digital data at the input of the receiver.

### Modeling the Equalizer

The next part of the system to be modeled is the equalizer, which attempts to return the signal lost due to the cable. The equalizer is broken into two parts -- an automatic gain control (AGC) and a high-frequency boost (HFB) section. The AGC section controls the value of the peak voltage so that it is within the dynamic range of the circuit. In other words, the AGC applies positive or negative flat gain to the incoming signal to maintain a constant peak voltage. The HFB section adds high-frequency gain to the incoming signal, which was lost because of the cable's low-pass filter attributes.

The AGC is a constant multiplier to the input signal. Depending on the peak voltage at the output of the AGC block, the input value from the control block will be changed to add more or less gain to the AGC input signal. To reduce the effect the AGC has on the HFB section, the AGC will change its gain at a different time than the HFB section. Allowing both the AGC and HFB to change at the same time can cause adaptation loop oscillations. Loop oscillations occur when the AGC says the signal should be lower and the HFB says the signal should be higher, or vice versa. When this occurs in adaptive equalizers, it is often difficult to adapt correctly and errors in received data occur. Some second-order effects of frequency response can be added, but shouldn't be necessary for the high-level model I've designed here. Listing Two is an example of the AGC MatLab code. The algorithm examines the received tones for zero crossings. When they occur, it waits four bit times, then takes a sample. When the preamble (a section at the beginning of the received signal that is of the same frequency as the data bits) is detected, all these samples are averaged, and the resulting peak is used to scale the received signal. To detect the preamble, it is necessary to realize that if a zero crossing occurs before the sample is taken, tones are done. Tones are 8 bits long, so when a zero crossing occurs before 8 bits, you assume the preamble has started.

The second step in equalization is the addition of high-frequency boost. This is necessary because the cable acts as a low-pass filter, causing the signal to lose more and more of its high-frequency components as it travels down the cable. The loss of high-frequency components makes transitions in the data stream occur more slowly. To fix the problem, you can add a filter that compensates for the cable loss by adding enough high-frequency boost to flatten out the loss curve over the desired frequency range. Figure 5 shows the dB loss of a long and short cable over frequency. You can see that the short cable does not have as much high-frequency loss as the long cable. An ideal filter would add enough gain to make the curve flat across the band of interest: 100 KHz-20 MHz, in this case. The responses in Figure 5 also include magnetics that add a 5th-order Butterworth with a cutoff at 20 MHz. These magnetics are part of the channel specification. It is a good idea to have linear phase in the equalizer filter to reduce the variance in the group delay. Group delay is a measure of the dispersion of different frequencies passing through a filter. For this model, group delay will not be a concern because the model that will be used for timing recovery is ideal.

I use the MatLab function FILTER() to model the filter. The first thing to worry about is the order of the filter and how many poles and zeros it should have. Based on experience, I'll use a 2-zero/3-pole filter. The selection of poles and zeros determines the amount of boost and the phase characteristic of the filter. For this article, I'll assume there is no noise in the channel. This lets me fix the values of the poles and zeros, and just add high-frequency boost to the signal when needed. In most communication systems, the channel is variable depending on the length and temperature; therefore, an adaptive equalizer is necessary to receive the data without errors. In this model, use of the preamble voltage value will determine if boost is necessary.

Picking the poles and zeros becomes a matter of flattening out the frequency response of the channel. Essentially, this means multiplying by one over the transfer function of the cable. For this case, assume that the transfer function is unknown, thus a brute-force method of trial and error can be used to vary all the poles and zeros while measuring the jitter/noise margin output. The set of poles and zeros that have the lowest jitter and the highest noise margin are selected.

At this point, I need to introduce the concept of an "eye diagram," a graph commonly used in communication systems to view the noise margin ("eye opening") at the sample point and the jitter at zero crossings. Figure 6 shows two eye diagrams, created by sending data down a 100-meter cable and overlapping the received signal every 100 nanoseconds. The idea behind this graph is that jitter can be measured by looking at the width of the zero crossings, and noise margin can be measured by looking at the height of the lowest signal, with respect to the maximum signal, at the sample point. The sample point is generally midway between two zero crossings. Looking at Figure 6, it is apparent that high-frequency boost is necessary to open the eye. Once the poles and zeros are determined, Example 1 completes the use of FILTER(), illustrating how easy it is to implement a filter using MatLab. The first line changes the form of the zeros and poles to a more-common form: the coefficients of a transfer function H(z) (see Figure 7). The second line calls the built-in MatLab function with the desired coefficients and the input data stream. This illustrates why MatLab is useful in system-level modeling -- the FILTER() function is already available and need not be implemented by the user. Listing Three is the MatLab code for the HFB section of the equalizer.

### Clock-Recovery Modeling

Now that the AGC and HFB section have been modeled, it is time to turn to clock-recovery modeling. This always seems to be one of the most frustrating parts of the model because there are many different ways to implement it, and phase-locked loops are forever causing headaches for designers. Many clock-recovery circuits use a PLL approach to recover the timing of the received signal and the data. While the idea of a PLL is conceptually simple, the implementation of the system is challenging. To the modeler, it comes down to accuracy versus simulation time: The more accuracy you need, the longer the simulation times need to be. A PLL's job is to lock the phase and frequency of a reference-frequency clock to that of an incoming signal, allowing data to be recovered without error. A typical way to do this is to use some preamble at the beginning of the received signal that is of the same frequency as the data bits. In the model presented, we indeed assume the preamble exists, and that it produces zero crossings at the data rate.

The approach I used to model the PLL required making a few assumptions. First, I assume that there is no drift in the edges during data transmission. In a real system, this is not the case, but if the number of data bits received between preambles is limited, this should not be a problem. The algorithm of the PLL in this system looks for zero crossings in the preamble. When they are found, it measures the time between the current crossing and the last crossing. When data starts, it takes the average of the zero-crossing times and uses this for the clock period of the data-recovery clock. The data-recovery clock is then synchronized to the transition in the startdata pattern. This way, the PLL model has accomplished both a frequency and phase lock to the data without going through the time-consuming complete PLL simulation. Listing Four is the MatLab code for the PLL implementation.

Another interesting approach is to use a Digital PLL that checks the distance between every edge in the received signal and the reference-clock edges during the data and the preamble. If the distance between edges is not zero, the PLL then changes to a different phase of the reference clock. For example, if the edge of the received signal occurs before the edge of the reference clock, then it picks a phase of the reference clock that makes the next cycle shorter. In essence, this speeds up the reference clock for a cycle, so that the next edges occur at the same time. A good patent on this subject is Patent No. 5,185,768 by Ferraiolo et al. from IBM, which describes a Digital PLL architecture similar to that described here.

### Testing the Model

What do you accomplish by modeling gain control, high-frequency boost, and a timing-recovery circuit? A simple test can demonstrate the usefulness of system-level modeling: Send the packet generated with the bit response from a long cable through the system model without the high-frequency boost. When this is done, errors occur between the data sent and data recovered. These errors are due, in part, to timing errors and amplitude errors. The system model can show these errors. The earlier that errors in the design are caught, the less costly they are to fix. System-level modeling helps to check errors before designers even start with the circuits. Listings Five and Six implement the bit error-rate checking and run the simulation.

### Conclusion

As the code listings illustrate, MatLab provides an easy-to-understand programming language that can be used for system modeling. With its interpretive language, quick changes can be made from the command line. The many built-in functions aid in the development of system models, and toolbox extensions can be added to increase the capabilities of MatLab programs. One new feature is the MatLab compiler, which takes the interpretive MatLab code, creates C versions of it, then compiles. These compiled versions of code execute many times faster, once again allowing the designer to view more architectures or complex simulations than before. By allowing the freedom to explore, MatLab is indeed a system-level programming tool with extensions for almost any system imaginable.

The Math Works
24 Prime Park Way
Natick, MA 01760-1500
508-647-7000
http://www.mathworks.com/

#### Listing One

```function [wave]=ddjshiftflipadd(bits, impulse, samplesperbit)% This function does what is called the Shift-Flip-Add of the impulse.
% It creates a convolution of the impulse and a bit pattern.

</p>
% Initialize
wavelength=ceil(length(bits)*samplesperbit+length(impulse));
wave(wavelength)=0;
position=0;               % This is the real position
rpos=0;                   % This is the rounded position used for array indexs

</p>
% Do for all bits in the input stream....
for i=1:length(bits)
if (bits(i)==1)
for j=1:length(impulse)
wave(rpos+j)=wave(rpos+j)+impulse(j);
end
end
% flip impulse for a zero
if (bits(i)==0)
for j=1:length(impulse)
wave(rpos+j)=wave(rpos+j)-impulse(j);
end
end
position=position+samplesperbit;
rpos=round(position);
end
% Because of bit response sampling and noise in the setup of the bench it is
% necessary to truncate the end of the resultant convolution.
biggest=0.05*max(wave);
k=wavelength;
while(abs(wave(k))<=biggest)
k=k-1;
end
wave=wave(1:k);
```

#### Listing Two

```function [pos, input, gain]=ddjagc(input)%
% ddjagc takes the analog packet input and looks at the begin of this packet
% for the tones. It then finds the fourth bit from each zero crossing and
% stores the received voltage. When preamble is detected by a zero crossing
% occuring before the fourth bit. The stored voltages are averaged and then
% used to compute the gain necessary to make the signal .5V

</p>
% setup values and initialize settings
samplebit=4;              % what bit to sample
samplesperbit=100/3;      % the number of data points in one bit
point=1;                  % starting point in data stream
stoppoint=length(input);  % ending point in data stream
counter=-1;               % counter to determine when the samplebit is found
gain=[];                  % used to store the voltages and the resultant gain

</p>
% These two variables are used to determine when sign changes
% have occurred. This is used to detect zero crossings.
olds=sign(input(point));
news=olds;
% Loop through the input
while(point<stoppoint)
% check both positive and negative edges to model rectification.
if(((olds<news)|(olds>news))&(counter<0))
olds=news;
counter=round(samplebit*samplesperbit);
elseif((olds>news)|(olds<news))
% we are in preamble stop looking at tones
gain=mean(gain);
input=gain*input;
pos=point;
return;
end
if(counter==0)
% maximum signal must be 0.5 volts the absolute value is needed to
% check both the positive and negative peaks of the tones.
gain=[gain 0.5/abs(input(point))];
end
% increment through the data stream.
point=point+1;
% give it a buffer zone for noise issues and unclean zero crossings.
if(abs(input(point))>0.1)
news=sign(input(point));
end
if(counter>=0)
counter=counter-1;
end
end
```

#### Listing Three

```function output=ddjhfb(input, pos)%
% ddjhfb implements the high frequency boost block of the equalizer model. It
% checks the peaks of the preamble signal and if the average of all measured
% peaks is less than .35V then gain is added. Otherwise no gain is added.

</p>
% setup values and initialize settings
baud_rate=30e6;                    % The data bit rate
samplesperbit=100/3;               % number of samples bit data bit
point=pos;                         % current position in the data input stream
stoppoint=length(input);           % last position in the data stream
olds=sign(input(pos));             % sign of last bit
news=olds;                         % sign of the current measurements
counter=round(0.5*samplesperbit);  % counter to determine when to sample peak
peaks=[];                          % vector that contains the peak voltages

</p>
% This is a 2 zero/3 pole filter
z=[3.0 -10]'*(-2e6*pi);            % zeros of the filter
p=[10 29 26]'*(-2e6*pi);           % poles of the filter

</p>
% the minus sign is because of the 180 degree phase shift of the filter
k=-(p(1)*p(2)*p(3))/(z(1)*z(2));   % DC gain of the filter

</p>
% Create filter function for equalizer. Here we use the bilinear transform to
% change from the analog domain in which the filter is designed to the
% digital domain that MatLab uses for implementing filters.
[zd, pd, kd]=bilinear(z, p, k, baud_rate*samplesperbit);
[b,a]=zp2tf(zd,pd,kd);

</p>
% loop through the input stream
while(point<stoppoint)
if(((olds<news)|(olds>news))&(counter<0))
olds=news;
counter=round(0.5*samplesperbit);
elseif((olds==news)&(counter<-0.5*samplesperbit))
% preamble is over. Filter the input if necessary.
peaks=mean(peaks);
if(peaks<0.35)
% filter the input
output=filter(b, a, input);
else
output=input;
end
return;
end
if(counter==0)
% maximum signal must be .5 volts. absolute value is used to check both
% the positive and negative peaks of preamble.
peaks=[peaks abs(input(point))];
end
% increment through the data stream.
point=point+1;
% give it a buffer zone for noise issues and unclean zero crossings.
if(abs(input(point))>.1)
news=sign(input(point));
end
counter=counter-1;
end

</p>

<a href="#rl3">Back to Article</A>
</P>
<H4><a name="l4">Listing Four</H4>

function [data, pos]=ddjpll(input, pos)
%
% ddjpll attempts to model the pll of the receiver block. This model
% is not really a pll but a model that is similar in function. This
% model measures the length of the preamble bits and determines what
% the clockperiod should be from those measurements. Preamble is a
% 101010... pattern so there is a zero crossing every bit time.

</p>
samplesperbit=100/3;        % the number of samples per bit
point=pos;                  % starting point in the input stream
datayet=0;                  % flag that is set to 1 when data is reached
olds=sign(input(point));    % the polarity of the old data bit
news=olds;                  % the polarity of the current data bit
nextcheck=0;                % counter that determines when next crossing occurs
clockperiod=0;              % average of all the times between zero crossings
period=[];                  % vector containing times between zero crossings

</p>
% loop through input stream until data is found
while(~datayet)
% don't get the new sign until it is expected
if(nextcheck<0)
news=sign(input(point));
end
% if the old sign and new sign do not equal then a zero
% crossing occured so measure the time for the bit.
if(news~=olds)
period=[period (samplesperbit-(nextcheck+1))];
nextcheck=round(samplesperbit);
olds=news;
end
% if a zero crossing has not occurrred within two bit times we are in data.
if(nextcheck<-samplesperbit)
datayet=1;
pos=point;
end
nextcheck=nextcheck-1;
point=point+1;
end
% determine the clock period
clockperiod=mean(period(2:end));
% look for the transition in the startdata word
% then use this transition to start data recovery
databegin=0;
olds=sign(input(point));
while(~databegin)
if(olds~=sign(input(point)))
databegin=1;
end
point=point+1;
end
% shift to the peaks for data sampling
point=point+clockperiod/2;
% now recover data using a simple zero crossing slicer model.
data=[];
while(point<length(input)-clockperiod)
if(input(round(point))>=0)
data=[data 1];
else
data=[data 0];
end
point=point+clockperiod;
end
```

#### Listing Five

```function errs=ddjber(dpacket, decodedata)%
% ddjber is a function that compares the original data packet
% with the decoded data at the output of the receiver block.
errs=0;
% strip the first 4 bits off decodedata because they
% are part of the startdata word
decodedata=decodedata(5:end);
point=1;
while(point<length(dpacket))
if(dpacket(point)~=decodedata(point))
errs=errs+1;
end
point=point+1;
end
```

#### Listing Six

```function ddjsystem%
% This is the main function that models our communication system.
% Every block of the system is called from this program.
% clear any figures open
close all;

% definitions
samplesperbit=100/3; % the number of points in a transmitted bit
% first create the sample packet to be sent
tone=[1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0];
tone=[tone tone tone tone tone];
preamble=[1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0];
startdata=[1 1 1 1 0 0 0 0];
data=[1 0 1 0 1 1 0 1 1 0 0 1 0 1 0 0 0 1 1 0 1 0 1 1 0 1 0 0 1 0 0 1];
data=[data data data];
enddata=[0 0 0 0 1 1 1 1];
% the digital packet
dpacket=[tone preamble startdata data enddata];
% create the analog packet by using the shift-flip-add routine
% load the bit response data
plot(apacket)
zoom on;
% Now lock the AGC value on the tone portion of the waveform
% determine the gain needed to maintain voltage regulation,
% multiply the apacket by gain factor and then return the new apacket.
[pos, apacket, gain]=ddjagc(apacket);
hold on;
plot(apacket,'r');
disp(['preamble start position ' num2str(pos) ' agc gain ' num2str(gain)]);
% now lock to the peaks of the preamble to set the HFB gain.
% Send input through the filter and return the apacket.
apacket=ddjhfb(apacket, pos);
plot(apacket,'k');
grid on;
% now lock clock recovery circuit to preamble and decode digital data.
[decodedata, pos]=ddjpll(apacket, pos);
disp(['data start position ' num2str(pos)]);
% finally compare the decoded data with dpacket.
errors=ddjber(data, decodedata);
```

DDJ

### More Insights

 To upload an avatar photo, first complete your Disqus profile. | View the list of supported HTML tags you can use to style comments. | Please read our commenting policy.