CN103064093A - 一种gps接收机中ldpc码辅助的迭代载波同步方法 - Google Patents

Abstract

一种GPS接收机中LDPC码辅助的迭代载波同步方法，属数字通信中信道编码及载波同步技术领域。LDPC码的理论性能很好，但是它对同步的要求很苛刻，在有频偏及相偏的条件下，其性能会急剧恶化。本发明方法联合译码和估计，在LDPC译码器中增加了对载波同步参数的估计,并用译码中产生的软信息估计参数。译码器每进行一次译码迭代，便对参数进行一次估计，估计时要用到LDPC译码中产生的软信息，然后再用估计的参数对译码器的输入进行修正，修正的输入经过LDPC码译码处理又可得到新的软信息，这些软信息又可用来得到新的估计参数，通过多次迭代，该方法可使参数的估计值收敛到实际取值，从而使LDPC码接近理论上的性能。

Description

一种GPS接收机中LDPC码辅助的迭代载波同步方法

技术领域

本发明涉及一种GPS接收机中LDPC码辅助的迭代载波同步方法，属数字通信中信道编码及载波同步技术领域。

背景技术

GPS是美国军方研发的卫星导航系统，可以在各种天气下精确地提供地表附近任意位置精确的时间和位置信息。该系统由空间段，控制段和用户段三部分组成。空间段至少包括24颗卫星。卫星使用码分复用(CDMA)技术在两个频率上广播测距码和导航数据，即GPS系统只使用两个频率：L1(1275.42MHz)和L2(1227.6MHz)。导航数据使GPS接收机能够确定卫星发送信号时的位置，测距码使接收机能够计算出信号传输时间。接收机根据卫星发送信号的位置和信号传输到接收机的时间就可计算出当前所在位置。

IS-GPS-800规范定义了L1链路上的民用L1信号(L1C)。L1C信号由L1C_P和L1C_D两部分组成。L1C_D又分为子帧，帧和超帧。帧分成3个长度不等的子帧。为了向用户广播完整的消息，需发送多个帧，即超帧。子帧2和子帧3分别包含600比特和274数据。两个子帧都要用码率为1/2的LDPC编码器编码，得到总长度为1748个符号序列。这1748个符号接着送到交织器进行交织，1748个符号按行写入38行46的矩阵，然后按列读出。交织后的输出进行二进制偏移载波调制(BOC(1,1))，即先以1.023Mbps的码片速率进行直接序列扩频(DSSS)，再用1.023MHz的方波副载波调制扩频后的信号，副载波的作用是改变发射信号的频谱。最后再将基带BOC(1,1)信号通过BPSK调制到L1射频链路上传输。

LDPC(Low-Density Parity-Check)码是一种纠错能力极好的线性分组码，目前已被很多标准采用，比如IEEE802.16e、IEEE802.11n、DVB-S2和GPS。虽然理论上LDPC码具有接近香农限的性能，但是这是在理想相干检测的前提下得出来的，实际应用中载波同步不理想会引进频偏和相偏，而LDPC码译码算法对频偏和相偏很敏感，这使得实际系统中LDPC码的性能和理论性能有很大的差距。从理论上分析相偏和频偏对LDPC码译码器的影响很困难，目前还没有这方面的文献，但是通过仿真可以知道，相偏和频偏对LDPC译码器的影响很严重。

鉴于LDPC码对载波同步不理想引起的频偏和相偏很敏感，GPS中提高载波同步的精度对整个系统的性能都有很大的提高。

发明内容

为了克服载波同步不理想引起的频偏和相偏对LDPC码迭代译码器的影响，本发提供了一种GPS接收机中LDPC码辅助的迭代载波同步方法。该方法将译码和同步相结合，在标准LDPC译码的基础上，增加了对同步参数估计的步骤。在每次迭代译码中，用LDPC译码中的判决消息估计相偏和频偏，再用估计的相偏和频偏对译码器变量节点的对数似然比消息进行修正，通过LDPC译码算法和同步算法的相互影响、相互修正达到精确同步、提高译码性能的目的。

为了实现上述发明目的，本发明采用的技术方案如下：

一种GPS接收机中LDPC码辅助的迭代载波同步方法，应用于GPS接收机中的LDPC码译码器，以实现精确估计同步不理想引起的频偏和相偏，提高整个系统性能的目的；预先设C是由m×n维校验矩阵H＝{h_j,k}确定的长为n的二元LDPC码；R_j表示同第j个校验节点相连的变量节点的集合，即R_j＝{k|h_k,j＝1}，R_j\k表示从R_j中除去第k个变量节点的集合；C_k表示同第k个变量节点相连的校验节点的集合，即C_k＝{j|h_k,j＝1}，C_k\j表示从C_k中除去第j个校验节点的集合；a表示码C中的码字经过BPSK调制后的发送信息序列，r表示a经过信道后的接收信息序列，也即LDPC译码器的输入，r中第k个分量对应码C的第k个变量节点；p_k(0)和p_k(1)分别表示第k个变量节点为0和1的概率，L(P_k)表示第k个变量节点的初始对数似然比消息；该方法的步骤如下：

1）初始化

假设BPSK将0映射到-1、1映射到1，则发送序列经过均值为0方差为σ²白高斯噪声信道后

r_k是r的第k个分量，所以第k个变量节点的初始对数似然比消息为

$L\left(P_k\right)=\ln \frac{p_k\left(0\right)}{p_k\left(1\right)}={|r_k-1|}^2-{|r_k+1|}^2$

2）校验节点消息更新

校验节点收集与它相邻的变量节点的消息；LDPC码译码器采用偏移最小和算法；当LDPC码译码器使用最小和译码时，令L'(r_jk)表示第j个校验节点传给第k个变量节点的对数似然比消息；最小和译码算法中校验节点消息更新的公式如下：

