%Simulation of ADPLL (Fig. 10.18)
clear

%********** ENTER PARAMETERS BELOW
%Initial State and Simulation Time
	Fin = 400000; % Inputfrequency (Hz)
	Tmaximum = .6E-4; % Simulation duration (seconds)
	Tminimum = 0; % Start of Plots (seconds)
	Noff = 0; % add offset to produce N2
	N2 = 200; % ³ Noff, initial input to OA (indicates frequency)
	Nout(1) = 895; % Initial output count (indicates phase)
	Phin = 0.5; % clock cycles when first input transition occurs
%Configuraton
	q = 4; % OA shifted 2^-q to N2
	nv = 12; % OA capacity, Nv = 2^nv
	Fclock = 8.192e6; % OA clock
%**********

%Derived
	Tclock = 0; % start time
	Twottq = 2^(-q);
	Nv = 2^12;
	Fstep = Fclock/Nv; % minimum step size
	N2ss = Fin/Fstep; % N2 for steady state
	NatS = Nv - 1 - (N2ss- Noff)/Twottq; % Nout at steady state
	Fmax = Fstep*(Noff + floor((Nv-1)*2^(-q))); % frequency with maximum N2 (max freq)
	Fmax = Fmax*(1-Fmax/(2*Fstep*Nv)); % to allow excursion (N2/2) about average NatS
	Fmin = Fclock*2^(-q-1); % minimum for stability
	Tin = Phin/Fclock; % time of first input transition 
	Nmax = ceil(Tmaximum*Fclock); % number of samples in simulation
	if Nmax > 512 % check for fewer than half maximum array elements
		fprintf('\n%g samples. Must be fewer than 513 in MATLAB 3.5.  Reduce Tmaximum.\n',Nmax)
	else % continue simulation
	 Tmax = Nmax/Fclock;
	 Nmin = floor(Tminimum*Fclock); % sample number at start of plot
	 Tmin = Nmin/Fclock;
%Output
	 fprintf('\nFclock = %g Hz; Fin = %g Hz (period %g sec.)\n',Fclock,Fin,1/Fin);
	 fprintf('Initial:      N2 = %g; Nout = %g\n',N2,Nout(1));
	 fprintf('Steady State: N2 = %g; Nout = %g at Input Transition\n',N2ss,NatS);
	 fprintf('Nv = 2^%g = %g; Noff = %g\n',nv,Nv,Noff);
	 fprintf('q = %g so 2^(-q) = %g = 1/%g\n',q,Twottq,1/Twottq);
	 fprintf('Minimum step is %g Hz\n',Fstep);
	 fprintf('Allowed range: %g Hz < Fin < appx. %g Hz\n',Fmin,Fmax);
	
%Simulation

 j=1; k=1; Iover(1) = 1;
 for i=1:Nmax,
	 Tclock = Tclock + 1/Fclock; %increment clock time; time at start of Nout(i+1)
 	if Tclock > Tin % an input transition has occurred (before this clock)
	 	N2 = Noff + floor((Nv - 1 - Nout(i))*2^(-q)); % Nout (before clock) was loaded in register,
			 	% becoming N2 after complementing, decrementing and offsetting
	 	Tin = Tin + 1/Fin; % next input transition
		 Isamp(j) = i; % index of output giving new N2 value 
	 	N2i(j) = N2; % new N2
		 j = j+1;
 	end; % if Tclock
	 Nnext= Nout(i) + N2; % tentative Nout after clock
	 if Nnext >= Nv, % overflow of OA occurs
		 Nout(i+1) = rem(Nnext,Nv); % effect of overflow
		 Iover(k+1) = i; % index of Nout just before overflow
	 	Finst(k+1) = Fclock/(Iover(k+1) - Iover(k)); % Instantaneous frequency
		 		% (averaged over one output period)
	 	k = k+1;
	 else,
	 	Nout(i+1) = Nnext;
 	end % If Nnext
 end; %for

%Plotting

 clg;	
 t = linspace(Tmin,Tmax+1/Fclock,Nmax+2-Nmin); % times at clock instants
 for m = 1:i-Nmin, % change sampled points into steps
	 T(2*m-1) = t(m);
	 T(2*m) = t(m+1);
 	No(2*m-1) = Nout(m+Nmin);
	 No(2*m) = Nout(m+Nmin);
 end; % m
 axis([Tmin Tmax 0 Nv])
 subplot(211); plot(T,No) % stepped Nout plot
 title('Nout with sampled values shown')
 xlabel('seconds')
 for m=1:j-1 % form matrix of sampled values and instantaneous frequencies
	 zz(m) = Nout(Isamp(m));
 end % m
 hold on; plot((Isamp-.5)/Fclock, zz,'g*'); hold off % plot values of Nout  sampled for N2

 for m=1:j-2 % form matrix of instantaneous frequencies
 	Is(2*m-1) = Isamp(m);
 	Is(2*m) = Isamp(m+1);
	 Fi(2*m-1) = N2i(m)*Fclock/Nv;
 	Fi(2*m) = Fi(2*m-1);
 end % m

 fmean = mean(Fi);
 subplot(212); plot([Is(1) Is(2*j-4)]/Fclock, [fmean fmean],'b');
 hold on; plot(Is/Fclock,Fi,'r'); hold off % instantaneous frequency vs. time
 vv = axis; % get auto axis ranges
 axis([Tmin Tmax vv(3) vv(4)]) % use same x values as in first plot
 title('Instanteous & Mean Frequencies')
 xlabel('seconds @ mid-period')
 ylabel('Hz')

 fprintf('Mean frequency over %g seconds: %g Hz (Fin = %g Hz)\n\n',(Iover(k-1)-Iover(2))/Fclock,fmean,Fin)
		 % average frequency over the simulation time
	end % simulation
