%Synthesizer Simulation
%Time Responses and Phase Plane
%Sample & Hold PD, 0-1 cycle range
%Finite efficiency, ramp flat at recycle and sample
%Approximations re recycle shape, recycle during sample, min. sample spacing
clear;

% MODIFY PARAMETERS BELOW & GO
%*****************************
Fref = 10000;  % Reference frequency in Hz
N1 = 1000;  % divide number during simulation
N0 = 1120;  % initial divide number
PNL = .47;  % cycle final PD output, Phe
Kp = 4;     % PD gain in V/cycle
KLF = 5e4;  % Filter DC gain
Fp1 = 1e-3; % Filter first pole in Hz
Fz = 100;   % Filter zero in Hz
Fp2 = 1000; % Filter second pole in Hz
Kv = 1e6;   % VCO gain in Hz/V
RCYL = 0.1; % ramp recycle time in cycles
FLAT = 0.05;% flat part of ramp during sample in cycles
Eff = .8;   % PD Efficiency
Tend = 4.8e-3;  % Length of simulation.  Student Matlab Fref*Tend ² 1k (v3.5) or 8k (v4) or 16k (v5)
truncatePh = 1; % Pha displays: 0 = multicycle  1 =  mod 1 cycle
Punits = 1;     % Plot Phase Units: 0 = radians  1 = cycles  2 = degrees
Funits = 1;     % Plot Frequency Units: 0 = radians/sec   1 = Hz
% Available Outputs
%  op=1:  Phase Plane = Freq. vs. Pha
%  op=2:  Phase Plane + Phase Error vs. Time
%  op=3:  Frequency + PD Output vs. Time
%  op=4:  All of above = 4 plots
op = 4; % Choose Output
ACC = 1e-8; % phase error that ends iteration
%******************************
FNL = N1*Fref; % final frequency
F0 = N0*Fref;  % initial frequency
K = Kp*KLF*Kv; 
Pha0 = PNL+(F0 - FNL)/K; % initial phase
Tref = 1/Fref;
fprintf('\nLoop with 0-1 cycle Sample-and-Hold Phase Detector');
fprintf('\nRamp flat for %g cycle after recycle and %g cycle during sample',RCYL,FLAT)
fprintf('\nPD Efficiency = %g',Eff);
fprintf('\nKp = %g V/cycle; KLF = %g; Kv = %g Hz/V',Kp,KLF,Kv);
fprintf('\nFp1 = %g Hz; Fz = %g Hz; Fp2 = %g Hz',Fp1,Fz,Fp2);
fprintf('\nFref = %g Hz; N goes from %g to %g',Fref, N0, N1);
fprintf('\nFrequency from %g Hz to %g Hz;', F0,FNL);
fprintf('\nPD Out from %g cycle to %g cycle\n', Pha0,PNL);

twopi = 2*pi;
wp1 = twopi * Fp1;
wp2 = twopi * Fp2; % add additional poles below this line & add names to p below
wz = twopi * Fz;
p = [-wp1,-wp2]; % form a vector of open-loop poles
z = [-wz]; % form a vector of open-loop zeros
Kx = wp1*wp2*K/wz;

%Polynominals for open loop gain:
 [N,Df] = zp2tf(z,p,Kx); % convert zeros, poles, and gain factor to open-loop numerator and denominator.
 sd = size(Df); % array size; sd(2) = number of columns = number of poles
 D = Df; D(sd(2)+1) = 0; % multiply denominator by s by shifiting left (adding 0 to right)
 [a,b,c,d] = tf2ss(N,D); % elements of dynamic matrix equation
 [af,bf,cf,df] = tf2ss(N,Df);

clf; % clear previous graphs; for v.3.5, use clg #########

%    INITIAL CONDITIONS

x0 = [0; 0; 0];
xf0 = [0; 0];
yywas = 0;
TT(1) = 0; % initial sample time
Pha(1) = Pha0; % initial phase
Phe(1) = Pha0; % initial PD output
FF(1) = F0; % initial frequency
i = 1;
flrPhOld = 0; smplAgain = 0;
T = Tref; % initial approximation of sample period

%    BEGIN SIMULATION