$L^{\text{\′}}\left(r_{\mathrm{jk}}\right)=\underset{k^{\′}\∈R_k\backslash j}{\mathrm{\Π}}\mathrm{sign}\left(\text{L}\left(q_{k^{\′}j}\right)\right)\·\underset{k^{\′}\∈R_k\backslash j}{\min }\left(\right|L\left(q_{k^{\′}j}\right)\left|\right)$

其中 $\mathrm{sign}\left(x\right)=\left\{\begin{array}{cc}1,& x\&GreaterEqual;0\\ 0,& x<0\end{array}\right.$ 是符号函数，符号min(·)表示求最小值，L(q_k′j)表示第k′个变量节点传给第j个校验节点的对数似然比消息，符号k'∈R_k\j表示k'是集合R_k\j的元素；在最小和算法中引进偏移因子β，便得到偏移最小和算法，其中β的取值介于0和1之间，用L(r_jk)表示偏移最小和算法中第j个校验节点传给第k个变量节点的对数似然比消息；偏移最小和算法中校验节点消息更新的公式如下：

L(r_jk)＝sign(L'(r_jk))·max(|L'(r_jk)|-β，0)

其中符号max(·)表示求最大值；不同LDPC码对应的偏移因子取值不同，GPS系统使用了两种LDPC码，码长分别为1200和578，这两种码的偏移因子分别通过仿真来求得，即让偏移因子在0和1间取不同值，根据仿真结果选择使误码率最低的那个偏移因子；

3）计算判决消息并做判决

变量节点收集与它相邻的校验节点的消息及来自信道的初始消息；用L(q_k)表示第k个变量节点的判决消息，则

$L\left(q_k\right)=L\left(P_k\right)+\underset{j^{\&prime;}\text{\&Element;}{C}_k}{\mathrm{\&Sigma;}}L\left(r_{j^{\&prime;}k}\right)$

如果L(q_k)＞0，则判定

否则判定

表示码字中第k个变量节点的判决输出；

4）判断译码是否结束

满足以下条件之一即表示译码结束：

A．

H是LDPC码的校验矩阵，

表示

的转置；

B．已达到指定的最大迭代次数，

否则转入下一步；

5）若译码没有结束，则对同步参数进行一次估计；

用b＝[Δf,Δθ]表示同步不理想引起的频偏和相偏，其中Δf和Δθ分别表示频偏和相偏，使用最大对数似然概率的方法估计b，并用期望-最大（EM）算法求解该最大对数似然概率估计问题，EM迭代求解中要用到LDPC译码中产生的软信息；

b的最大对数似然概率估计是指找到使对数似然概率密度函数lnp(r|b)最大的b作为其估计值，记作

即

符号

表示求使括号内参数取最大值的b，EM算法通过引进先验分布已知的隐藏变量来达到简化计算的目的；取发送信息序列a作为先验分布已知的随机变量，则标准EM算法求解最大对数似然估计的迭代计算公式表示为

$b^{(n+1)}=\underset{b}{\arg \max }\left\{\underset{a}{\mathrm{\&Sigma;}}p\left(a\right|r,b^{\left(n\right)})\ln p\left(r\right|a,b)\right\}---\left(1\right)$

其中b⁽ⁿ⁾表示第n次迭代中的b估计值，b⁽ⁿ⁺¹⁾表示在b的当前估计值为b⁽ⁿ⁾时下一次的估计值，p(r|a,b)表示已知a和b的条件下r的条件概率密度函数，p(a|r,b⁽ⁿ⁾)表示已知a和b⁽ⁿ⁾的条件下r的条件概率密度函数，

是求和符号；对于GPS系统，在扩频序列和方波副载波同步理想，即接收端可以精确恢复它们的定时信息的情况下，整个系统可简化为BPSK调制系统，在已知发送信息序列a和同步参数b＝[Δf,Δθ]的条件下，译码器输入r的对数条件概率表示为

$\ln p\left(r\right|a,b)=-\frac{1}{{2\mathrm{\&sigma;}}^2}\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{|a_k{e}^{j(2\mathrm{\&pi;\&Delta;fkT}+\mathrm{\&Delta;\&theta;})}-r_k|}^2$ $\left(2\right)$

$=-\frac{1}{2{\mathrm{\&sigma;}}^2}\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}({\left|a_k\right|}^2+{\left|r_k\right|}^2)+\frac{1}{{\mathrm{\&sigma;}}^2}\mathrm{Re}\left\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{a}_k^{*}{e}^{-j(2\mathrm{\&pi;\&Delta;fkT}+\mathrm{\&Delta;\&theta;})}\right\}$

其中L是码字的长度，T是发送符号的持续时间，a_k和r_k分别是发送和接收信息序列中第k个分量的取值，σ²是噪声方差，Re{·}表示取复数的实部，

表示a_k的复共轭；公式（2）分成两部分，前一部分与估计参数b＝[Δf,Δθ]无关，它只影响公式（1）所能取的最大值，而不能影响使（1）式取得最大值的参数b，所以这一部分忽略；同理，公式（2）后一部分的比例因子也去掉，所以对（1）式求解简化为

$b^{(n+1)}=\underset{b}{\arg \max }\left\{\underset{a}{\mathrm{\&Sigma;}}p\left(a\right|r,b^{\left(n\right)})\mathrm{Re}\right\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{a}_k^{*}{e}^{-j(2\mathrm{\&pi;\&Delta;fkT}+\mathrm{\&Delta;\&theta;})}\left\}\right\}$ $\left(3\right)$

