#!/<yourPath>/perl -w

use POSIX;

# rev 1.1           added estimate of phase margin. Modified from Mansuri, JSSC, Nov. 2002, p. 1377
# rev 1.2, 7/21/04, added Gardner stability criterion and thermal noise for resistor and capacitors
# rev 1.3  7/23/04  corrected reporting of magnitude in dB. forgot to divide log(mag) by log (10)
$rev = "1.3";

printf(STDOUT "\n# plloop.pl rev <%s>. Copyright 2004, Dennis Fischette\n# Generates PLL transfer functions and Loop Equations\n",$rev); 
#printf(STDOUT "\n# plloop.pl rev <%s>.\n# Generates PLL transfer functions and Loop Equations\n",$rev); 

if (($#ARGV != 10) && ($#ARGV != 12)) {
    printf STDERR "$0 <transferFunctionOutputPrefix> <kvco> <icp> <res> <c1> <c2> <fbdiv> <fbdly> <tref> <ig> <rmsPctPerJit> [<tau3> <tau4>]\n";
    printf STDERR "\n\tkvco:  VCO gain (Hz/V). e.g. 1e9\n";
    printf STDERR "\ticp:   charge-pump current (Amp). e.g. 30e-6\n";
    printf STDERR "\tres:   low-pass filter resistance (Ohm). e.g. 3000\n";
    printf STDERR "\tc1:    low-pass filter large capacitor (Farad). e.g. 100e-12\n";
    printf STDERR "\tc2:    low-pass filter small (smoothing) capacitor (Farad). e.g. 7.5e-12\n";
    printf STDERR "\tfbdiv: VCO feedback divider. e.g. 10\n";
    printf STDERR "\tfbdly: feedback path delay. e.g. 5e-9\n";
    printf STDERR "\ttref:  reference period (Sec). e.g. 20e-9\n";
    printf STDERR "\tig:    total gate-leakage on low-pass filter capacitors (Amp)\n";
    printf STDERR "\trmsPctPerJit: estimated std deviation of VCO period jitter as percentage of VCO period. e.g. 0.4\n";
    printf STDERR "\ttau3:  (optional) time-constant (Sec) of parasitic pole between LPF and VCO V-to-I e.g. 3e-9. Default=0\n";
    printf STDERR "\t\tmodels additional LPF ripple filter\n";
    printf STDERR "\ttau4:  (optional) time-constant *(Sec) of parasitic pole between parasitic pole3 and VCO e.g. 1e-9. Default=0\n";
    printf STDERR "\t\tmodels VCO V-to-I filter\n";
    exit(1);
}
$PI = 3.141592654;
$Kboltzmann = 1.38e-23;
$temperature = 373.0;

$tau3 = 0;
$tau4 = 0;

$outputfilePrefix = $ARGV[0];
$kv = $ARGV[1];
$icp = $ARGV[2];
$res = $ARGV[3];
$c1 = $ARGV[4];
$c2 = $ARGV[5];
$div = $ARGV[6];
$dly = $ARGV[7];
$tpfd = $ARGV[8];
$igate_leakage = $ARGV[9];
$jitterPctRms = $ARGV[10]/100.0;
if ($#ARGV == 12) {
    $tau3 = $ARGV[11];
    $tau4 = $ARGV[12];
}

$outputfile = sprintf("%s.%s",$outputfilePrefix,"out");
$outputfile_closed_mag = sprintf("%s.%s",$outputfilePrefix,"closedMag");
$outputfile_closed_phs = sprintf("%s.%s",$outputfilePrefix,"closedPhs");
$outputfile_open_mag = sprintf("%s.%s",$outputfilePrefix,"openMag");
$outputfile_open_phs = sprintf("%s.%s",$outputfilePrefix,"openPhs");

if (!open(OUTPUT,"> $outputfile")) {
    printf STDERR "error: cannot open file <%s> for write\n",$outputfile;
    exit(1);
}
if (!open(OUTPUT_CLOSED_MAG,"> $outputfile_closed_mag")) {
    printf STDERR "error: cannot open file <%s> for write\n",$outputfile_closed_mag;
    exit(1);
}
if (!open(OUTPUT_CLOSED_PHS,"> $outputfile_closed_phs")) {
    printf STDERR "error: cannot open file <%s> for write\n",$outputfile_closed_phs;
    exit(1);
}
if (!open(OUTPUT_OPEN_MAG,"> $outputfile_open_mag")) {
    printf STDERR "error: cannot open file <%s> for write\n",$outputfile_open_mag;
    exit(1);
}
if (!open(OUTPUT_OPEN_PHS,"> $outputfile_open_phs")) {
    printf STDERR "error: cannot open file <%s> for write\n",$outputfile_open_phs;
    exit(1);
}

$period = $tpfd/$div;
$vcoFreq = 1/$period;
$jitterRms = $period * $jitterPctRms;

$rc1 = $res * $c1;
$rc2 = $res * $c2;
$wn = sqrt($kv*$icp/($div*$c1));
$wn_in_hz = 0.5 * $wn/ $PI;
$damping = 0.5 * $res * $c1 * $wn;
$loopgain = $kv * $icp / $div;
($wn_upper_limit_in_hz,$res_upper_limit) = &calc_discrete_time_stability_limit($res,$c1,$wn,$tpfd);    

# transfer function gain, phase margin
@one = (1,0);
@neg_one = (-1,0);
#@two = (2,0);
#@wz_c = (1/$rc1,0);
#@wn_c = ($wn,0);
#@damping_c = ($damping,0);
#@wn_sqr = &mul_c(\@wn_c,\@wn_c);

# 2nd-order closedloop
# @h2 = @wn_c*@wn_c * (@one + @s*@wz_c)/ (@s*@s + @two*@s*@damping_c*@wn_c + @wn_c*@wn_c);

$rc1_c[0] = $rc1;
$rc1_c[1] = 0;
$ct_c[0] = $c1+$c2;
$ct_c[1] = 0;
$c2_c[0] = $c2;
$c2_c[1] = 0;
$k_c[0] = $loopgain;
$k_c[1] = 0;
$tau3_c[0] = $tau3;
$tau3_c[1] = 0;
$tau4_c[0] = $tau4;
$tau4_c[1] = 0;
$dly_c[0] = $dly;
$dly_c[1] = 0;

printf(OUTPUT "# plloop.pl rev <%s>. Copyright 2004, Dennis Fischette\n# PLL transfer functions\n",$rev); 
printf(OUTPUT_CLOSED_MAG "# plloop.pl rev <%s>. Copyright 2004, Dennis Fischette\n# PLL transfer functions\n",$rev); 
printf(OUTPUT_CLOSED_PHS "# plloop.pl rev <%s>. Copyright 2004, Dennis Fischette\n# PLL transfer functions\n",$rev); 
printf(OUTPUT_OPEN_MAG "# plloop.pl rev <%s>. Copyright 2004, Dennis Fischette\n# PLL transfer functions\n",$rev); 
printf(OUTPUT_OPEN_PHS "# plloop.pl rev <%s>. Copyright 2004, Dennis Fischette\n# PLL transfer functions\n",$rev); 
printf(OUTPUT "%s\n", "# note: transfer functions assume continuous-time behavior (f(ref) >> PLL bandwidth)"); 
printf(OUTPUT "%s\n", "# note: transfer functions are 3rd-order or 5th-order"); 
printf(OUTPUT "#\t\t\t\t%s\t\t\t%s\n", "closed-loop", "open-loop"); 
printf(OUTPUT "%s\t%s\t\t%s\t%s\t%s\t%s\n", "#log(f)" ,"f(Hz)", 
                     "Gain(dB)", "Phase (degrees)", 
                     "Gain(dB)", "Phase (degrees)"); 

printf(OUTPUT_CLOSED_MAG "%s\n", "# Closed-loop Magnitude transfer function"); 
printf(OUTPUT_CLOSED_MAG "%s\t%s\n", "#log(f)","Gain(dB)"); 
printf(OUTPUT_CLOSED_PHS "%s\n", "# Closed-loop Phase transfer function"); 
printf(OUTPUT_CLOSED_PHS "%s\t%s\n", "#log(f)","Phase(degrees)"); 
printf(OUTPUT_OPEN_MAG "%s\n", "# Open-loop Magnitude transfer function"); 
printf(OUTPUT_OPEN_MAG "%s\t%s\n", "#log(f)","Gain(dB)"); 
printf(OUTPUT_OPEN_PHS "%s\n", "# Open-loop Phase transfer function"); 
printf(OUTPUT_OPEN_PHS "%s\t%s\n", "#log(f)","Phase(degrees)"); 

$f = 1;
$prev_freq2 = $f;
$points_per_dec = 100;
$max_magnitude_closedloop_dB = -1e6;
$magnitude_closedloop_dB = -1e6;
$magnitude_openloop_dB = -1e6;
$phase_openloop = 0.0;
for ($decade=0; $decade < 10; $decade++) {
  for ($points=1; $points < 9; $points = $points + (1/$points_per_dec)) { 
    $prev_freq2 = $freq2;
    $freq2 = $points*$f; 
    $omega = 2*$PI*$freq2; 
    @s     = (0.0, $omega);   

    # num3g = ((K/s)*(1+s*rc1)) * (e**(-s*Tdly)) 
    @num3g = &mul_c(\@s,\@rc1_c);
    @num3g = &add_c(\@one,\@num3g);
    @num3g = &mul_c(\@k_c,\@num3g);
    @num3g = &div_c(\@num3g,\@s);
    @temp2 = &mul_c(\@s,\@dly_c);
    @temp2 = &mul_c(\@neg_one,\@temp2);
    @temp2 = &exp_c(\@temp2);
    @num3g = &mul_c(\@temp2,\@num3g);
    # denom3g = (s*(ct+(s*rc1*c2))) * (1+s*Tau3) * (1+s*Tau4) 
    @denomTau3 = &mul_c(\@s,\@tau3_c);
    @denomTau3 = &add_c(\@one,\@denomTau3);
    @denomTau4 = &mul_c(\@s,\@tau4_c);
    @denomTau4 = &add_c(\@one,\@denomTau4);
    @temp2 = &mul_c(\@rc1_c,\@c2_c);
    @temp2 = &mul_c(\@s,\@temp2);
    @temp2 = &add_c(\@ct_c,\@temp2);
    @denom3g = &mul_c(\@s,\@temp2);
    @denom3g = &mul_c(\@denom3g,\@denomTau3);
    @denom3g = &mul_c(\@denom3g,\@denomTau4);

    # open-loop
    @g3 = &div_c(\@num3g,\@denom3g);

    # closed-loop
    @denom3 = &add_c(\@one,\@g3);
    @h3 = &div_c(\@g3,\@denom3);

    $prev_magnitude_openloop_dB = $magnitude_openloop_dB;
    $prev_magnitude_closedloop_dB = $magnitude_closedloop_dB;
    $prev_phase_openloop = $phase_openloop;
    $magnitude_closedloop_dB = 20.0 * (log(&mag_c(\@h3)) / log(10));
    $magnitude_openloop_dB   = 20.0 * (log(&mag_c(\@g3)) / log(10));
    $phase_openloop = &ang_c_degrees(\@g3);

    if ($magnitude_closedloop_dB > $max_magnitude_closedloop_dB) {
        $max_magnitude_closedloop_dB = $magnitude_closedloop_dB;
        $peaking_frequency = $freq2;
    }
    if (($prev_magnitude_closedloop_dB > 0) && ($magnitude_closedloop_dB <= 0)) {
        ($bandwidth_closedloop_0dB,$intr_fraction) = &interpolate($prev_freq2,$prev_magnitude_closedloop_dB,$freq2,$magnitude_closedloop_dB,0);
    }
    if (($prev_magnitude_closedloop_dB > -3) && ($magnitude_closedloop_dB <= -3)) {
        ($bandwidth_closedloop_minus3dB,$intr_fraction) = &interpolate($prev_freq2,$prev_magnitude_closedloop_dB,$freq2,$magnitude_closedloop_dB,-3);
    }
    if (($prev_magnitude_openloop_dB > 0) && ($magnitude_openloop_dB <= 0)) {
        ($bandwidth_openloop_0dB,$intr_fraction) = &interpolate($prev_freq2,$prev_magnitude_openloop_dB,$freq2,$magnitude_openloop_dB,0);
        $phase_margin = 180.0 + $prev_phase_openloop + $intr_fraction*($phase_openloop - $prev_phase_openloop);
    }
    if ((($points_per_dec * $points) % ($points_per_dec/10)) == 0) { 
        printf(OUTPUT "%.3f\t%.2e\t%12.3f\t%12.3f\t%12.3f\t%12.3f\n", log($freq2)/log(10), $freq2, 
                         $magnitude_closedloop_dB, &ang_c_degrees(\@h3), $magnitude_openloop_dB, $phase_openloop);
        printf(OUTPUT_CLOSED_MAG "%.3f\t%12.3f\n", log($freq2)/log(10), $magnitude_closedloop_dB);
        printf(OUTPUT_CLOSED_PHS "%.3f\t%12.3f\n", log($freq2)/log(10), &ang_c_degrees(\@h3));
        printf(OUTPUT_OPEN_MAG "%.3f\t%12.3f\n", log($freq2)/log(10), $magnitude_openloop_dB);
        printf(OUTPUT_OPEN_PHS "%.3f\t%12.3f\n", log($freq2)/log(10),$phase_openloop);
    }
  }
  $f = $f*10.0;
}

$inverse_time_constant = 2.0 * $wn * $damping;
$inverse_time_constant_in_hz = 0.5 * $inverse_time_constant/ $PI;
$osr =  1.0/($tpfd*$bandwidth_openloop_0dB);

if ($damping >= 1.0) {
    $jitter_longterm = $jitterRms * sqrt($vcoFreq/$wn_in_hz);
}
else {
    $jitter_longterm = $jitterRms * sqrt($vcoFreq/($wn_in_hz*$damping));
}
$dvc2_jlt = $icp*$jitter_longterm/$c2;
$freq_modulation_jlt = $kv * $dvc2_jlt;
$period_modulation_jlt = $period * ($freq_modulation_jlt/$vcoFreq);

$igate_phase_error = $tpfd * ($igate_leakage/$icp);
$dvc2_ig = $icp*$igate_phase_error/$c2;
$freq_modulation_ig = $kv * $dvc2_ig;
$period_modulation_ig = $period * ($freq_modulation_ig/$vcoFreq);

$estimated_phase_margin = (180/$PI) * (atan(2*$PI*$bandwidth_openloop_0dB * $rc1) 
    - atan(2*$PI*$bandwidth_openloop_0dB * $rc2) - atan(2*$PI*$bandwidth_openloop_0dB * $tau3)
    - atan(2*$PI*$bandwidth_openloop_0dB * $tau4) - 2*$PI*$bandwidth_openloop_0dB * $dly);

($thermal_noise_res,$thermal_noise_c1,$thermal_noise_c2) = &calc_thermal_noise($res,$c1,$c2,$bandwidth_openloop_0dB,$vcoFreq);    

printf STDOUT "\nplloop.pl rev <%s>. Copyright 2004, Dennis Fischette\n\n",$rev; 
printf STDOUT "Basic Loop Equations (Assumes 2nd-order loop with no feedback delay)\n";
printf STDOUT "------------------------------------------------------------------------------\n";
printf STDOUT "\nVCO period: %.3e sec\n", $period;
printf STDOUT "damping(D): %.2f\n", $damping;
printf STDOUT "natural frequency(wn): %.3e Hz\n", $wn_in_hz;
printf STDOUT "1/(Tau/2) (=2*D*wn): %.3e Hz\n", $inverse_time_constant_in_hz;
printf STDOUT "\tTau = loop time constant\n";
printf STDOUT "osr: %.3f\n", $osr;
printf STDOUT "\tosr = f(ref)/BW   where BW is unity-gain bandwidth of openloop transfer function\n";
printf STDOUT "trc2/tpfd %.2f\n", $rc2/$tpfd;
printf STDOUT "\tratio of R*C2 to reference period\n\n";

printf STDOUT "Gardner Discrete-Time Stability Limit (assumes 2nd-order loop - no C2 lpf cap)\n";
printf STDOUT "------------------------------------------------------------------------------\n";
printf STDOUT "Rlpf(max) = %.1f Ohm \t wn(max) = %.3e Hz\n\n", $res_upper_limit, $wn_upper_limit_in_hz;    

printf STDOUT "Equation-based phase margin estimate: accounts for all parasitic poles and delays\n";
printf STDOUT "----------------------------------------------------------------------------------\n";
printf STDOUT "estimated phase margin: %.2f degrees\n\n", $estimated_phase_margin;

printf STDOUT "RMS Thermal Noise on LPF resistors and capacitors (assumes f(vco) is upper limit of noise BW)\n";
printf STDOUT "--------------------------------------------------------------------------------------------\n";
printf STDOUT "Vrms(R) = %.3e V \t Vrms(C1) = %.3e V \t Vrms(C2) = %.3e V\n\n", 
               $thermal_noise_res,$thermal_noise_c1,$thermal_noise_c2;

if ($igate_leakage != 0.0) {
printf STDOUT "Estimated Phase Error due to Capacitor Leakage and Modulation due to Phase Error Correction:\n";
printf STDOUT "--------------------------------------------------------------------------------------------\n";
    printf STDOUT "gate-leakage static phase error: %.3e sec\n", $igate_phase_error;
    printf STDOUT "igate-modulated dV(c2): %.3e V\t\t-> assumes trc2 >> phase error\n", $dvc2_ig;
    printf STDOUT "VCO period modulation (igate): %.3e sec\t-> assumes trc2 >> phase error\n\n", $period_modulation_ig;
}

printf STDOUT "Estimated Accumulated Phase Error due to VCO Jitter and Modulation due to Phase Error Correction:\n";
printf STDOUT "-------------------------------------------------------------------------------------------------\n";
printf STDOUT "phase error: +/- %.3e sec\t\t\t-> due to accumulated random VCO jitter.\n", $jitter_longterm;
printf STDOUT "\tFor damping >= 1, jitter assumed to accumulate randomly for time=(1/wn_in_hz)\n";
printf STDOUT "\tFor damping < 1, jitter assumed to accumulate randomly for time=(1/(damping*wn_in_hz))\n";
printf STDOUT "dVctl(c2) modulation: %.3e V\t\t-> due to instantaneous response to phase error.\n", $dvc2_jlt;
printf STDOUT "\tAssumes trc2 >> phase error.\n";
printf STDOUT "VCO period modulation: +/- %.3e sec\t-> due to instantaneous response to phase error.\n", $period_modulation_jlt;
printf STDOUT "\tAssumes trc2 >> phase error.\n\n";

printf STDOUT "Closed-Loop Transfer Function Summary\n";
printf STDOUT "--------------------------------------------------------------------------------------------\n";
printf STDOUT "Peaking: %.3e dB at frequency\t %.2e Hz\n", $max_magnitude_closedloop_dB, $peaking_frequency; 
printf STDOUT "Bandwidth(0dB): %.3e Hz\n", $bandwidth_closedloop_0dB; 
printf STDOUT "Bandwidth(-3dB): %.3e Hz\n\n", $bandwidth_closedloop_minus3dB; 

printf STDOUT "Open-Loop Transfer Function Summary\n";
printf STDOUT "--------------------------------------------------------------------------------------------\n";
printf STDOUT "Bandwidth(0dB): %.3e Hz\n", $bandwidth_openloop_0dB; 
printf STDOUT "PhaseMargin: %.3e degrees\n", $phase_margin; 

close(OUTPUT);
close(OUTPUT_CLOSED_MAG);
close(OUTPUT_CLOSED_PHS);
close(OUTPUT_OPEN_MAG);
close(OUTPUT_OPEN_PHS);
exit(0);

sub add_c {
  my($a) = $_[0];
  my($b) = $_[1];

  my($real) = 0.0; 
  my($imag) = 0.0;

  $real = $a->[0] + $b->[0];
  $imag = $a->[1] + $b->[1];

  return($real, $imag);
  
}  ## sum_c


sub sub_c {
  my($a) = $_[0];
  my($b) = $_[1];

  my($real) = 0.0; 
  my($imag) = 0.0;

  $real = $a->[0] - $b->[0];
  $imag = $a->[1] - $b->[1];

  return($real, $imag);
  
}  ## sum_c


sub mul_c {
  my($a) = $_[0];
  my($b) = $_[1];

  my($real) = 0.0; 
  my($imag) = 0.0;

  $real = $a->[0]*$b->[0] - $a->[1]*$b->[1];
  $imag = $a->[0]*$b->[1] + $a->[1]*$b->[0];

  return($real, $imag);

}  ## mul_c


sub div_c {
  my($a) = $_[0];
  my($b) = $_[1];

  my($real) = 0.0; 
  my($imag) = 0.0;

  $real =       ($a->[0]*$b->[0] + $a->[1]*$b->[1])/($b->[0]*$b->[0]+$b->[1]*$b->[1]);
  $imag = (0.0 - $a->[0]*$b->[1] + $a->[1]*$b->[0])/($b->[0]*$b->[0]+$b->[1]*$b->[1]);

  return($real, $imag);

}  ## div_c


sub mag_c {
  my($a) = $_[0];

  return(sqrt($a->[0]*$a->[0] + $a->[1]*$a->[1]));

}  ## mag_c


sub ang_c {
  my($a) = $_[0];

  return(atan2($a->[1], $a->[0]));
 
}  ## ang_c

sub ang_c_degrees {
  my($a) = $_[0];

  return((180/$PI)*atan2($a->[1], $a->[0]));
 
}  ## ang_c

sub exp_c {
  my($a) = $_[0];

  my($real) = 0.0; 
  my($imag) = 0.0;

  $real = exp($a->[0]) * cos($a->[1]);
  $imag = exp($a->[0]) * sin($a->[1]);

  return($real, $imag);

}  ## exp_c

sub interpolate
{
    my($x1,$y1,$x2,$y2,$ytarg) = @_;
    my($xinterp,$fraction);
    $fraction = ($ytarg - $y1) / ($y2 - $y1);
    $xinterp = $x1 + $fraction * ($x2 - $x1);
    return($xinterp,$fraction);
}

sub calc_discrete_time_stability_limit
{
    my($r,$c1,$wn,$tref) = @_;    
    my($wref) = (2.0*$PI)/$tref;
    my($wn_upper_limit);
    my($wn_upper_limit_in_hz);
    my($res_upper_limit);

    # 2nd-order discrete-time assumptions - Gardner(1980), fig 4

    # wn^2 < wref^2 / (pi*(R*C1*wref + pi)
    $wn_upper_limit = sqrt(($wref*$wref) / ($PI*($r*$c1*$wref + $PI)));
    $wn_upper_limit_in_hz = $wn_upper_limit / (2.0 * $PI);

    # R < wref/(pi*C1*wn^2) - pi/(wref*RC1 + pi)
    $res_upper_limit = $wref/ ($PI*$c1*$wn*$wn) - $PI/($wref*$r*$c1 + $PI);

    return($wn_upper_limit_in_hz,$res_upper_limit);
}

sub calc_thermal_noise
{
    my($r,$c1,$c2,$wc_in_hz,$fvco) = @_;    
    my($dV_r,$dV_c1,$dV_c2);
    my($deltaFreq) = $fvco - $wc_in_hz;

    # capacitor thermal noise: V^2 = kT/C  Vrms^2
    $dV_c1 = sqrt($Kboltzmann * $temperature / $c1);
    $dV_c2 = sqrt($Kboltzmann * $temperature / $c2);

    # resistor thermal noise: V^2/f = 4kTR Vrms^2/Hz
    $dV_r = sqrt(4.0 * $Kboltzmann * $temperature * $r * $deltaFreq);     

    return($dV_r,$dV_c1,$dV_c2);
}