一个简单的python数字锁相环

一个简单的数字锁相环demo：

import matplotlib
matplotlib.use('Qt5Agg')

import numpy as np
import matplotlib.pyplot as plt

#parameters
phase_offset      = 0.00    #carrier phase offset
frequency_offset  = 0.30    #carrier frequency offset
wn                = 0.01    #pll bandwidth
zeta              = 0.707   #pll damping factor
K                 = 1000    # pll loop gain
n                 = 400   #number of samples

#generate loop filter parameters (active PI design)
t1 = K/(wn*wn)    #tau_1
t2 = 2*zeta/wn    #tau_2

#feed-forward coefficients (numerator)
b0 = (4*K/t1)*(1.+t2/2.0)
b1 = (8*K/t1)
b2 = (4*K/t1)*(1.-t2/2.0)

#feed-back coefficients (denominator)
a1 = -2.0
a2 =  1.0

#print filter coefficients (as comments)
print b0,b1,b2
print a1,a2

#filter buffer
v0=0.0
v1=0.0
v2=0.0

#initialize states
phi     = phase_offset   #input signal's initial phase
phi_hat = 0.0            #PLL's initial phase
n =500

delta_phi = np.zeros(n)

#run basic simulation
for i in range(n):

    #compute input sinusoid and update phase
    x = np.exp(1j*phi)
    phi += frequency_offset

    #compute PLL output from phase estimate
    y = np.exp(1j * phi_hat)

    #compute error estimate
    delta_phi[i] = np.angle(x * np.conjugate(y))#cargf( x * conjf(y) )carg is a standard library function which calculates the argument (phase angle) of a complex number.

    #push result through loop filter, updating phase estimate

    #advance buffer
    v2 = v1  #shift center register to upper register
    v1 = v0  #shift lower register to center register

    #compute new lower register
    v0 = delta_phi[i] - v1*a1 - v2*a2

    #compute new output
    phi_hat = v0*b0 + v1*b1 + v2*b2

    phi_hat = np.mod(phi_hat, 2 * np.pi)

print delta_phi

plt.plot(delta_phi)
plt.grid()
plt.show()

可以看出迭代近200次后收敛。