while TT(i)<Tend,
	u1(1) = Phe(i) - Pha0; % input to sampler at start of period = held value
	u1(2) = u1(1); %  constant Hold output -> same at beginning and end of period
	y = 0; % phase change at divider output (initially 0 to insure loop starts)
	while abs(y-1)>ACC, % iteration to find T corresponding to one cycle at divider output
		[Phi, Gamma] = C2D(a,b,T); % matrices for use in phase solution eqution at T
		x = ltitr(Phi,Gamma,u1',x0); % state variable solution (for phase) at each value of t (only 2 in this case)
		yy = c*x'; % corresponding outputs (d is null for polynominals above)
					% x is transpose of X in PLB Eq. (6.33); yy is lower element in PLB Eq. (6.45)
					% x and yy have a row for each time increment
		y = (yy(2)-yywas + T*F0)/N1; % ÖN output phase change since last cycle
		T = T/y; % adjusting toward y = 1
	end; % while abs(y-1)
	yywas = yy(2);
	TT(i+1) = T + TT(i);
	Ph = Pha(i) + T*Fref - 1;
	flrPh = floor(Ph); % closest integer below Ph (0 in steady state lock if Pha truncated)
	if flrPh < flrPhOld, % another sample occurred on same ramp
		smplAgain = smplAgain + 1; % add one to number of extra samples on this ramp
	else
		smplAgain = 0; % reset count - this sample was on next ramp.
	end
	flrPhOld = flrPh;	
	Phe(i+1) = Ph - flrPh; % truncate Phe (Phe is mod one cycle)
	if truncatePh,
		Pha(i+1) = Phe(i+1); % Pha truncated too
	else
		Pha(i+1) = Ph; % Pha not truncated
	end; % truncate Pha (show phase mod 1 cycle)
	if Phe(i+1) < RCYL, % sampel during finite recycle time
		Phe(i+1) = 0; % exponential dependence on time into RCYL can be introduced here for better accuracy
	else
		Phe(i+1) = Phe(i+1) - RCYL; % ramp does not start until RCYL into period
	end; 
	Phe(i+1) = Phe(i+1) + FLAT*smplAgain; % subtract time ramp was stopped during earlier samples
	% *********RESTRICTION: divider output period must not be less than FLAT
	Phe(i+1) = Phe(i)*(1-Eff) + Phe(i+1)*Eff; % efficiency [Eq. (5.16)]
	[Phif, Gammaf] = C2D(af,bf,T); % matrices for use in frequency solution equation
	xf = ltitr(Phif,Gammaf,u1',xf0); % state variable solution (for frequency) at each value of T
	yyf = cf*xf'; % similar to yy above
	FF(i+1) = yyf(2) + F0; % add frequency offset for total frequency at sample time
	i = i+1;
	x0 = x(2,:)'; %initial state variables for phase equation next period same as final values this period
	xf0 = xf(2,:)'; %initial state variables for frequency equation next period same as final values this period
end; %while TT

% END SIMULATION.  The rest is plotting.

if Punits == 1,
	Pscale = 1;
elseif Punits == 2,
	Pscale = 360;
else,
	Pscale = 2*pi;
end % if Punits

if Funits == 1,
	Fscale = 1;
else
	Fscale = 2*pi;
end % if Funits

if op ~= 3, % selecting phase-plane plot
	if op==1,
		subplot(111),
	elseif op==2,
		subplot(211),
	elseif op==4,
		subplot(221),
	end % if op==1
		plot(Pscale*Pha,Fscale*FF,'r')
		title('Phase Plane, Freq. vs. Pha at samples')
	grid
	if Punits == 1, % Set display of phase units
		xlabel('phase error (cycles)');
	elseif Punits == 2,
		xlabel('phase error (degrees)');
	else,
		xlabel('phase error (radians)')
	end % if Punits
	if Funits == 1, % Set display of frequency units
		ylabel('freq (Hz)');
	else,
	 ylabel('freq (rad/sec)');
	end % if Funits
end % if op~=3

if (op == 2)|(op == 4), % selecting phase error vs. time
	if op == 2,
		subplot (212),
	else,
		subplot (223),
	end % if op==2
	plot(TT,Pscale*Pha,'g')
	title('Pha, Phase Error at samples')
	grid
	xlabel('seconds')
	if Punits == 1, % Set display of phase units
		ylabel('phase (cycles)');
	elseif Punits == 2,
		ylabel('degrees');
	else,
		ylabel('radians');
	end % if Punits
end % if ()|()

if (op == 3)|(op == 4), 
	if op == 3, % selecting frequency vs. time
		subplot(211);
	else,
		subplot(222),
	end % if op==3 
	plot(TT,Fscale*FF,'b')
	title('Frequency Out at samples')
	grid
	xlabel('seconds');
	if Funits == 1, % Set display of frequency units
		ylabel('Hz');
	else,
	 ylabel('radians / second');
	end % if Funits

	if op == 3, % selecting PD output vs. time
		subplot(212);
	else,
		subplot(224),
	end % if op==3
	mx = length(Phe);
	TT(mx+1) = TT(mx)  + Tref;
	for i = 1:mx % form a new vector giving values of Phe before and after step change
		T(2*i-1) = TT(i);
		T(2*i) = TT(i+1);
		Ph(2*i-1) = Phe(i);
		Ph(2*i) = Phe(i);
	end; %i 
	plot(T,Ph,'m')
	title('Phe, PD Output/(Kp cycle)')
	grid
	xlabel('seconds')
end % if ()|()


