CN103595405B - 一种基于粒子滤波的锁相环实现方法 - Google Patents

Abstract

本发明涉及一种基于粒子滤波的锁相环实现方法，属于通信信号处理技术领域。本发明方法，具有比传统锁相环应用范围更加广泛的优势，可应用于非线性，非高斯的场合。

Description

一种基于粒子滤波的锁相环实现方法

技术领域

本发明涉及一种基于粒子滤波的锁相环实现方法，属于通信信号处理技术领域。

背景技术

锁相环技术应用于对宇宙中飞行目标的跟踪和遥测。锁相环以锁定输入载波信号的相位和频率为目标的一种载波环实现形式，然而传统锁相环中，环路滤波器采用的是低通滤波器，属于线性器件，只能解决线性问题，在本发明中采用粒子滤波器替代低通滤波器，可以使得锁相环的应用范围更加广泛

粒子滤波可以解决非线性，非高斯***的问题，因此应用范围广，日渐成为学术界的研究热点。近年来，粒子滤波逐渐引起信号处理和通信领域的重视。而本发明将粒子滤波方法应用于锁相环中来提高锁相环的性能，扩大锁相环的应用范围。

发明内容

本发明的目的是为改善现有锁相环不能应用于非线性非高斯的情况下的缺陷，提出一种基于粒子滤波的锁相环实现方法，扩大锁相环能够使用的范围。

本发明的目的是通过以下技术方案实现的。

本发明的一种基于粒子滤波的锁相环实现方法，其实现步骤如下：

步骤1、将输入的数字信号u_i输入进入鉴相器，鉴相器另一个端口输入的是上一时刻锁相环输出的信号，记为u_o。经鉴相器输出的结果为u_d。若为初始状态，则u_o可以为锁相环前端的器件粗略估计的信号。

步骤2、产生粒子。

按照设定的概率分布π（π通常选取一个均值为0，方差很大的Gaussian分布,可以依据具体情况而定，这种情况下，一般视为均匀分布）,采样得到k时刻粒子集i=1,2,...,N，其中上标i表示样本序号，i=1,2,…,N，下标k为符号周期序号，k=1,2,…。

步骤3、确定状态转移方程与观测方程。

粒子滤波的状态方程如式（1）所示：

x_k+1=Φx_k+Ψe_k+Γu_k（1）

其中， $x_k=\left[\begin{array}{c}{\stackrel{\·\·}{\mathrm{\θ}}}_k\\ {\stackrel{\·}{\mathrm{\θ}}}_k\end{array}\right],$ 一般x_k代表此时滤波器的状态； $\mathrm{\Φ}=\left[\begin{array}{cc}-1& 0\\ a_1& -1\end{array}\right],$ 其中a₁=T（T为环路滤波器更新时间间隔）； $\mathrm{\Ψ}=\left[\begin{array}{c}{a}_2\\ a_2+a_3\end{array}\right],$ 其中a₂=ω³ _np·T·k（ω_np为环路等效噪声的1.2倍，k为根据实际情况设定的系数），a₃=2ω² _np·T·k； $\mathrm{\Γ}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right];$ $u_k=\left[\begin{array}{c}{\stackrel{\·\·}{u}}_k\\ {\stackrel{\·}{u}}_k\end{array}\right]$ 和表示的是在环路滤波过程中产生的噪声。

粒子滤波的观测方程如式（2）所示：

$y_{k+1}=\left[\begin{array}{cc}0& 1\end{array}\right]{\hat{x}}_{k+1}+a_4{\hat{e}}_k+v_{k+1}---\left(2\right)$

其中，a₄=2ω_np·k；v_k+1表示在观测过程中所产生的噪声。

步骤4、粒子滤波初始化

步骤4.1根据式（1）k=0时，可以设x₀=0，将步骤2里产生的粒子集合里的所有粒子带入到方程（1）中，得到N个粒子；

步骤4.2，记i=1时产生的N个粒子样本值为,其中每一个粒子的重要性权值为

步骤4.3，将步骤4.1输出的每个粒子重要性权值进行归一化

$w_1^{\left(i\right)}=\frac{w_1^{*\left(i\right)}}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_1^{*\left(i\right)}}=\frac{1}{N}---\left(3\right)$

步骤4.4，计算估计值和

