逐幸存路径处理(PSP)算法

混合信号处理的发展历程大致可以分为四个阶段：
一、采用过采样的方法，达到类似多通道接收效果，借助已有的多通道处理方式进行分离
二、通过变换域方法实现混合信号分离
三、通过最大似然序列联合估计方法分离混合信号
四、基于深度神经网络的解调算法

混合信号盲分离的核心是从一路接收信号中估计出两路信号的符号序列。传统的处理思想是先实现混合信号的分离，再进行解调。该类方法的核心思想是利用分量信号某个参数的差异构造正定条件，寻求在某个变换域上“正交”，或采用多通道分离方法实现分离，然而该类方法的性能往往受多种参数的影响，且分离后还需要解调，存在错误传播现象，不能满足实际需求。

鉴于此，跳过分离环节，直接对接收混合信号进行符号序列的盲检测，将信号调制参数的影响等效成一信道滤波器，恢复出通信双方发送的符号序列，从而避免了先分离后解调类处理方法带来的不足，这就是联合信道估计与序列检测的思想，可以用粒子滤波、逐幸存路径算法进行实现。

逐幸存路径处理(Per-Survivor Processing, PSP)

1.信号模型及PSP算法描述

对于一个接收机接收到的一个包含两个QPSK信号的混合信号，等效到基带后可以表示为
$$\begin{array}{r}y(t)=h_1{e}^{j(\Delta {\omega }_{1}t+{\theta }_1)}{x}_1(t)+h_2{e}^{j(\Delta {\omega }_{2}t+{\theta }_2)}{x}_2(t)+v(t)\end{array}\text{\hspace{1em}(1)}$$其中 $h_i$ 是两个调制信号的实际振幅衰减，$\Delta {\omega }_i$ 是两个信号的频率偏移，${\theta }_i$ 是初始相位，$v(t)$ 是加性高斯白噪声。中间信号 $x_i(t)$ 可以定义为
$$\begin{array}{rl}{x}_i(t)=& \sum _{n=-\infty }^{\infty }{s}_n^{(i)}{g}_i(t-nT+{\tau }_i(t))\end{array}\text{\hspace{1em}(2)}$$其中，其中 $s_n^{(i)}$ 是具有单位能量（QPSK调制）、各符号出现概率等概、独立同分布的随机序列。$g_i(t)$ 是由成型滤波器、信道滤波器和匹配滤波器组成的等效信道滤波器的脉冲响应，$0\le {\tau }_i(t)<t$是两个调制信号的相对时延，$T$ 是符号周期。

以 $1/T$ 的速率采样公式(1)中的信号，假设等效信道滤波器的脉冲响应持续有限时间 $2LT$，从 $(1-L)T$ 到 $LT$，可以得到采样后的序列$$\begin{array}{r}{y}_k=h_1{e}^{j(\Delta {\omega }_{1}kT+{\theta }_1)}{x}_k^{(1)}+h_2{e}^{j(\Delta {\omega }_{2}kT+{\theta }_2)}{x}_k^{(2)}+v_k\end{array}\text{\hspace{1em}(3)}$$其中$y_k=y(kT),v_k=v(kT)$，且
$$\begin{array}{r}{x}_k^{(i)}=x_i(kT)=\sum _{n=1-L}^L{s}_{k+n}^{(i)}{g}_i(-nT+{\tau }_k^{(i)})\end{array}\text{\hspace{1em}(4)}$$其中 ${\tau }_k^{(i)}={\tau }_i(kT),k=0,1,...,K-1$，
通过定义 $2L\times 1$ 源符号向量${\text{s}}_k^{(i)}={\left[s_{k-L+1}^{(i)},s_{k-L+2}^{(i)},\dots ,s_{k+L}^{(i)}\right]}^{\mathrm{T}}$ 以及对应的CIR(channel impulse responses, 广义的信道脉冲响应)向量$\begin{array}{c}{\text{f}}_k^{(i)}=h_i{e}^{j(\Delta {\omega }_{i}kT+{\theta }_i)}.\end{array}{\left[g_i((L-1)T+{\tau }_k^{(i)}),\dots ,g_i(-LT+{\tau }_k^{(i)})\right]}^T$。通过 ${\text{s}}_k^{(1)},{\text{s}}_k^{(2)}$ 和 ${\text{f}}_k^{(1)},{\text{f}}_k^{(2)}$ ，可以构造向量${\text{s}}_k=\left[\begin{array}{c}{\text{s}}_k^{(1)}\\ {\text{s}}_k^{(2)}\end{array}\right]$ 和 ${\text{f}}_k=\left[\begin{array}{c}{\text{f}}_k^{(1)}\\ {\text{f}}_k^{(2)}\end{array}\right]$，并将公式(3)改写成向量形式$$y_k={\text{f}}_k^T{\text{s}}_k+v_k\text{\hspace{1em}(5)}$$
上述是关于接收信号的描述，对于混合信号的接收，我们的最终目的是恢复出发送的信息序列，即得到发送的符号序列 $\begin{array}{rl}{s}_{0:K-1}^{(i)}& =\end{array}\left\{s_0^{(i)},s_1^{(i)},\dots ,s_{K-1}^{(i)}\right\}$，可用信息是接收到的信息序列 $y_{0\text{:}K-1}=\{y_0,y_1,\dots ,y_{K-1}\}$。如果CIR $f_k$ 已知，那么可以用序列检测，如使用 Viterbi 算法。显然，对于盲信号分离，$f_k$ 是未知的，需要将 $f_k$ 和符号序列 $s_{0\text{:}K-1}^{(i)}$一起进行联合估计。