$=\underset{b}{\arg \max }\left\{\mathrm{Re}\right\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\left[a_k^{*}\right|r,b^{\left(n\right)}\left]e^{-j2\mathrm{\&pi;\&Delta;fkT}}{e}^{-\mathrm{j\&Delta;\&theta;}}\right\}\}$

其中

表示在给定r和b⁽ⁿ⁾的条件下

的条件期望；复数实部的最大值不能超过其模的最大值，且只有复数为实数时，实部的最大值和模的最大值等价，求公式（3）复数实部的最大值分两步：

A、求出使复数

模取最大值的Δf；

B、Δθ取复数 $\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\left[a_k^{*}\right|r,b^{\left(n\right)}]e^{-j2\mathrm{\&pi;\&Delta;fkT}}$ 的幅角，即旋转 $\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\left[a_k^{*}\right|r,b^{\left(n\right)}]e^{-j2\mathrm{\&pi;\&Delta;fkT}}$ 使其与实轴重合，此时模的最大值也就是实部的最大值；

根据上述分析，EM迭代算法中频偏及相偏的估计公式表示为

$\text{\&Delta;}\hat{f}=\underset{\mathrm{\&Delta;f}}{\max }\left|\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\right[a_k^{*}|r,b^{\left(n\right)}]e^{-j2\mathrm{\&pi;\&Delta;fkT}}|$

$\text{\&Delta;}\widehat{\mathrm{\&theta;}}=\mathrm{angle}\left\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\right[a_k^{*}|r,b^{\left(n\right)}]e^{-j2\mathrm{\&pi;\&Delta;}\hat{f}\mathrm{kT}}\}$

其中angle{·}表示取复数的幅角，

和分别表示第n次EM迭代时Δf和Δθ的估计值，即对于BPSK调制，a中各信息比特取1或-1，即a_k＝1或a_k＝-1，因此

$E_a\left[a_k^{*}\right|r,b^{\left(n\right)}]=p(a_k=1|r,b^{\left(n\right)})-p(a_k=-1|r,b^{\left(n\right)})---\left(4\right)$

p(a_k＝1|r,b⁽ⁿ⁾)和p(a_k＝-1|r,b⁽ⁿ⁾)分别表示已知r和b⁽ⁿ⁾的条件下a_k取1和-1的概率，所以估计Δf和Δθ需知道a中信息比特a_k的联合统计p(a_k＝1|r,b⁽ⁿ⁾)和p(a_k＝-1|r,b⁽ⁿ⁾)，当迭代次数n足够大时，第n次迭代中产生的判决信息L(q_k)和a_k的联合统计满足如下关系

$L\left(q_k\right)\text{=ln}\frac{p(a_k=-1|r,b^{\left(n\right)})}{p(a_k=1|r,b^{\left(n\right)})}$

根据上式及p(a_k＝-1|r,b⁽ⁿ⁾)+p(a_k＝1|r,b⁽ⁿ⁾)=1得到联合统计

$p(a_k=1|r,b^{\left(n\right)})=\frac{e^{L\left(q_k\right)}}{1+e^{L\left(q_k\right)}},$ $p(a_k=-1|r,b^{\left(n\right)})=\frac{1}{1+e^{L\left(q_k\right)}}$

代入公式(4)得到

$E_a\left[a_k^{*}\right|{r,b}^{\left(n\right)}]=\frac{e^{L\left(q_k\right)}-1}{e^{L\left(q_k\right)}+1}$

上式就把LDPC迭代译码中获得的判决信息L(q_k)应用到了同步参数频偏和相偏的估计；上述同步参数频偏和相偏的估计分两步进行，即：

A、估计频偏

$\mathrm{\&Delta;}\hat{f}=\underset{\mathrm{\&Delta;f}}{\max }|\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k\&CenterDot;\frac{e^{L\left(q_k\right)}-1}{e^{L\left(q_k\right)}+1}\&CenterDot;e^{-j2\mathrm{\&pi;\&Delta;fkT}}|$

B、估计相偏

$\mathrm{\&Delta;}\widehat{\mathrm{\&theta;}}=\mathrm{angle}\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k\&CenterDot;\frac{e^{L\left(q_k\right)}-1}{e^{L\left(q_k\right)}+1}\&CenterDot;e^{-j2\mathrm{\&pi;\&Delta;}\hat{f}\mathrm{kT}}\}$

6）修正接收信息序列r

用当前频偏和相偏估计值对接收信息序列r进行修正；用r′表示修正后的接收信息序列，其第k个分量用r′_k表示，r_k表示r中第k个分量；修正公式为

$r_k^{\&prime;}=r_k{e}^{-j(\mathrm{\&Delta;}\widehat{\mathrm{\&theta;}}+2\mathrm{\&pi;\&Delta;}\hat{f}\mathrm{kT})}$

7）计算变量节点初始对数似然比消息

用修正后的接收信息序列r′计算各变量节点的初始对数似然比消息；在信道噪声方差为σ²的情况下，

所以第k个变量节点的初始对数似然比消息为

$L\left(P_k\right)=\ln \frac{p_k\left(0\right)}{p_k\left(1\right)}={|r_k^{\&prime;}-1|}^2-{|r_k^{\&prime;}+1|}^2$

8）变量节点消息更新

变量节点收集与它相邻的校验节点及信道的消息，

$L\left(q_{\mathrm{kj}}\right)=L\left(P_k\right)+\underset{j^{\&prime;}\&Element;C_k\backslash j}{\mathrm{\&Sigma;}}\&CenterDot;L\left(r_{j^{\&prime;}k}\right)$

其中L(q_kj)表示第k个变量节点传给第j个校验节点的对数似然比消息，L(P_k)是来自信道的初始对数似然比消息，Lr_j′k)是第j′个校验节点传给第k个变量节点的对数似然比消息，j'∈C_k\j表示j'是集合C_k\j的元素；

