Modeling with MatLab
by Mark Weaver



Example 1:

[b,a]=zp2tf([z1 z2],[p1 p2 p3];
output=filter(b,a,input);


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.

% 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

% 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

% 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

% 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
   % begin AGC adapt
   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.

% 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

% 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

% 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

% 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);

% 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

Listing Four
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.

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

% 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 data/bit_response;    
% load the bit response data
apacket=ddjshiftflipadd(dpacket, bit_response2, samplesperbit);
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);
disp(['received errors ' num2str(errors)]);



1