通过实验发现，当载波频率之间的位移相对较大时，$f_k$ 随时间 $k$ 的变化太快，现有的算法如 最小均方算法(LMS) 和递归最小二乘(RLS)无法跟踪。在这种情况下，必须修改信号模型以保持最终形式的CIR尽可能恒定。

假设载波频偏被估计为$\widehat{\Delta {\omega }_i}$，有一个小误差$\widetilde{\Delta {\omega }_i}$，即
$$\Delta {\omega }_i=\widehat{\Delta {\omega }_i}+\widetilde{\Delta {\omega }_i}\text{\hspace{1em}(6)}$$将公式(6)代入到公式(3)，将频偏估计 $\widehat{\Delta {\omega }_i}$ 与源符号 ${\text{s}}_k^{(i)}$ 结合，在CIR中只留下误差 $\widetilde{\Delta {\omega }_i}$，可以得到修正后的模型为
$$y_k={\underline{\text{f}}}_k^T{\underline{\text{s}}}_k+v_k\text{\hspace{1em}(7)}$$其中 ${\underline{\text{s}}}_k=\left[\begin{array}{c}{\underline{\text{s}}}_k^{(1)}\\ {\underline{\text{s}}}_k^{(2)}\end{array}\right]$，${\underline{\text{f}}}_k=\left[\begin{array}{c}{\underline{\text{f}}}_k^{(1)}\\ {\underline{\text{f}}}_k^{(2)}\end{array}\right]$，且
$${\underline{\text{s}}}_k^{(i)}=e^{j\widehat{\Delta {\omega }_i}kT}\cdot {\left[s_{k-L+1}^{(i)},\dots ,s_{k+L}^{(i)}\right]}^T\text{\hspace{1em}(8)}$$
$$\begin{array}{rl}{\underline{\text{f}}}_k^{(i)}& =h_i{e}^{j(\widetilde{\Delta {\omega }_i}kT+{\theta }_i)}.\\ & {\left[g_i((L-1)T+{\tau }_k^{(i)}),\dots ,g_i(-LT+{\tau }_k^{(i)})\right]}^T\end{array}\text{\hspace{1em}(9)}$$现在可以观察到，新的CIR ${\underline{\text{f}}}_k$不再随时间k快速变化，可以用一个简单的LMS算法来处理。