9）返回到步骤2)更新校验节点消息，进行下一次迭代。

上述BPSK全称为Binary Phase Shift Keying，意为二进制移相键控法。它是把模拟信号转换成数据值的转换方式之一。是利用偏离相位的复数波浪组合来表现信息键控移相方式的一种。

上述LDPC码是一种线性分组奇偶校验码,是麻省理工学院Robert Gal lager于1962年在博士论文中提出的一种具有稀疏校验矩阵的分组纠错码，具有译码复杂度低、结构灵活、性能优良等特点。GPS系统使用了两种LDPC码，码长分别为1200和548，这两种LDPC码的校验矩阵可在“Global Positioning Systems Wing(GPSW)System Engineering&IntegrationInterface Specification IS-GPS-800Navstar GPS Space Segment/User Segment L1CInterfaces”中找到。

上述EM算法是最大期望算法（Expectation-maximization algorithm，又译为期望最大化算法），是在概率模型中寻找参数最大似然估计或者最大后验估计的算法，其中概率模型依赖于无法观测的隐藏变量。最大期望经常用在机器学习和计算机视觉的数据聚类领域。

本发明方法的原理如下:该方法以迭代的方式对同步参数进行估计，并将迭代过程和LDPC的译码相结合，即每进行一次LDPC译码迭代的同时，也进行一次参数估计，估计参数时要用到LDPC译码时产生的软信息，估计出的参数再对译码器的输入进行修正，LDPC译码器又根据修正后的输入重新计算信道的初始对数似然比消息，通过这种互相影响的方式，达到用LDPC码辅助载波同步的目的。

本发明方法的优点如下：

本发明使用LDPC迭代译码中产生的软信息来辅助GPS系统中的载波同步，将对同步参数的估计放在LDPC码的译码器，从而可以在解调器不需做任何修改的条件下和任何传统同步算法相结合，对传统同步算法估计的频偏和相偏进行精细估计。

具体实施方式

下面结合实施例对本发明作进一步说明，但不限于此。

实施例：

一种GPS接收机中LDPC码辅助的迭代载波同步方法，应用于GPS接收机中的LDPC码译码器，以实现精确估计同步不理想引起的频偏和相偏，提高整个系统性能的目的；预先设C是由m×n维校验矩阵H＝{h_j,k}确定的长为n的二元LDPC码；R_j表示同第j个校验节点相连的变量节点的集合，即R_j＝{k|h_k,j＝1}，R_j\k表示从R_j中除去第k个变量节点的集合；C_k表示同第k个变量节点相连的校验节点的集合，即C_k＝{j|h_k,j＝1}，C_k\j表示从C_k中除去第j个校验节点的集合；a表示码C中的码字经过BPSK调制后的发送信息序列，r表示a经过信道后的接收信息序列，也即LDPC译码器的输入，r中第k个分量对应码C的第k个变量节点；p_k(0)和p_k(1)分别表示第k个变量节点为0和1的概率，L(P_k)表示第k个变量节点的初始对数似然比消息；该方法的步骤如下：

1）初始化

假设BPSK将0映射到-1、1映射到1，则发送序列经过均值为0方差为σ²白高斯噪声信道后

r_k是r的第k个分量，所以第k个变量节点的初始对数似然比消息为

$L\left(P_k\right)=\ln \frac{p_k\left(0\right)}{p_k\left(1\right)}={|r_k-1|}^2-{|r_k+1|}^2$

2）校验节点消息更新

校验节点收集与它相邻的变量节点的消息；LDPC码译码器采用偏移最小和算法；当LDPC码译码器使用最小和译码时，令L'(r_jk)表示第j个校验节点传给第k个变量节点的对数似然比消息；最小和译码算法中校验节点消息更新的公式如下：

$L^{\text{\&prime;}}\left(r_{\mathrm{jk}}\right)=\underset{k^{\&prime;}\&Element;R_k\backslash j}{\mathrm{\&Pi;}}\mathrm{sign}\left(\text{L}\left(q_{k^{\&prime;}j}\right)\right)\&CenterDot;\underset{k^{\&prime;}\&Element;R_k\backslash j}{\min }\left(\right|L\left(q_{k^{\&prime;}j}\right)\left|\right)$

