%Synthesizer Simulation
%Time Responses and Phase Plane
%Sample & Hold PD, 0-1 cycle range
%  Matrix starts at 0 initial conditions

clear;

% MODIFY PARAMETERS BELOW & GO
%*****************************
Fref = 10000;  % Reference frequency in Hz
N1 = 1000;  % divide number during simulation
N0 = 750;  % prior divide number
PNL = .5;  % radian final phase, PHA
Kp = 3;  % PD gain, V/cycle (modify other units appropriately if A/cycle used)
KLF = 60000;  % Filter DC gain
Fp1 = .01;  % Filter first pole in Hz
Fz = 300;  % Filter zero in Hz
Fp2 = 3000;  % Filter second pole in Hz
Kv = 1e6;  % VCO gain, Hz/V
Tend = 4.8e-3; % Length of simulation.  Student Matlab Fref*Tend ² 1k (v3.5) or 8k (v4) or 16k (v5)
truncatePh = 1; % Phase Plane Plot uses 1 = Phe, 0 = Pha
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 = 3; % Choose Output
ACC = 1e-8; % phase error that ends iteration
%******************************
FNL = N1*Fref; F0 = N0*Fref; K = Kp*KLF*Kv; % final frequency, initial frequency, K
Pha0 = PNL+(F0 - FNL)/K; % initial phase.  Æ¿ * K = ÆF, Æ from initial to final
Tref = 1/Fref;  % reference period
fprintf('\nLoop with 0 - 1 cycle Sample-and-Hold Phase Detector');
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('\nPhase Error from %g cycle to %g cycle\n', Pha0,PNL);
twopi = 2*pi;
wp1 = twopi * Fp1; % first nonzero pole in rad/sec
wp2 = twopi * Fp2; % second nonzero pole in rad/sec
wz = twopi * Fz; % zero in rad/sec
Kx = wp1*wp2*K/wz; % some constants to be used
wps = wp1 + wp2;
wpm = wp1*wp2;

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

%Polynominals for open loop gain:
	N = Kx*[1    wz]; % numerator
	D = [1    wps   wpm    0 ]; % denominator for phase
	Df = [1    wps    wpm]; % denominator for frequency

[a,b,c,d] = tf2ss(N,D); % elements of dynamic matrix equation for phase
x0 = [0; 0; 0]; % zero initial conditions for phase and derivatives
yywas = 0;
[af,bf,cf,df] = tf2ss(N,Df); % elements of dynamic matrix equation for freq.
xf0 = [0; 0]; % zero initial conditions for frequency and derivative

%    BEGIN SIMULATION

TT(1) = 0;
Pha(1) = Pha0; % phase
Phe(1) = Pha0; % PD output
FF(1) = F0; % frequency
i = 1;
T = Tref;
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,
		[Phi, Gamma] = C2D(a,b,T); % matrices for use in solution eqution
		x = ltitr(Phi,Gamma,u1',x0); % state variable solution for 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); % for use in computing next point
	TT(i+1) = T + TT(i); % Time of sample
	Ph = Pha(i) + T*Fref - 1; % Phase at sample
	Pha(i+1) = Ph;
	Phe(i+1) = Ph - floor(Ph); % limit range of this variable to 0-1 cycle
	[Phif, Gammaf] = C2D(af,bf,T); % for use in frequency solution
	xf = ltitr(Phif,Gammaf,u1',xf0); % for use in frequency solution
	yyf = cf*xf'; % frequency less offset
	FF(i+1) = yyf(2) + F0; % frequency out
	i = i+1;
	x0 = x(2,:)'; % set initial x for next iteration as current x
	xf0 = xf(2,:)'; % set initial x for next iteration as current x
end; %while TT

% END SIMULATION.  The rest is plotting.

%  Set multipliers for chosen units
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

%  Create each of the 4 plot types and locate them depending on what combination was chosen.

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
	if truncatePh,
		plot(Pscale*Phe,Fscale*FF,'r')
		title('Phase Plane, Freq. vs. Phe')
	else
		plot(Pscale*Pha,Fscale*FF,'r')
		title('Phase Plane, Freq. vs. Pha')
	end; %if truncate
	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')
	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')
	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 
	plot(TT,Phe)
	title('Phe, PD Output/(Kp cycle)')
	grid
	xlabel('seconds')
end % if ()|()