${\hat{x}}_1=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_1^{\left(i\right)}{w}_1^{\left(i\right)}---\left(4\right)$

${\hat{e}}_1=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\hat{e}}_1^{\left(i\right)}{w}_1^{\left(i\right)}---\left(5\right)$

步骤5、根据步骤4中得到的估计值和和观测观测方程得到这一时刻的输出值y₁。

步骤6、将步骤5的输出值输入至压控振荡器，将输出结果作为锁相环的结果，同时反馈给鉴相器。

步骤7、进行步骤1。

步骤8、粒子滤波：

步骤8.1、重复步骤2，得到粒子集合，根据式（1）可以将粒子集合里所有粒子带入到方程（1）中，得到N个粒子；

步骤8.2、计算第k时刻里，每一个粒子的重要性权值；

$w_k^{*\left(i\right)}=w_{k-1}^{*\left(i\right)}\exp [-\frac{1}{{2\mathrm{\σ}}^2}{(y_k^{\left(i\right)}-\left[\begin{array}{cc}0& 1\end{array}\right]x_k^i-k_4{e}_{k-1}^i)}^2]---\left(6\right)$

步骤8.3、权值归一化，当各个粒子权值计算完成后进行权值归一化；

${\mathrm{\ω}}_k^{\left(i\right)}={\mathrm{\ω}}_k^{*\left(i\right)}/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{*\left(i\right)}---\left(7\right)$

步骤8.4计算估计值和

${\hat{x}}_k=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_k^{\left(i\right)}{w}_k^{\left(i\right)}---\left(8\right)$

${\hat{e}}_k=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\hat{e}}_k^{\left(i\right)}{w}_k^{\left(i\right)}---\left(9\right)$

步骤9、根据步骤8中得到的估计值和和观测观测方程得到这一时刻的输出值y_k，并将结果输入至压控振荡器，将其结果作为锁相环的输出，并反馈给鉴相器，并重复执行步骤7～步骤9。

有益效果

本发明方法，具有比传统锁相环应用范围更加广泛的优势，可应用于非线性，非高斯的场合。

附图说明

图1是本发明的方法粒子滤波部分的冲击响应图。

具体实施方式

下面结合附图和实施例对本发明做进一步说明。

实施例

一种基于粒子滤波的锁相环实现方法，其实现步骤如下：

步骤1、将输入的数字信号u_i=sin(2π×20t)输入进入鉴相器，鉴相器另一个端口输入的是上一时刻锁相环输出的信号，记为u_o，作为初始时刻，没有上一时刻的输出信号u_o，这是可根据前端的收集到的信息作为u_o，我们在这里设u_o=cos(2π×19.5t)。经鉴相器输出的结果为

步骤2、产生粒子。

按照设定的概率分布π（π通常选取一个均值为0，方差很大的Gaussian分布,可以依据具体情况而定，我们在这里设定方差为1000）,采样得到k时刻粒子集i=1,2,...,N，其中上标i表示样本序号，i=1,2,…,N，这里N我们取值为100，即取1024个粒子；下标k为符号周期序号，k=1,2,…。

步骤3、确定状态转移方程与观测方程。

粒子滤波的状态方程如式（1）所示：

x_k+1=Φx_k+Ψe_k+Γu_k（1）

其中， $x_k=\left[\begin{array}{c}{\stackrel{\·\·}{\mathrm{\θ}}}_k\\ {\stackrel{\·}{\mathrm{\θ}}}_k\end{array}\right],$ 一般x_k代表此时滤波器的状态； $\mathrm{\Φ}=\left[\begin{array}{cc}-1& 0\\ a_1& -1\end{array}\right],$ 其中a₁=T（T为环路滤波器更新时间间隔在这里，我们设这个值为0.01ms）； $\mathrm{\Ψ}=\left[\begin{array}{c}{a}_2\\ a_2+a_3\end{array}\right],$ 其中a₂=ω³ _np·T·k（ω_np为环路等效噪声的1.2倍，在这里我们设环路等效噪声带宽为50Hz，则ω_np=60Hz，k为根据实际情况设定的系数，我们将k值设为1），a₃=2ω² _np·T·k； $\mathrm{\Γ}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right];$ $u_k=\left[\begin{array}{c}{\stackrel{\·\·}{u}}_k\\ {\stackrel{\·}{u}}_k\end{array}\right],$ 和表示的是在环路滤波过程中产生的噪声，我们设定和都服从高斯分布N（0,0.01²）。