其中 $\mathrm{sign}\left(x\right)=\left\{\begin{array}{cc}1,& x\&GreaterEqual;0\\ 0,& x<0\end{array}\right.$ 是符号函数，符号min(·)表示求最小值，L(q_k′j)表示第k′个变量节点传给第j个校验节点的对数似然比消息，符号k'∈R_k\j表示k'是集合R_k\j的元素；在最小和算法中引进偏移因子β，便得到偏移最小和算法，其中β的取值介于0和1之间，用L(r_jk)表示偏移最小和算法中第j个校验节点传给第k个变量节点的对数似然比消息；偏移最小和算法中校验节点消息更新的公式如下：

L(r_jk)＝sign(L'(r_jk))·max(|L'(r_jk)|-β，0)

其中符号max(·)表示求最大值；不同LDPC码对应的偏移因子取值不同，GPS系统使用了两种LDPC码，码长分别为1200和578，这两种码的偏移因子分别通过仿真来求得，即让偏移因子在0和1间取不同值，根据仿真结果选择使误码率最低的那个偏移因子；

3）计算判决消息并做判决

变量节点收集与它相邻的校验节点的消息及来自信道的初始消息；用L(q_k)表示第k个变量节点的判决消息，则

$L\left(q_k\right)=L\left(P_k\right)+\underset{j^{\&prime;}\text{\&Element;}{C}_k}{\mathrm{\&Sigma;}}L\left(r_{j^{\&prime;}k}\right)$

如果L(q_k)＞0，则判定

否则判定

表示码字中第k个变量节点的判决输出；

4）判断译码是否结束

满足以下条件之一即表示译码结束：

A．H是LDPC码的校验矩阵，

表示

的转置；

B．已达到指定的最大迭代次数，

否则转入下一步；

5）若译码没有结束，则对同步参数进行一次估计；

用b＝[Δf,Δθ]表示同步不理想引起的频偏和相偏，其中Δf和Δθ分别表示频偏和相偏，使用最大对数似然概率的方法估计b，并用期望-最大（EM）算法求解该最大对数似然概率估计问题，EM迭代求解中要用到LDPC译码中产生的软信息；

b的最大对数似然概率估计是指找到使对数似然概率密度函数ln p(r|b)最大的b作为其估计值，记作

即

符号

表示求使括号内参数取最大值的b，EM算法通过引进先验分布已知的隐藏变量来达到简化计算的目的；取发送信息序列a作为先验分布已知的随机变量，则标准EM算法求解最大对数似然估计的迭代计算公式表示为

$b^{(n+1)}=\underset{b}{\arg \max }\left\{\underset{a}{\mathrm{\&Sigma;}}p\left(a\right|r,b^{\left(n\right)})\ln p\left(r\right|a,b)\right\}---\left(5\right)$

其中b⁽ⁿ⁾表示第n次迭代中的b估计值，b⁽ⁿ⁺¹⁾表示在b的当前估计值为b⁽ⁿ⁾时下一次的估计值，p(r|a,b)表示已知a和b的条件下r的条件概率密度函数，p(a|r,b⁽ⁿ⁾)表示已知a和b⁽ⁿ⁾的条件下r的条件概率密度函数，

是求和符号；对于GPS系统，在扩频序列和方波副载波同步理想，即接收端可以精确恢复它们的定时信息的情况下，整个系统可简化为BPSK调制系统，在已知发送信息序列a和同步参数b＝[Δf,Δθ]的条件下，译码器输入r的对数条件概率表示为

$\ln p\left(r\right|a,b)=-\frac{1}{{2\mathrm{\&sigma;}}^2}\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{|a_k{e}^{j(2\mathrm{\&pi;\&Delta;fkT}+\mathrm{\&Delta;\&theta;})}-r_k|}^2$ $\left(6\right)$

$=-\frac{1}{2{\mathrm{\&sigma;}}^2}\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}({\left|a_k\right|}^2+{\left|r_k\right|}^2)+\frac{1}{{\mathrm{\&sigma;}}^2}\mathrm{Re}\left\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{a}_k^{*}{e}^{-j(2\mathrm{\&pi;\&Delta;fkT}+\mathrm{\&Delta;\&theta;})}\right\}$

其中L是码字的长度，T是发送符号的持续时间，a_k和r_k分别是发送和接收信息序列中第k个分量的取值，σ²是噪声方差，Re{·}表示取复数的实部，

表示a_k的复共轭；公式（2）分成两部分，前一部分与估计参数b＝[Δf,Δθ]无关，它只影响公式（1）所能取的最大值，而不能影响使（1）式取得最大值的参数b，所以这一部分忽略；同理，公式（2）后一部分的比例因子

也去掉，所以对（1）式求解简化为

$b^{(n+1)}=\underset{b}{\arg \max }\left\{\underset{a}{\mathrm{\&Sigma;}}p\left(a\right|r,b^{\left(n\right)})\mathrm{Re}\right\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{a}_k^{*}{e}^{-j(2\mathrm{\&pi;\&Delta;fkT}+\mathrm{\&Delta;\&theta;})}\left\}\right\}$ $\left(7\right)$

$=\underset{b}{\arg \max }\left\{\mathrm{Re}\right\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\left[a_k^{*}\right|r,b^{\left(n\right)}\left]e^{-j2\mathrm{\&pi;\&Delta;fkT}}{e}^{-\mathrm{j\&Delta;\&theta;}}\right\}\}$

其中表示在给定r和b⁽ⁿ⁾的条件下

的条件期望；复数实部的最大值不能超过其模的最大值，且只有复数为实数时，实部的最大值和模的最大值等价，求公式（3）复数实部的最大值分两步：

A、求出使复数

模取最大值的Δf；

B、Δθ取复数 $\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\left[a_k^{*}\right|r,b^{\left(n\right)}]e^{-j2\mathrm{\&pi;\&Delta;fkT}}$ 的幅角，即旋转 $\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\left[a_k^{*}\right|r,b^{\left(n\right)}]e^{-j2\mathrm{\&pi;\&Delta;fkT}}$ 使其与实轴重合，此时模的最大值也就是实部的最大值；

根据上述分析，EM迭代算法中频偏及相偏的估计公式表示为

$\mathrm{\&Delta;}\hat{f}=\underset{\mathrm{\&Delta;f}}{\max }|\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a[a_k^{*}\left|{r,b}^{\left(n\right)}\right]e^{-j2\mathrm{\&pi;\&Delta;fkT}}|$

$\mathrm{\&Delta;}\widehat{\mathrm{\&theta;}}=\mathrm{angle}\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a{\left[a_k^{*}\right|{r,b}^{\left(n\right)}]e}^{-j2\mathrm{\&pi;\&Delta;}\hat{f}\mathrm{kT}}\}$

其中angle{·}表示取复数的幅角，

和

分别表示第n次EM迭代时Δf和Δθ的估计值，即对于BPSK调制，a中各信息比特取1或-1，即a_k＝1或a_k＝-1，因此

$E_a\left[a_k^{*}\right|r,b^{\left(n\right)}]=p(a_k=1|r,b^{\left(n\right)})-p(a_k=-1|r,b^{\left(n\right)})---\left(8\right)$