定义 ${\varphi }_k=\left({\underline{s}}_k^{(1)},{\underline{s}}_k^{(2)}\right)$，定义在 $k$ 时刻信道状态为 μ k = ϕ k − L + 2 : k + L = { ϕ k − L + 2 , … , ϕ k + L } \mu_k=\phi_{k-L+2:k+L}=\{\phi_{k-L+2},\ldots,\phi_{k+L}\} μk​=ϕk−L+2:k+L​={ϕk−L+2​,…,ϕk+L​}，当 ${\varphi }_{k\text{+}L}$ 可得到时，状态从 ${\mu }_{k-1}$ 转移到 ${\mu }_k$，输出 $y_k$。这种转移可以写成
$${\mu }_{k-1}\stackrel{({\varphi }_{k+L},y_k)}{⟶}{\mu }_k$$ 对于QPSK调制，符号有四种可能。当成型滤波采用升余弦脉冲时，$g_i(t)$ 随时间 $t$ 衰减非常快，当 L 取 2 时信道可以描述为 $4^{2(2\times 2-1)}=4096$ 个状态的网格。为了进一步降低复杂性，引入DFSE(decision feedback sequence estimation, 判决反馈序列估计)检测源符号。为了保持分离性能，我们将最不显著的抽头作为唯一的反馈抽头。因此，我们将 $k$ 时刻的状态重新定义为 ${\mu }_k={\varphi }_{k+1:k+2}$。现在状态数减少到256，方便更快速地进行符号检测。

PSP算法实施如下：
1.初始化：令 $k=0$，对于每个状态 ${\underline{\mu }}_0$，用某个特定值初始化 ${\underline{\text{f}}}_0$，保存幸存路径 $p_0=({\underline{\mu }}_0)$ 及其对应的度量 $\Gamma ({\underline{\mu }}_0)=0$。
2.符号检测(DFSE)：对第 $k$ 步的分支，度量评估所有可能的状态转换，根据$$\begin{array}{r}\lambda ({\underline{\mu }}_k\to {\underline{\mu }}_{k+1})={\left|e({\underline{\mu }}_k\to {\underline{\mu }}_{k+1})\right|}^2\end{array}\text{\hspace{1em}(10)}$$其中
$$\begin{array}{r}e({\underline{\mu }}_k\to {\underline{\mu }}_{k+1})=y_k-{\underline{\text{f}}}_k^T{\underline{\text{s}}}_k({\underline{\mu }}_k\to {\underline{\mu }}_{k+1})\end{array}\text{\hspace{1em}(11)}$$其中 ${\underline{\text{s}}}_k({\underline{\mu }}_k\to {\underline{\mu }}_{k+1})$ 是相对于 ${\underline{\mu }}_k\to {\underline{\mu }}_{k+1}$ 状态转移的符号向量。对于所有后续状态 $\underline{\mu }{}_{k+1}$，累积的度量是通过对当前状态 ${\underline{\mu }}_k$ 进行最小化来确定的。
$$\begin{array}{r}\Gamma ({\underline{\mu }}_{k+1})=\underset{k}{\min }\left(\Gamma ({\underline{\mu }}_k)+\lambda ({\underline{\mu }}_k\to {\underline{\mu }}_{k+1})\right)\end{array}\text{\hspace{1em}(12)}$$最后，根据 $p_{k+1}=(p_k,{\underline{\mu }}_{k+1})$，通过合并满足（12）的转移路径来扩展终止于当前状态的幸存节点。
3.参数估计：对于转移路径 ${\underline{\mu }}_k\to {\underline{\mu }}_{k+1}$ 来扩展当前的幸存节点，CIR向量根据下式更新
$${\underline{\text{f}}}_{k+1}={\underline{\text{f}}}_k+\gamma \cdot e({\underline{\mu }}_k\to {\underline{\mu }}_{k+1})\cdot {\underline{\text{s}}}_k^{\ast }({\underline{\mu }}_k\to {\underline{\mu }}_{k+1})\text{\hspace{1em}(13)}$$其中 $\ast$ 表示共轭，$\gamma$ 是控制调整速率的尺度因子。
4.令 $k=k+1$，当 $k=K$时结束算法，否则回到步骤2.

与部分算法需要估计参数 $h_i,\Delta {\omega }_i,{\theta }_i,{\tau }_k^{(i)}$ 的方法不同，PSP算法仅估计CIR向量 ${\underline{\text{f}}}_k^{(i)}$

2.初始化及定时问题

通过对PSP算法的描述，可以看出其中各参数的初始化是非常重要的。好的初始化可以确保结果收敛并加速迭代过程，坏的初始化可能导致结果不收敛。为了获得 ${\underline{\text{f}}}_0$ 真实值的好的近似，根据公式(9)所述需要分别估计 ${\text{h}}_i,{\theta }_i,g_i(-nT+{\tau }_k^{(i)})$。