粒子滤波的观测方程如式（2）所示：

$y_{k+1}=\left[\begin{array}{cc}0& 1\end{array}\right]{\hat{x}}_{k+1}+a_4{\hat{e}}_k+v_{k+1}---\left(2\right)$

其中，a₄=2ω_np·k；v_k+1表示在观测过程中所产生的噪声，我们设定v_k+1服从高斯分布N（0,0.01²）。

步骤4、粒子滤波初始化

步骤4.1根据式（1）k=0时，可以设x₀=0，将步骤2里产生的粒子集合里的所有粒子带入到方程（1）中，得到1024个粒子；

步骤4.2，记i=1时产生的1024个粒子样本值为,其中每一个粒子的重要性权值为

步骤4.3，将步骤4.1输出的每个粒子重要性权值进行归一化

$w_1^{\left(i\right)}=\frac{w_1^{*\left(i\right)}}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_1^{*\left(i\right)}}=\frac{1}{1024}---\left(3\right)$

步骤4.4，计算估计值和

${\hat{x}}_1=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_1^{\left(i\right)}{w}_1^{\left(i\right)}---\left(4\right)$

${\hat{e}}_1=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\hat{e}}_1^{\left(i\right)}{w}_1^{\left(i\right)}---\left(5\right)$

步骤5、根据步骤4中得到的估计值和和观测观测方程得到这一时刻的输出值 $y_1=\frac{1}{2}\sin (2\mathrm{\π}\×0.5t).$

步骤6、将步骤5的输出值输入至压控振荡器，将输出结果作为锁相环的结果，同时反馈给鉴相器输出结果。

步骤7、进行步骤1。

步骤8、粒子滤波。

步骤8.1、重复步骤2，得到粒子集合，根据式（1）可以将粒子集合里所有粒子带入到方程（1）中，得到1024个粒子；

步骤8.2、计算第k时刻里，每一个粒子的重要性权值；

$w_k^{*\left(i\right)}=w_{k-1}^{*\left(i\right)}\exp [-\frac{1}{{2\mathrm{\σ}}^2}{(y_k^{\left(i\right)}-\left[\begin{array}{cc}0& 1\end{array}\right]x_k^i-k_4{e}_{k-1}^i)}^2]---\left(6\right)$

步骤8.3、权值归一化，当各个粒子权值计算完成后进行权值归一化；

${\mathrm{\ω}}_k^{\left(i\right)}={\mathrm{\ω}}_k^{*\left(i\right)}/\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_k^{*\left(i\right)}---\left(7\right)$

步骤8.4计算估计值和

${\hat{x}}_k=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{x}_k^{\left(i\right)}{w}_k^{\left(i\right)}---\left(8\right)$

${\hat{e}}_k=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\hat{e}}_k^{\left(i\right)}{w}_k^{\left(i\right)}---\left(9\right)$

步骤9、根据步骤8中得到的估计值和和观测观测方程得到这一时刻的输出值y_k，并将结果输入至压控振荡器，将其结果作为锁相环的输出，并反馈给鉴相器，并重复执行步骤7～步骤9，最后会输出与输入端相位差一致的信号，结果为cos(2π×20t)。

Claims (1)

1.一种基于粒子滤波的锁相环实现方法，其特征在于步骤如下：

步骤1、将输入的数字信号u_i输入进入鉴相器，鉴相器另一个端口输入的是上一时刻锁相环输出的信号，记为u_o；经鉴相器输出的结果为u_d；若为初始状态，则u_o为锁相环前端的器件粗略估计的信号；

步骤2、产生粒子；

按照设定的概率分布π,π选取一个均值为0，方差很大的Gaussian分布,依据具体情况而定，这种情况下，视为均匀分布,采样得到k时刻粒子集其中上标i表示样本序号，i＝1,2,…,N，下标k为符号周期序号，k＝1,2,…；

步骤3、确定状态转移方程与观测方程；

粒子滤波的状态方程如式(1)所示：

x_k+1＝Φx_k+Ψe_k+Γu_k(1)