p(a_k＝1|r,b⁽ⁿ⁾)和p(a_k＝-1|r,b⁽ⁿ⁾)分别表示已知r和b⁽ⁿ⁾的条件下a_k取1和-1的概率，所以估计Δf和Δθ需知道a中信息比特a_k的联合统计p(a_k＝1|r,b⁽ⁿ⁾)和p(a_k＝-1|r,b⁽ⁿ⁾)，当迭代次数n足够大时，第n次迭代中产生的判决信息L(q_k)和a_k的联合统计满足如下关系

$L\left(q_k\right)\text{=ln}\frac{p(a_k=-1|r,b^{\left(n\right)})}{p(a_k=1|r,b^{\left(n\right)})}$

根据上式及p(a_k＝-1|r,b⁽ⁿ⁾)+p(a_k＝1|r,b⁽ⁿ⁾)=1得到联合统计

$p(a_k=1|r,b^{\left(n\right)})=\frac{e^{L\left(q_k\right)}}{1+e^{L\left(q_k\right)}},$ $p(a_k=-1|r,b^{\left(n\right)})=\frac{1}{1+e^{L\left(q_k\right)}}$

代入公式(4)得到

$E_a\left[a_k^{*}\right|{r,b}^{\left(n\right)}]=\frac{e^{L\left(q_k\right)}-1}{e^{L\left(q_k\right)}+1}$

上式就把LDPC迭代译码中获得的判决信息L(q_k)应用到了同步参数频偏和相偏的估计；上述同步参数频偏和相偏的估计分两步进行，即：

A、估计频偏

$\mathrm{\&Delta;}\hat{f}=\underset{\mathrm{\&Delta;f}}{\max }|\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k\&CenterDot;\frac{e^{L\left(q_k\right)}-1}{e^{L\left(q_k\right)}+1}\&CenterDot;e^{-j2\mathrm{\&pi;\&Delta;fkT}}|$

B、估计相偏

$\mathrm{\&Delta;}\widehat{\mathrm{\&theta;}}=\mathrm{angle}\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k\&CenterDot;\frac{e^{L\left(q_k\right)}-1}{e^{L\left(q_k\right)}+1}\&CenterDot;e^{-j2\mathrm{\&pi;\&Delta;}\hat{f}\mathrm{kT}}\}$

6）修正接收信息序列r

用当前频偏和相偏估计值对接收信息序列r进行修正；用r′表示修正后的接收信息序列，其第k个分量用r′_k表示，r_k表示r中第k个分量；修正公式为

$r_k^{\&prime;}=r_k{e}^{-j(\mathrm{\&Delta;}\widehat{\mathrm{\&theta;}}+2\mathrm{\&pi;\&Delta;}\hat{f}\mathrm{kT})}$

7）计算变量节点初始对数似然比消息

用修正后的接收信息序列r′计算各变量节点的初始对数似然比消息；在信道噪声方差为σ²的情况下，

所以第k个变量节点的初始对数似然比消息为

$L\left(P_k\right)=\ln \frac{p_k\left(0\right)}{p_k\left(1\right)}={|r_k^{\&prime;}-1|}^2-{|r_k^{\&prime;}+1|}^2$

8）变量节点消息更新

变量节点收集与它相邻的校验节点及信道的消息，

$L\left(q_{\mathrm{kj}}\right)=L\left(P_k\right)+\underset{j^{\&prime;}\&Element;C_k\backslash j}{\mathrm{\&Sigma;}}L\left(r_{j^{\&prime;}k}\right)$

其中L(q_kj)表示第k个变量节点传给第j个校验节点的对数似然比消息，L(P_k)是来自信道的初始对数似然比消息，L(r_j'k)是第j′个校验节点传给第k个变量节点的对数似然比消息，j'∈C_k\j表示j'是集合C_k\j的元素；

9）返回到步骤2)更新校验节点消息，进行下一次迭代。

Claims (1)

1.一种GPS接收机中LDPC码辅助的迭代载波同步方法，应用于GPS接收机中的LDPC码译码器，以实现精确估计同步不理想引起的频偏和相偏，提高整个系统性能的目的；预先设C是由m×n维校验矩阵H＝{h_j,k}确定的长为n的二元LDPC码；R_j表示同第j个校验节点相连的变量节点的集合，即R_j＝{k|h_k,j＝1}，R_j\k表示从R_j中除去第k个变量节点的集合；C_k表示同第k个变量节点相连的校验节点的集合，即C_k＝{j|h_k,j＝1}，C_k\j表示从C_k中除去第j个校验节点的集合；a表示码C中的码字经过BPSK调制后的发送信息序列，r表示a经过信道后的接收信息序列，也即LDPC译码器的输入，r中第k个分量对应码C的第k个变量节点；p_k(0)和p_k(1)分别表示第k个变量节点为0和1的概率，L(P_k)表示第k个变量节点的初始对数似然比消息；该方法的步骤如下：

1）初始化

假设BPSK将0映射到-1、1映射到1，则发送序列经过均值为0方差为σ²白高斯噪声信道后

r_k是r的第k个分量，所以第k个变量节点的初始对数似然比消息为

$L\left(P_k\right)=\ln \frac{p_k\left(0\right)}{p_k\left(1\right)}={|r_k-1|}^2-{|r_k+1|}^2$

2）校验节点消息更新

校验节点收集与它相邻的变量节点的消息；LDPC码译码器采用偏移最小和算法；当LDPC码译码器使用最小和译码时，令L'(r_jk)表示第j个校验节点传给第k个变量节点的对数似然比消息；最小和译码算法中校验节点消息更新的公式如下：