还有一个问题是定时，当进行连续分离时，定时是需要考虑的。在许多情况下，符号周期 $T$ 是未知的，需要估计，因此误差是不可避免的。当由于时间增加而累积的定时误差超过一个符号周期（例如T）时，主抽头将滑向区域 $[1-L，L]$ 的任一边缘，甚至滑出区域，这导致信号重构的不准确，进而导致分离性能的下降。此外，还需要在精确点对接收到的混合信号进行采样，以满足条件 ${\tau }_k^{(1)}+{\tau }_k^{(2)}=T$，从而实现最佳的信号分离[1]。

不失一般性，可以假定载波频率偏移只有符号不同，即 $\begin{array}{r}\Delta {\omega }_1=-\Delta {\omega }_2\end{array}$.

2.1 幅度估计

根据两个源序列内部和之间的独立性，忽略噪声的影响，将接收到的混合信号 $y_k$ 的平均功率可以从（3）写成
$$\begin{array}{rl}E(|y_k{|}^2)=& h_1^2\sum _{n=1-L}^L{g}_1^2(-nT+{\tau }_k^{(1)})+\\ & h_2^2\sum _{n=1-L}^L{g}_2^2(-nT+{\tau }_k^{(2)})\\ \approx & 2{\bar{h}}^2{P}_0\end{array}\text{\hspace{1em}(14)}$$其中 ${\bar{h}}^2=(h_1^2+h_2^2)/2,P_0=g_i^2(0)$。当两个混合信号有相近的功率时，根据公式(14)幅度 $h_i$ 可以被估计为$$\begin{array}{r}{\hat{h}}_i=\sqrt{\frac{1}{2K{P}_0}\sum _{k=0}^{K-1}|y_k{|}^2}\end{array}\text{\hspace{1em}(15)}$$

2.2 初始相位估计

构造辅助信号
$$\begin{array}{rl}{Z}_k^{(i)}& =y_k{e}^{-j\Delta {\omega }_{i}kT}\\ & =h_i{e}^{j{\theta }_i}{x}_k^{(i)}+h_l{e}^{j((\Delta {\omega }_l-\Delta {\omega }_i)kT+{\theta }_l)}{x}_k^{(l)}\end{array}\text{\hspace{1em}(16)}$$其中 i , l ∈ { 1 , 2 } , i ≠ l i,l\in\{1,2\},i\neq l i,l∈{1,2},i=l，可以得到其四次形式如下
$$E\{(Z_k^{(i)}{)}^4\}=h_i^{4}E\{(x_k^{(i)}{)}^4\}{e}^{j4{\theta }_i}\text{\hspace{1em}(17)}$$初始相位 ${\theta }_i$ 可以被估计为$$\widehat{{\theta }_i}=\frac{1}{4}\text{arg}\left[E\{(Z_k^{(i)}{)}^4\}\right]\text{\hspace{1em}(18)}$$这里存在一个 $k\pi /2,k=0,1,2,3$ 的模糊，但通常不太重要。

2.3 脉冲响应初始化和定时估计

通过 2.1 和 2.2可以估计出$h_i,{\theta }_i$，如果再获得时间延迟 ${\tau }_k(i)$，我们可以正确地初始化 $g_i(-nT+{\tau }_k^{(i)})$，从而得到CIR各参数。然而 ${\tau }_k^{(i)}$ 的估计是较困难的。其中一种简单的做法是将其设置为中间值。

对于定时，Gardner定时同步算法广泛地应用在单调制信号解调中，其可以通过以下方式定位不受码间串扰(ISI)影响的最佳内插时刻：
$$t_k=t_{k-1}+\hat{T}+\eta \hat{T}{e}_k\text{\hspace{1em}(19)}$$ 其中 $t_k=kT-{\xi }_k$ 表示第 $k$ 次插值的时刻，$\hat{T}$ 是符号周期的估计值，$\eta$是归一化的步长参数，Gardner定时误差检测器的输出值可以写成
$$\begin{array}{rl}{e}_k& =\Re \left\{\left[y(t_{k-1})-y(t_k)\right]\cdot y^{\ast }(\frac{t_{k-1}+t_k}{2})\right\}\end{array}\text{\hspace{1em}(20)}$$其中 $\Re$ 意味着取实部。

