def PLL_cbb(x,fs,loop_type,Kv,fn,zeta):
T = <NUM_LIT:1>/float(fs)<EOL>Kv = <NUM_LIT:2>*np.pi*Kv <EOL>if loop_type == <NUM_LIT:1>:<EOL><INDENT>K = <NUM_LIT:2>*np.pi*fn <EOL><DEDENT>elif loop_type == <NUM_LIT:2>:<EOL><INDENT>K = <NUM_LIT:4> *np.pi*zeta*fn <EOL>tau2 = zeta/(np.pi*fn)<EOL><DEDENT>elif loop_type == <NUM_LIT:3>:<EOL><INDENT>K = Kv <EOL>tau1 = K/((<NUM_LIT:2>*np.pi*fn)^<NUM_LIT:2>);<EOL>tau2 = <NUM_LIT:2>*zeta/(<NUM_LIT:2>*np.pi*fn)*(<NUM_LIT:1> - <NUM_LIT:2>*np.pi*fn/K*<NUM_LIT:1>/(<NUM_LIT:2>*zeta))<EOL><DEDENT>else:<EOL><INDENT>print('<STR_LIT>')<EOL><DEDENT>filt_in_last = <NUM_LIT:0>; filt_out_last = <NUM_LIT:0>;<EOL>vco_in_last = <NUM_LIT:0>; vco_out = <NUM_LIT:0>; vco_out_last = <NUM_LIT:0>;<EOL>vco_out_cbb = <NUM_LIT:0><EOL>n = np.arange(len(x))<EOL>theta_hat = np.zeros(len(x))<EOL>ev = np.zeros(len(x))<EOL>phi = np.zeros(len(x))<EOL>for k in  range(len(n)):<EOL><INDENT>phi[k] = np.imag(x[k] * np.conj(vco_out_cbb))<EOL>pd_out = phi[k]<EOL>gain_out = K/Kv*pd_out <EOL>if loop_type == <NUM_LIT:2>:<EOL><INDENT>filt_in = (<NUM_LIT:1>/tau2)*gain_out<EOL>filt_out = filt_out_last + T/<NUM_LIT:2>*(filt_in + filt_in_last)<EOL>filt_in_last = filt_in<EOL>filt_out_last = filt_out<EOL>filt_out = filt_out + gain_out<EOL><DEDENT>elif loop_type == <NUM_LIT:3>:<EOL><INDENT>filt_in = (tau2/tau1)*gain_out - (<NUM_LIT:1>/tau1)*filt_out_last<EOL>u3 = filt_in + (<NUM_LIT:1>/tau2)*filt_out_last<EOL>filt_out = filt_out_last + T/<NUM_LIT:2>*(filt_in + filt_in_last)<EOL>filt_in_last = filt_in<EOL>filt_out_last = filt_out<EOL><DEDENT>else:<EOL><INDENT>filt_out = gain_out;<EOL><DEDENT>vco_in = filt_out<EOL>if loop_type == <NUM_LIT:3>:<EOL><INDENT>vco_in = u3<EOL><DEDENT>vco_out = vco_out_last + T/<NUM_LIT:2>*(vco_in + vco_in_last)<EOL>vco_in_last = vco_in<EOL>vco_out_last = vco_out<EOL>vco_out = Kv*vco_out <EOL>vco_out_cbb = np.exp(<NUM_LIT>*vco_out)<EOL>ev[k] = vco_in<EOL>theta_hat[k] = vco_out<EOL><DEDENT>return theta_hat, ev, phi<EOL>
Baseband Analog PLL Simulation Model

:param x: input phase deviation in radians
:param fs: sampling rate in sample per second or Hz
:param loop_type: 1, first-order loop filter F(s)=K_LF; 2, integrator
            with lead compensation F(s) = (1 + s tau2)/(s tau1),
            i.e., a type II, or 3, lowpass with lead compensation
            F(s) = (1 + s tau2)/(1 + s tau1)
:param Kv: VCO gain in Hz/v; note presently assume Kp = 1v/rad
            and K_LF = 1; the user can easily change this
:param fn: Loop natural frequency (loops 2 & 3) or cutoff
            frequency (loop 1)
:param zeta: Damping factor for loops 2 & 3
:return: theta_hat = Output phase estimate of the input theta in radians,
         ev = VCO control voltage,
         phi = phase error = theta - theta_hat

Mark Wickert, April 2007 for ECE 5625/4625
Modified February 2008 and July 2014 for ECE 5675/4675
Python version August 2014