$L^{\text{\&prime;}}\left(r_{\mathrm{jk}}\right)=\underset{k^{\&prime;}\&Element;R_k\backslash j}{\mathrm{\&Pi;}}\mathrm{sign}\left(\text{L}\left(q_{k^{\&prime;}j}\right)\right)\&CenterDot;\underset{k^{\&prime;}\&Element;R_k\backslash j}{\min }\left(\right|L\left(q_{k^{\&prime;}j}\right)\left|\right)$

其中 $\mathrm{sign}\left(x\right)=\left\{\begin{array}{cc}1,& x\&GreaterEqual;0\\ 0,& x<0\end{array}\right.$ 是符号函数，符号min(·)表示求最小值，L(q_k′j)表示第k′个变量节点传给第j个校验节点的对数似然比消息，符号k'∈R_k\j表示k'是集合R_k\j的元素；在最小和算法中引进偏移因子β，便得到偏移最小和算法，其中β的取值介于0和1之间，用L(r_jk)表示偏移最小和算法中第j个校验节点传给第k个变量节点的对数似然比消息；偏移最小和算法中校验节点消息更新的公式如下：

L(r_jk)＝sign(L'(r_jk))·max(|L'(r_jk)|-β，0)

其中符号max(·)表示求最大值；不同LDPC码对应的偏移因子取值不同，GPS系统使用了两种LDPC码，码长分别为1200和578，这两种码的偏移因子分别通过仿真来求得，即让偏移因子在0和1间取不同值，根据仿真结果选择使误码率最低的那个偏移因子；

3）计算判决消息并做判决

变量节点收集与它相邻的校验节点的消息及来自信道的初始消息；用L(q_k)表示第k个变量节点的判决消息，则

$L\left(q_k\right)=L\left(P_k\right)+\underset{j^{\&prime;}\text{\&Element;}{C}_k}{\mathrm{\&Sigma;}}L\left(r_{j^{\&prime;}k}\right)$

如果L(q_k)＞0，则判定

否则判定

表示码字中第k个变量节点的判决输出；

4）判断译码是否结束

满足以下条件之一即表示译码结束：

A．

H是LDPC码的校验矩阵，表示

的转置；

B．已达到指定的最大迭代次数，

否则转入下一步；

5）若译码没有结束，则对同步参数进行一次估计；

用b＝[Δf,Δθ]表示同步不理想引起的频偏和相偏，其中Δf和Δθ分别表示频偏和相偏，使用最大对数似然概率的方法估计b，并用期望-最大（EM）算法求解该最大对数似然概率估计问题，EM迭代求解中要用到LDPC译码中产生的软信息；

b的最大对数似然概率估计是指找到使对数似然概率密度函数lnp(r|b)最大的b作为其估计值，记作

即

符号

表示求使括号内参数取最大值的b，EM算法通过引进先验分布已知的隐藏变量来达到简化计算的目的；取发送信息序列a作为先验分布已知的随机变量，则标准EM算法求解最大对数似然估计的迭代计算公式表示为

$b^{(n+1)}=\underset{b}{\arg \max }\left\{\underset{a}{\mathrm{\&Sigma;}}p\left(a\right|r,b^{\left(n\right)})\ln p\left(r\right|a,b)\right\}---\left(1\right)$

其中b⁽ⁿ⁾表示第n次迭代中的b估计值，b⁽ⁿ⁺¹⁾表示在b的当前估计值为b⁽ⁿ⁾时下一次的估计值，p(r|a,b)表示已知a和b的条件下r的条件概率密度函数，p(a|r,b⁽ⁿ⁾)表示已知a和b⁽ⁿ⁾的条件下r的条件概率密度函数，

是求和符号；对于GPS系统，在扩频序列和方波副载波同步理想，即接收端可以精确恢复它们的定时信息的情况下，整个系统可简化为BPSK调制系统，在已知发送信息序列a和同步参数b＝[Δf,Δθ]的条件下，译码器输入r的对数条件概率表示为

$\ln p\left(r\right|a,b)=-\frac{1}{{2\mathrm{\&sigma;}}^2}\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{|a_k{e}^{j(2\mathrm{\&pi;\&Delta;fkT}+\mathrm{\&Delta;\&theta;})}-r_k|}^2$ $\left(2\right)$

$=-\frac{1}{2{\mathrm{\&sigma;}}^2}\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}({\left|a_k\right|}^2+{\left|r_k\right|}^2)+\frac{1}{{\mathrm{\&sigma;}}^2}\mathrm{Re}\left\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{a}_k^{*}{e}^{-j(2\mathrm{\&pi;\&Delta;fkT}+\mathrm{\&Delta;\&theta;})}\right\}$

其中L是码字的长度，T是发送符号的持续时间，a_k和r_k分别是发送和接收信息序列中第k个分量的取值，σ²是噪声方差，Re{·}表示取复数的实部，

表示a_k的复共轭；公式（2）分成两部分，前一部分与估计参数b＝[Δf,Δθ]无关，它只影响公式（1）所能取的最大值，而不能影响使（1）式取得最大值的参数b，所以这一部分忽略；同理，公式（2）后一部分的比例因子

也去掉，所以对（1）式求解简化为

$b^{(n+1)}=\underset{b}{\arg \max }\left\{\underset{a}{\mathrm{\&Sigma;}}p\left(a\right|r,b^{\left(n\right)})\mathrm{Re}\right\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{a}_k^{*}{e}^{-j(2\mathrm{\&pi;\&Delta;fkT}+\mathrm{\&Delta;\&theta;})}\left\}\right\}$ $\left(3\right)$

$=\underset{b}{\arg \max }\left\{\mathrm{Re}\right\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\left[a_k^{*}\right|r,b^{\left(n\right)}\left]e^{-j2\mathrm{\&pi;\&Delta;fkT}}{e}^{-\mathrm{j\&Delta;\&theta;}}\right\}\}$