对于两个QPSK信号的单通道接收混合信号，将公式(1)，带入到公式(20)，并且令
$h_1=h_2=h,\Delta {\omega }_1=-\Delta {\omega }_2=\Delta \omega ,g_1(t)=g_2(t)=g(t),{\tau }_i(t_{k-1})={\tau }_i(t_k)={\tau }_i((t_{k-1}+t_k)/2)={\tau }_i,{\xi }_{k-1}={\xi }_k=\xi$ 得到
E ( e k ) = h 2 cos ( Δ ω T 2 ) ⋅ { ∑ n = − ∞ ∞ g ( − n T − T 2 − ξ + τ 1 ) ⋅ [ g ( − n T − T − ξ + τ 1 ) − g ( − n T − ξ + τ 1 ) ] + ∑ n = − ∞ ∞ g ( − n T − T 2 − ξ + τ 2 ) ⋅ [ g ( − n T − T − ξ + τ 2 ) − g ( − n T − ξ + τ 2 ) ] } (21) \begin{aligned}E(e_k)&=h^2\text{cos}(\frac{\Delta\omega T}{2})\cdot\\ &\{\sum_{n=-\infty}^\infty g(-nT-\frac{T}{2}-\xi+\tau_1)\cdot[g(-nT-T-\xi+\tau_1)-g(-nT-\xi+\tau_1)]\\ +&\sum_{n=-\infty}^\infty g(-nT-\frac{T}{2}-\xi+\tau_2)\cdot[g(-nT-T-\xi+\tau_2)-g(-nT-\xi+\tau_2)]\} \end{aligned}\tag{21} E(ek​)+​=h2cos(2ΔωT​)⋅{n=−∞∑∞​g(−nT−2T​−ξ+τ1​)⋅[g(−nT−T−ξ+τ1​)−g(−nT−ξ+τ1​)]n=−∞∑∞​g(−nT−2T​−ξ+τ2​)⋅[g(−nT−T−ξ+τ2​)−g(−nT−ξ+τ2​)]}​(21)可以将符号定时误差定义为
$$d=\xi -\frac{{\tau }_1+{\tau }_2}{2}\text{\hspace{1em}(22)}$$显然，利用升余弦的轴对称特性，$E(e_k)$是关于 $d$ 的奇函数。因此，可以认为，Gardner算法是适合用在混合信号上的，当以上假定全满足时，算法将收敛于 $({\tau }_1+{\tau }_2)/2$，即源信号最佳插值点的中点。显然，这些假设并不能全满足，当条件未完全满足时，例如当 $h_1\ne h_2,g_1(t)\ne g_2(t)$ 同时满足或满足其中一个时，仍可以期望Gardner算法收敛到源信号两个最佳内插时刻之间的某一点。

通过Gardner算法后，时延 ${\tau }_k^{(1)}$ 和 ${\tau }_k^{(2)}$ 将落于 $[-T/4,T/4)$ 且符号相反。为了满足最佳性能条件，在采样的时候可以延迟 $T/2$，然后可以使用区域 $[T/4,3T/4)$ 的中点 $T/2$ 来初始化等效信道滤波器的脉冲响应。

3.仿真结果

对于两路 QPSK 信号，采用升余弦成型，滚降系数为 0.33，两路幅度为 1，频偏设置为 $\pm \frac{1}{100}T$，可以得到结论，PSP 算法非常接近信道已知的 Viterbit 算法性能。

总结

PSP 算法就是 Viterbi + 信道反馈估计，实现了时变信道下最大似然序列检测，在信道截短长度足够长时，其性能逼近最大似然解，但其复杂度随着信号阶数和信道截短长度爆炸上升。

参考文献

[1]Tu S , Hui Z , Na G . Single-channel blind separation of Two QPSK signals using per-survivor processing[C]// IEEE Asia Pacific Conference on Circuits & Systems. IEEE, 2009.
[2]熊涛. 同频混合信号单通道盲分离关键技术研究[D].信息工程大学,2015.