其中， $X_k=\left[\begin{array}{c}{\stackrel{\·\·}{\mathrm{\θ}}}_k\\ {\stackrel{\·}{\mathrm{\θ}}}_k\end{array}\right],$ x_k代表此时滤波器的状态； $\mathrm{\Φ}=\left[\begin{array}{cc}-1& 0\\ {\mathrm{\α}}_1& -1\end{array}\right],$ 其中a₁＝T,T为环路滤波器更新时间间隔； $\mathrm{\Ψ}=\left[\begin{array}{c}{a}_2\\ a_2+a_3\end{array}\right],$ 其中a₂＝ω³ _np·T·k,ω_np为环路等效噪声的1.2倍，k为根据实际情况设定的系数，a₃＝2ω² _np·T·k； $\mathrm{\Γ}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right];u_k=\left[\begin{array}{c}{\stackrel{\·\·}{u}}_k\\ {\stackrel{\·}{u}}_k\end{array}\right]$ 和表示的是在环路滤波过程中产生的噪声；

粒子滤波的观测方程如式(2)所示：

$y_{k+1}=\left[\begin{array}{cc}0& 1\end{array}\right]{\hat{x}}_{k+1}+a_4{\hat{e}}_k+v_{k+1}---\left(2\right)$

其中，a₄＝2ω_np·k；v_k+1表示在观测过程中所产生的噪声；

步骤4、粒子滤波初始化

步骤4.1根据式(1)k＝0时，设x₀＝0，将步骤2里产生的粒子集合里的所有粒子带入到方程(1)中，得到N个粒子；

步骤4.2，记i＝1时产生的N个粒子样本值为其中每一个粒子的重要性权值为

步骤4.3，将步骤4.1输出的每个粒子重要性权值进行归一化

$w_1^{\left(i\right)}=\frac{w_1^{*\left(i\right)}}{\underset{i=1}{\overset{N}{\Σ}}{w}_1^{*\left(i\right)}}=\frac{1}{N}---\left(3\right)$

步骤4.4，计算估计值和

${\hat{x}}_1=\underset{i=1}{\overset{N}{\Σ}}{x}_1^{\left(i\right)}{w}_1^{\left(i\right)}---\left(4\right)$

${\hat{e}}_1=\underset{i=1}{\overset{N}{\Σ}}{\hat{e}}_1^{\left(i\right)}{w}_1^{\left(i\right)}---\left(5\right)$

步骤5、根据步骤4中得到的估计值和和观测观测方程得到这一时刻的输出值y₁；

步骤6、将步骤5的输出值输入至压控振荡器，将输出结果作为锁相环的结果，同时反馈给鉴相器；

步骤7、进行步骤1；

步骤8、粒子滤波：

步骤8.1、重复步骤2，得到粒子集合，根据式(1)将粒子集合里所有粒子带入到方程(1)中，得到N个粒子；

步骤8.2、计算第k时刻里，每一个粒子的重要性权值；

$w_k^{*\left(i\right)}=w_{k-1}^{*\left(i\right)}\exp \[-\frac{1}{2{\mathrm{\σ}}^2}{\left(y_k^{\left(i\right)}-\left[\begin{array}{cc}0& 1\end{array}\right]x_k^i-k_4{e}_{k-1}^i\right)}^2\]---\left(6\right)$

步骤8.3、权值归一化，当各个粒子权值计算完成后进行权值归一化；

${\mathrm{\ω}}_k^{\left(i\right)}={\mathrm{\ω}}_k^{*\left(i\right)}/\underset{i=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_k^{*\left(i\right)}---\left(7\right)$

步骤8.4计算估计值和

${\hat{x}}_k=\underset{i=1}{\overset{N}{\Σ}}{x}_k^{\left(i\right)}{w}_k^{\left(i\right)}---\left(8\right)$

${\hat{e}}_k=\underset{i=1}{\overset{N}{\Σ}}{\hat{e}}_k^{\left(i\right)}{w}_k^{\left(i\right)}---\left(9\right)$

步骤9、根据步骤8中得到的估计值和和观测观测方程得到这一时刻的输出值y_k，并将结果输入至压控振荡器，将其结果作为锁相环的输出，并反馈给鉴相器，并重复执行步骤7～步骤9。