其中

表示在给定r和b⁽ⁿ⁾的条件下

的条件期望；复数实部的最大值不能超过其模的最大值，且只有复数为实数时，实部的最大值和模的最大值等价，求公式（3）复数实部的最大值分两步：

A、求出使复数

模取最大值的Δf；

B、Δθ取复数 $\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\left[a_k^{*}\right|r,b^{\left(n\right)}]e^{-j2\mathrm{\&pi;\&Delta;fkT}}$ 的幅角，即旋转 $\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\left[a_k^{*}\right|r,b^{\left(n\right)}]e^{-j2\mathrm{\&pi;\&Delta;fkT}}$ 使其与实轴重合，此时模的最大值也就是实部的最大值；

根据上述分析，EM迭代算法中频偏及相偏的估计公式表示为

$\text{\&Delta;}\hat{f}=\underset{\mathrm{\&Delta;f}}{\max }\left|\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\right[a_k^{*}|r,b^{\left(n\right)}]e^{-j2\mathrm{\&pi;\&Delta;fkT}}|$

$\text{\&Delta;}\widehat{\mathrm{\&theta;}}=\mathrm{angle}\left\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k{E}_a\right[a_k^{*}|r,b^{\left(n\right)}]e^{-j2\mathrm{\&pi;\&Delta;}\hat{f}\mathrm{kT}}\}$

其中angle{·}表示取复数的幅角，

和

分别表示第n次EM迭代时Δf和Δθ的估计值，即

对于BPSK调制，a中各信息比特取1或-1，即a_k＝1或a_k＝-1，因此

$E_a\left[a_k^{*}\right|r,b^{\left(n\right)}]=p(a_k=1|r,b^{\left(n\right)})-p(a_k=-1|r,b^{\left(n\right)})---\left(4\right)$

p(a_k＝1|r,b⁽ⁿ⁾)和p(a_k＝-1|r,b⁽ⁿ⁾)分别表示已知r和b⁽ⁿ⁾的条件下a_k取1和-1的概率，所以估计Δf和Δθ需知道a中信息比特a_k的联合统计p(a_k＝1|r,b⁽ⁿ⁾)和p(a_k＝-1|r,b⁽ⁿ⁾)，当迭代次数n足够大时，第n次迭代中产生的判决信息L(q_k)和a_k的联合统计满足如下关系

$L\left(q_k\right)\text{=ln}\frac{p(a_k=-1|r,b^{\left(n\right)})}{p(a_k=1|r,b^{\left(n\right)})}$

根据上式及p(a_k＝-1|r,b⁽ⁿ⁾)+p(a_k＝1|r,b⁽ⁿ⁾)=1得到联合统计

$p(a_k=1|r,b^{\left(n\right)})=\frac{e^{L\left(q_k\right)}}{1+e^{L\left(q_k\right)}},$ $p(a_k=-1|r,b^{\left(n\right)})=\frac{1}{1+e^{L\left(q_k\right)}}$

代入公式(4)得到

$E_a\left[a_k^{*}\right|{r,b}^{\left(n\right)}]=\frac{e^{L\left(q_k\right)}-1}{e^{L\left(q_k\right)}+1}$

上式就把LDPC迭代译码中获得的判决信息L(q_k)应用到了同步参数频偏和相偏的估计；上述同步参数频偏和相偏的估计分两步进行，即：

A、估计频偏

$\mathrm{\&Delta;}\hat{f}=\underset{\mathrm{\&Delta;f}}{\max }|\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k\&CenterDot;\frac{e^{L\left(q_k\right)}-1}{e^{L\left(q_k\right)}+1}\&CenterDot;e^{-j2\mathrm{\&pi;\&Delta;fkT}}|$

B、估计相偏

$\mathrm{\&Delta;}\widehat{\mathrm{\&theta;}}=\mathrm{angle}\{\underset{k=1}{\overset{L}{\mathrm{\&Sigma;}}}{r}_k\&CenterDot;\frac{e^{L\left(q_k\right)}-1}{e^{L\left(q_k\right)}+1}\&CenterDot;e^{-j2\mathrm{\&pi;\&Delta;}\hat{f}\mathrm{kT}}\}$

6）修正接收信息序列r

用当前频偏和相偏估计值对接收信息序列r进行修正；用r′表示修正后的接收信息序列，其第k个分量用r′_k表示，r_k表示r中第k个分量；修正公式为

$r_k^{\&prime;}=r_k{e}^{-j(\mathrm{\&Delta;}\widehat{\mathrm{\&theta;}}+2\mathrm{\&pi;\&Delta;}\hat{f}\mathrm{kT})}$

7）计算变量节点初始对数似然比消息

用修正后的接收信息序列r′计算各变量节点的初始对数似然比消息；在信道噪声方差为σ²的情况下，

所以第k个变量节点的初始对数似然比消息为

$L\left(P_k\right)=\ln \frac{p_k\left(0\right)}{p_k\left(1\right)}={|r_k^{\&prime;}-1|}^2-{|r_k^{\&prime;}+1|}^2$

8）变量节点消息更新

变量节点收集与它相邻的校验节点及信道的消息，

$L\left(q_{\mathrm{kj}}\right)=L\left(P_k\right)+\underset{j^{\&prime;}\&Element;C_k\backslash j}{\mathrm{\&Sigma;}}L\left(r_{j^{\&prime;}k}\right)$

其中L(q_kj)表示第k个变量节点传给第j个校验节点的对数似然比消息，L(P_k)是来自信道的初始对数似然比消息，L(r_j'k)是第j′个校验节点传给第k个变量节点的对数似然比消息，j'∈C_k\j表示j'是集合C_k\j的元素；

9）返回到步骤2)更新校验节点消息，进行下一次迭代。