%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Function: PLL.m
%Programmed by: Qian Wang
%Date: September 20, 2006
%Version: 1.0
%Description: This simulation program is to realize Phase-Looked Loops.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [vco_out,pd_out,lf_out,phi_vco] = PLL(sig_in,fs,num,den,K_vco,K_pd,tyofpd,f_vco,M,N)
%-----------------------------------------------------------
%                      Input parameters
% sig_in:   Input siganl (vector)
% fs:       Sampling rate
% num:      Numerators of Loop filter transfer function in Laplace form (vector)
% den:      Denominators of Loop filter transfer function in Laplace form (vector)
% K_vco:    VCO gain (rad/(V*s))
% K_pd:     Phase detector gain (V/rad)
% tyofpd:   Type of phase detector (1 = Multiplier; 2 = XOR)
% f_vco:    VCO center frequency
% M:        Frequency divider
% N:        Frequency divider
%                      Output parameters
% voc_out:  Output signal of VCO
% pd_out:   Output signal of phase detector
% lf_out:   Output signal of loop filter
% phi_vco:  VCO phase vector
%-----------------------------------------------------------

%********************** Preparation part ***************************

npts = length(sig_in);                 %Number of simulation points
T    = 1/fs;                           %Sample time
phi_vcoinitial = 2*pi*rand;          %Initialize the total phase before the PLL algorithm is started

%Initialize the VCO output before the PLL algorithm is started
vco_outinitial = sin(phi_vcoinitial);
%  vco_outinitial = square(phi_vcoinitial);

%Preallocating arrays
phi_vco = zeros(1,npts+1);             %Initialize VCO phase
vco_out = zeros(1,npts+1);             %Initialize VCO output
pd_out  = zeros(1,npts+1);             %Initialize phase detector output
lf_out  = zeros(1,npts+1);             %Initialize loop filter output

[numd,dend] = bilinear(num,den,fs);    %Transforms analog filters into their discrete equivalents
[Numd,Dend] = eqtflength(numd,dend);   %Equalize lengths of transfer function's numerator and denominator
[a,b,c,d]   = tf2ss(Numd,Dend);        %Convert transfer function loop filter parameters to state-space form

mn = size(a);
x  = zeros(mn(2),npts+1);              %Initialize state variable

%************************** Main loop part **************************

h = waitbar(0,'Please wait...');

for i = 1:npts
    waitbar(i/npts,h);

    %***************** Phase detector ********************

    %Initialize the start point at t = 0
    if i == 1
        phi_vco(i) = phi_vcoinitial;
        vco_out(i) = vco_outinitial;

        %Choose the type of phase detector
        if tyofpd == 1
            pd_out(i) = K_pd*sig_in(i)*vco_out(i);
        elseif tyofpd == 2
            pd_out(i) = -sign(sig_in(i)*vco_out(i));
        end
    else
        %Choose the type of phase detector
        if tyofpd == 1
            pd_out(i) = K_pd*sig_in(i)*vco_out(i);
        elseif tyofpd == 2
            pd_out(i) = -sign(sig_in(i)*vco_out(i));
        end
    end

    %******************** Loop filter *********************

    %Using state-space method to calculate loop filter output
    x(:,i+1)  = a*x(:,i) + b*pd_out(i);
    lf_out(i) = c*x(:,i) + d*pd_out(i);

    %******************** VCO *********************
    %Calculate VCO phase
    phi_vco(i+1) = phi_vco(i) + (2*pi*f_vco + K_vco*lf_out(i))*T;

    %Calculate VCO output
    vco_out(i+1) = sin(phi_vco(i+1));
    %         vco_out(i+1) = square(phi_vco(i+1));
  
end

close(h);

%%%%%%%%%%%%%%%%%%%%%% end of file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%