高速移动场景包含车对车通信、高速铁路通信，是无线通信技术的重要场景之一. 传统的正交频分复用（orthogonal frequency division multiplexing, OFDM）一直是大多数无线通信系统的波形选择，然而高移动性引起的多普勒频移会导致OFDM调制技术性能下降. 为了解决上述问题，正交时频空（orthogonal time frequency space, OTFS）调制作为新的二维波形调制技术应运而生. 与OFDM仅将信号调制在频域的方式不同，OTFS将信息符号在时延-多普勒域复用^[1]. OTFS调制将时频域的时变信道转换为时延-多普勒域的准静态信道，发射符号在时延-多普勒域中经历近乎恒定的信道增益，因此与OFDM相比，OTFS是更具前景的调制技术^[2].

Hadani等^[3]研究了OTFS在毫米波频率下的性能评估，表明大多情况下OTFS的误码率性能都优于OFDM的. Raviteja等^[4]以简单的矩阵形式推导出实际脉冲成形波形的时延-多普勒域输入输出关系，然后将其推广到任意波形，并得出如下结论：由于只有少数反射器参与信号传播，信道的时延-多普勒域表示是稀疏的，这大大降低了接收机的复杂性.

类似于OFDM多天线系统，将多输入多输出（multiple-inputs multiple-outputs, MIMO）与OTFS技术结合^[5-9]可以进一步提高频谱效率. 广义空间调制（generalized spatial modulation, GSM）被认为是MIMO技术的一个特例，它既继承了MIMO的优点，又避免了MIMO的不足（如天线间同步、信道间干扰以及需要多个射频链），并且GSM可以利用发射天线的空间位置来携带额外的比特，进一步提高系统的频谱效率. 目前已经有部分学者对GSM-OTFS结构做出相关研究，主要关于性能研究^[10-12]和信号检测算法设计^[13-14] 2方面. Zou等^[13]提出基于最小均方误差（minimum mean square error, MMSE）准则的决策反馈检测器. Wang等^[14]提出基于广义近似消息传递（generalized approximate message passing，GAMP）的迭代检测器，仿真结果表明，与传统的MMSE检测器相比，所提出的GAMP检测器具有更好的性能.

在享受GSM优点同时，GSM结构给信号检测算法设计引入的难度也不容忽视，消息传递类算法经常被用于解决此类问题. Fan等^[15]针对GSM结构基于分层消息传递算法设计了2种检测器，分别为阈值辅助的分层消息传递检测器和概率排序辅助的分层消息传递检测器. 同样是GSM结构，Fan等^[16]针对近似消息传递算法设计了基于概率排序的近似消息传递检测器. Wei等^[17]研究了近似消息传递（approximate message passing, AMP）算法在多个激活资源（天线或子载波）的索引调制（index modulation, IM）信号检测中的应用. Zhang等^[18]提出在期望传播算法框架下的联合信号检测算法，用于具有广义空间调制的上行链路大规模多用户MIMO系统. 期望传播（expectation propagation, EP）和置信传播（belief propagation, BP）算法已经被广泛应用于迭代接收机^[19]，在此启发下本研究提出混合期望传播和置信传播（hybrid expectation propagation and belief propagation, EP-BP）的联合信号检测方案，用于GSM-OTFS系统的激活天线组合索引检测和星座符号解调.

区别于大部分GSM结构下信号传递（message passing, MP）类算法^[15-17]先将消息传递过程中的概率分布近似为一元复高斯分布，再增加GSM约束的做法，在所提出的EP-BP信号检测算法中，取自发射符号集合的向量被视作矢量变量节点（variable node, VN）.

为了降低和积算法执行消息传递的计算复杂度^[20]，将离散概率分布近似为多元复高斯概率分布，从而将复杂的和积更新规则转换为均值向量和协方差矩阵的迭代传递. 通过迭代更新发射符号矢量的符号置信度来确定激活天线组合索引和星座符号. 为了进一步降低符号置信度计算须遍历发射符号集合带来的计算复杂度，设计了一种分步的两阶段混合期望传播与置信传播算法TS-EP-BP. 在TS-EP-BP的第1阶段，采用EP-BP算法估计激活天线组合的索引值. 在得到估计索引值之后，发射符号向量中不同激活天线上的星座符号可以看作独立分布，符号置信度的计算次数减少到调制星座的基数.

GSM-OTFS的系统模型如图1所示，系统共有$M$个子载波和$N$个时隙，以及${N_{\text{t}}}$根发射天线、${N_{\text{r}}}$根接收天线. 在每个传输时隙，系统仅选择${N_{\text{a}}}$根发送天线进行激活，其中$1 \leqslant {N_{\text{a}}} \leqslant {N_{\text{t}}}$，其余天线在该时隙保持沉默. 图中，灰色填充的天线为激活状态，白色填充的天线为静默状态.

在GSM-OTFS传输中，输入的比特信息被划分为$MN$个比特块，每个比特块有Q位信息比特作为选择激活天线组合索引和${N_{\text{t}}}$个${M_{ \text{mod} }}$阶调制星座符号索引，其中，$ Q = {Q_{\text{a}}}+{N_{\text{t}}}{Q_{\text{m}}} = \left\lfloor {{{\log }_2}\left( {{\mathrm{C}}_{{N_{\text{t}}}}^{{N_{\text{a}}}}} \right)} \right\rfloor + {N_{\text{t}}} \times {\log _2}\;{M_{ \text{mod} }} $，$\left\lfloor \bullet \right\rfloor $表示向下取整，$ {Q_{\text{a}}} $和$ {Q_{\text{m}}} $分别表示激活天线组合索引和星座图符号的比特数. 星座图集合$ {{{C}}} $的基为C_m，$ {C_{\text{m}}} = \left| {{{C}}} \right| = {2^{{Q_{\text{m}}}}} $. 每个比特块的激活天线组合索引均为一个${N_{\text{t}}} \times 1$的二元向量，其中元素1和0代表天线激活和非激活2种状态，激活天线组合索引矩阵$ {{\boldsymbol{A}}} $共有${C_{\text{a}}} = {2^{{Q_{\text{a}}}}}$个列向量，表示${C_{\text{a}}}$种激活方式. 为了最小化干扰，须对激活天线组合进行精巧设计使得各天线被等概激活. 以${N_{\text{t}}} = 4$、${N_{\text{a}}} = 2$为例，激活天线组合图样矩阵可以设计为

(1)$ {{\boldsymbol{A}}} = \left[ {\begin{array}{*{20}{c}} 1&1&0&0 \\ 1&0&1&0 \\ 0&1&0&1 \\ 0&0&1&1 \end{array}} \right]. $

激活天线组合索引矩阵在收发端均为已知信息，结合激活天线组合图样和调制星座图集合${{{C}}}$,发射符号集合用${{{S}}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{{C}}}$表示，其中的列向量可以表示为${{\boldsymbol{\beta}} } \in {{{S}}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{{C}}}$，且共有$ {{{C}}_{\text{c}}} = \left| {{{{S}}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{{C}}}} \right| = {C_{\text{a}}} \times {\left( {{C_{\text{m}}}} \right)^{{N_{\text{a}}}}} $个列向量. 以BPSK为例，${{{S}}}_{4,2}^{{{{{C}}}_{{\text{BPSK}}}}}$表达式如下：

(2)$ \begin{split} S_{4,2}^{{C_{{\rm{BPSK}}}}} = & \left\{ { \left[ {\begin{array}{*{20}{c}}{ - 1}\\{ - 1}\\0\\0\end{array}} \right],\left[ {\begin{array}{*{20}{c}}{ - 1}\\{ 1}\\0\\0\end{array}} \right],\left[ {\begin{array}{*{20}{c}}{ 1}\\{ - 1}\\0\\0\end{array}} \right],\left[ {\begin{array}{*{20}{c}}{ 1}\\{ 1}\\0\\0\end{array}} \right],} \right.\\&\quad \left[ {\begin{array}{*{20}{c}}{ - 1}\\0\\{ - 1}\\0\end{array}} \right],\left[ {\begin{array}{*{20}{c}}{ - 1}\\0\\{ 1}\\0\end{array}} \right],\left[ {\begin{array}{*{20}{c}}{ 1}\\0\\{ - 1}\\0\end{array}} \right],\left[ {\begin{array}{*{20}{c}}{ 1}\\0\\{ 1}\\0\end{array}} \right],\\&\quad \left[ {\begin{array}{*{20}{c}}0\\{ - 1}\\0\\{ - 1}\end{array}} \right],\left[ {\begin{array}{*{20}{c}}0\\{ - 1}\\0\\{ 1}\end{array}} \right],\left[ {\begin{array}{*{20}{c}}0\\{ 1}\\0\\{ - 1}\end{array}} \right],\left[ {\begin{array}{*{20}{c}}0\\{ 1}\\0\\{ 1}\end{array}} \right],\\&\quad \left. {\left[ {\begin{array}{*{20}{c}}0\\0\\{ - 1}\\{ - 1}\end{array}} \right],\left[ {\begin{array}{*{20}{c}}0\\0\\{ - 1}\\{ 1}\end{array}} \right],\left[ {\begin{array}{*{20}{c}}0\\0\\{ 1}\\{ - 1}\end{array}} \right],\left[ {\begin{array}{*{20}{c}}0\\0\\{ 1}\\{ 1}\end{array}} \right]} \right\}.\end{split} $

综上所述，第$k$个时延第$l$个多普勒的GSM符号可以表示为

(3)$ {{\boldsymbol{x}}}(k,l) = {[ \cdots ,0,{s_1},0, \cdots ,0,{s_2},0, \cdots ,0,{s_{{N_{\text{a}}}}},0, \cdots ]^{\text{T}}}. $

一个GSM符号中共有$ {N_{\text{a}}} $个非零符号，非零符号的数量和排列位置都与激活天线组合相对应，$k \in \left\{ {1,2, \cdots ,N} \right\}$，$l \in \left\{ {1,2, \cdots ,M} \right\}$. 为了得到每根发射天线上的信号，重新定义GSM符号为$ {{\boldsymbol{x}}}(k,l) = {[{x_1}(k,l),{x_2}(k,l), \cdots ,{x_{{N_{\text{t}}}}}(k,l)]^{\text{T}}} $.

在GSM-OTFS系统的发射端，对发射符号${{\boldsymbol{x}}_e}$进行逆辛有限傅里叶变换（inverse symplectic finite Fourier transform, ISFFT）操作，其中${{\boldsymbol{x}}_e}$表示第$e$根发射天线上的所有信号. 具体而言，${{\boldsymbol{x}}_e}$的第$k$个时延第$l$个多普勒上的ISFFT操作可以表示为

(4)$ \begin{split} {{{\boldsymbol{X}}}_e}[n,m] =& {\text{ISFFT}}\;({{{\boldsymbol{x}}}_e}[k,l])= \frac{1}{{\sqrt {MN} }} \times \\ &\sum\limits_{k = 0}^{N - 1} {\sum\limits_{l = 0}^{M - 1}} {{{{\boldsymbol{x}}}_e}[k,l]\exp \left({\rm{j}}2{\text{π}} \left(\frac{{nk}}{N} - \frac{{ml}}{M}\right)\right)} . \end{split} $

随后，对第$e$根发射天线上，第$n$个时隙和第$m$个子载波的时频域发射符号$ {{{\boldsymbol{X}}}_e}[n,m] $进行海森堡变换，转为第$e$根天线上的时域符号：

(5)$ {{\boldsymbol{s}}_e}(t) = \sum\limits_{n = 0}^{N - 1} {\sum\limits_{m = 0}^{M - 1} {{{{\boldsymbol{X}}}_e}[n,m]{g_{{\text{tx}}}}\;(t - n{T_{\text{s}}}){\exp\;\left[{{\rm{j}}2{\text{π}} m\Delta f(t - n{T_{\text{s}}})}\right]}} } . $

式中：${g}_{\text{tx}}\;(\cdot )$为发送脉冲，$ {T_{\text{s}}} $表示沿时间$t$的采样间隔，$\Delta f$表示频域子载波间隔. 时域信号经过信道到达接收端的过程可以由如下公式表示：

(6)$ {{\boldsymbol{r}} = {\boldsymbol{Hs}}+{\boldsymbol{n}}}{.} $

式中：r(t)表示接收信号矢量，$ {\boldsymbol{r}}(t) = [ {{\boldsymbol{r}}_1}(t),{{\boldsymbol{r}}_2}(t), \cdots , {{\boldsymbol{r}}_g}(t), \cdots ,{{\boldsymbol{r}}_{{N_{\text{r}}}}}(t) ]^{\text{T}} \in {{\mathbf{C}}^{MN{N_{\text{r}}} \times 1}} $，$ {{\boldsymbol{r}}_g}(t) \in {{\mathbf{C}}^{MN \times 1}} $，按照$N$个时隙的顺序排列，每个时隙包含$M$个符号. H为时域等效信道矩阵，$ {{\boldsymbol{H}}} \in {{\mathbf{C}}^{MN{N_{\text{r}}} \times MN{N_{\text{t}}}}} $，构成元素${{{\boldsymbol{H}}}_{e,g}} = \sum\limits_{p = 1}^P {{h_p}{\boldsymbol{\varPi}} _{MN}^{{l_p}}{\boldsymbol{\varDelta}} _{MN}^{{k_p}+{K_p}}} \in {{\mathbf{C}}^{MN \times MN}}$^[4]，表示第$e$根发射天线与第$g$根接收天线之间的信道矩阵，其中${{\boldsymbol{\varPi}} _{MN}} \in {{\mathbf{C}}^{MN \times MN}}$表示第1列为$\left[0,\;\;1,\;\; \cdots ,\;\;0\right]^{\rm{T}}_{MN}$的前向循环移位矩阵，${{\boldsymbol{\varDelta }}_{MN}} = {\text{diag}}\;[ 1,{{\rm{e}}^{2{\text{π}} {\rm{j}}/(MN)}}, \cdots , {{\rm{e}}^{2{\text{π}} {\rm{j}}(MN - 1)/(MN)}} ] \in {{\mathbf{C}}^{MN \times MN}}$表示对角矩阵，$P$为传播路径数，h_p、τ_p、v_p分别为第$p$条路径的路径增益、时延与多普勒分量，$ {\tau _p} = {l_p}/(M\Delta f) $，$ {v_p} = ( {k_p}+ {K_p} )/(N{T_{\mathrm{s}}}), $$ {l_p} $和$ {k_p} $分别表示对应信道的整数时延抽头和多普勒抽头，$ {K_p} $为分数多普勒抽头. s(t)表示发射信号矢量，$ {\boldsymbol{s}}(t) = \left[ {{\boldsymbol{s}}_1}(t),{{\boldsymbol{s}}_2}(t), \cdots ,{{\boldsymbol{s}}_e}(t), \cdots ,\right. \left.{{\boldsymbol{s}}_{{N_{\text{t}}}}}(t) \right]^{\text{T}} \in {{\mathbf{C}}^{MN{N_{\text{t}}} \times 1}} $，其中$ {{\boldsymbol{s}}_e}(t) \in {{\mathbf{C}}^{MN \times 1}} $，按照$N$个时隙的顺序排列，每个时隙包含$M$个符号. n(t)为噪声向量，${\boldsymbol{n}}(t) = [{{\boldsymbol{n}}_1}(t),{{\boldsymbol{n}}_2}(t), \cdots ,{{\boldsymbol{n}}_g}(t), \cdots , {{\boldsymbol{n}}_{{N_{\text{r}}}}}(t)]^{\text{T}} \in {{\mathbf{C}}^{MN{N_{\text{r}}} \times 1}}$，$ {{\boldsymbol{n}}_g}(t) \in {{\mathbf{C}}^{MN \times 1}} $，它的元素服从均值为0，方差为${\sigma ^2}$的复高斯分布.

在接收端进行相对应的逆变换，将第$g$根接收天线收到的时域信号$ {{\boldsymbol{r}}_g}(t) $通过维格纳变换映射到${{\boldsymbol{Y}}_g}(t,f)$上，对${{\boldsymbol{Y}}_g}(t,f)$进行离散化采样处理得到第$g$根接收天线上第$n$个时隙和第$m$个子载波的时频域信号${{\boldsymbol{Y}}_g}[n,m]$：

(7)$ {{{{\boldsymbol{Y}}}_g}(t,f) = \displaystyle\int_{ - \infty }^\infty {g_{{\text{rx}}}^*\;(t - n{T_{\text{s}}}){{\boldsymbol{r}}_g}(t){\exp\;\left[{ - {\rm{j}}2{\text{π}} f(t - n{T_{\text{s}}})}\right]}{\text{d}}t} ,} $

(8)$ {{{\boldsymbol{Y}}}_g}[n,m] = {{{\boldsymbol{Y}}}_g}(t,f)\left| {_{t = n{T_{\text{s}}},f = m\Delta f}} \right.. $

式中：${g}_{\text{rx}}\;(\cdot )$表示接收脉冲，$ {(\cdot )}^{*} $表示取共轭操作. 随后通过辛有限傅里叶变换（symplectic finite Fourier transform, SFFT）将$ {{{\boldsymbol{Y}}}_g}[n,m] $恢复到时延-多普勒域：

(9)$ {{\boldsymbol{y}}_g}[k,l] = \frac{1}{{\sqrt {MN} }}\sum\limits_{n = 0}^{N - 1} {\sum\limits_{m = 0}^{M - 1} {{{\boldsymbol{Y}}_g}[n,m]{\exp\;\left[{ - {\rm{j}}2{\text{π}} \left( {\dfrac{{nk}}{N} - \dfrac{{ml}}{M}} \right)}\right]}} } . $

式中：$l \in \left\{ {1,2, \cdots ,M} \right\}$，$k \in \left\{ {1,2, \cdots ,N} \right\}$.

经上述过程后，时延多普勒域上的输入输出关系可以由如下公式表示：

(10)$ {\text{vec}}\;({{\boldsymbol{y}}}){\mathbf{ = }}{{{\boldsymbol{H}}}_{{\text{DD}}}}{\text{vec}}\;({{\boldsymbol{x}}}){\mathbf+}{{\boldsymbol{w}}}. $

式中：$ {{\boldsymbol{y}}} $表示所有接收天线上总的时延多普勒域接收信号；$ {{{\boldsymbol{H}}}_{{\text{DD}}}} $为时延多普勒域等效信道矩阵，可以由$ {{{\boldsymbol{H}}}_{{\text{DD}}}} = ({{{\boldsymbol{F}}}_N} \otimes {{{\boldsymbol{I}}}_M}){{\boldsymbol{H}}}({{\boldsymbol{F}}}_N^{\text{H}} \otimes {{{\boldsymbol{I}}}_M}) $^[19]得到，其中$ {{{\boldsymbol{F}}}_N} $为经过归一化处理的$N$点离散傅里叶变换矩阵，$ {(\cdot )}^{\text{H}} $表示共轭转置操作，$ \otimes $表示克罗内克积；${{\boldsymbol{x}}}$表示所有发射天线上总的时延多普勒域发射信号；${{\boldsymbol{w}}}$表示时延多普勒域噪声向量；${\text{vec}}\;( \cdot )$表示对矩阵列向量化 .

基于式（10）的系统模型，后验概率$ p\left( {{{\boldsymbol{x}}}\left| {{\boldsymbol{y}}} \right.} \right) $的因式分解为

(11)$ p\left( {{{\boldsymbol{x}}}\left| {{\boldsymbol{y}}} \right.} \right) = p\left( {{\boldsymbol{x}}} \right)p\left( {{{\boldsymbol{y}}}\left| {{\boldsymbol{x}}} \right.} \right) = \prod\limits_{i = 1}^{MN{N_{\text{r}}}} {{f_i}({y_i}\left| {{\boldsymbol{x}}} \right.)} \prod\limits_{j = 1}^{MN} {p\left( {{{{\boldsymbol{x}}}_j}} \right)} . $

式中：${y_i}$表示接收信号${{\boldsymbol{y}}}$的第$i$个符号，${{{\boldsymbol{x}}}_j}$表示发射信号${{\boldsymbol{x}}}$的第$j$个GSM符号矢量.

似然函数表达式为

(12)$ {f_i}({y_i}\left| {\boldsymbol{x}} \right.) \propto \exp \left( - \frac{1}{{\sigma _i^2}}\left|{y_i} - \sum\limits_{j = 1}^{MN{N_t}} {{h_{i,j,e}}{{{x}}_{j}{(e)}}} \right|\right). $

式中：$ {h_{i,j,e}} $表示第$j$个GSM发射符号上第$e$根天线的元素与接收信号第$i$个元素之间的响应系数，对应式（10）等价于${{{H}}_{{\text{DD}}}}[i,(e - 1)MN+j]$.

先验概率的表达式为

(13)$ p({{\boldsymbol{x}}_j}) = \frac{1}{{{C_{\text{a}}}}}\sum\limits_{{\boldsymbol{\beta }} \in {{{S}}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{{C}}}} {\delta ({{\boldsymbol{x}}_j} - {{\boldsymbol{\beta}} })} . $

式中：$\delta $为冲击响应函数.

变分推断是常见的计算边缘概率的有效方法，核心思想是选择一个简单的分布$q$来近似难以处理的复杂概率分布$p$. 2个分布越接近，则Kullback-Leibler（KL）散度越小，表示为

(14)$ {{{\rm{Proj}}} _\varPhi }(p) = \mathop {\arg \min }\limits_{q \in \varPhi }\; {D_{{\rm{KL}}}}(p||q). $

式中：${D_{{\rm{KL}}}}(p||q)$为KL散度，表示2个概率分布的差异. BP算法是通过变量节点与因子节点之间迭代传递符号置信度对因子图进行概率推理的有效方法.

如图2所示给出了GSM-OTFS系统EP-BP因子图^[20]示例，其中包含变量节点（variable nodes, VNs）${{\boldsymbol{x}}_j}$、因子节点（factor nodes, FNs）${y_i}$和先验因子节点（prior factor nodes, PNs）${\varphi _i}$. 由于发射符号${{\boldsymbol{\beta}} } \in {{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}$均为向量，因子图模型中为矢量变量节点，采取多元复高斯概率密度函数（probability density function, PDF）近似. 须注意的是，OTFS的时延多普勒域信道矩阵$ {{\boldsymbol{H}}_{{\text{DD}}}} $具有带状的稀疏结构，在$ {{\boldsymbol{H}}_{{\text{DD}}}} $的每一行和每一列中只有少量的非零元素^[21]，因此只有当矢量VN和FN之间的信道系数不为0时，两节点之间才存在连接.

令$\mu _{{{\boldsymbol{x}}_j} \to {y_i}}^t({{\boldsymbol{x}}_j})$表示第$t$次迭代时第$j$个矢量VN向第$i$个FN传递的消息，$\mu _{{y_i} \to {{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j})$则是FN向VN传递的消息. 由和积算法可得，消息传递规则如下：

(15)$ \mu _{{{\boldsymbol{x}}_j} \to {y_i}}^t({{\boldsymbol{x}}_j}) = {\mu _{{\varphi _j} \to {{\boldsymbol{x}}_j}}}({{\boldsymbol{x}}_j})\prod\limits_{i' \ne i} {\mu _{{{{y}}_{i'}} \to {{\boldsymbol{x}}_j}}^{t - 1}({{\boldsymbol{x}}_j})} , $

(16)$ \mu _{{y_i} \to {{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j}) = \sum\limits_{{\boldsymbol{x}}\backslash {{\boldsymbol{x}}_j}} {{f_i}\left( {{y_i}|{{\boldsymbol{x}}}} \right)} \prod\limits_{j' \ne j} {\mu _{{{\boldsymbol{x}}_{j'}} \to {y_i}}^{t - 1}({{\boldsymbol{x}}_{j'}})} . $

式中：$ {\mu _{{\varphi _j} \to {{\boldsymbol{x}}_j}}}({{\boldsymbol{x}}_j}) = p({{\boldsymbol{x}}_j}) $表示第$j$个发射符号向量的先验概率，$ {{\boldsymbol{x}}}\backslash {{\boldsymbol{x}}_j} $表示在$ {{\boldsymbol{x}}} $中除去$ {{\boldsymbol{x}}_j} $.

接下来以第$j$个发射符号向量的检测为切入点，在计算符号置信度时，第$j$个矢量VN会收集所有相连接的节点传递来的信息，即第$j$个PN以及该VN连接的所有FN. 符号置信度的计算为所有传递到此VN节点消息的乘积，归一化计算如下：

(17)$ \mu _{{{\boldsymbol{x}}_j}}^t\left( {{{\boldsymbol{x}}_j} = {{\boldsymbol{\beta}} }} \right) = \frac{{{\mu _{{\varphi _j} \to {{\boldsymbol{x}}_j}}}({{\boldsymbol{x}}_j} = {{\boldsymbol{\beta}} })\mu _{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t({{\boldsymbol{x}}_j} = {{\boldsymbol{\beta}} })}}{{\displaystyle\sum\limits_{{{\boldsymbol{\beta}} }' \in {{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}} {{\mu _{{\varphi _j} \to {{\boldsymbol{x}}_j}}}({{\boldsymbol{x}}_j} = {{\boldsymbol{\beta}} }')\mu _{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t({{\boldsymbol{x}}_j} = {{\boldsymbol{\beta}} }')} }}. $

式中：$ \mu _{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t({{\boldsymbol{x}}_j}) $是所有FN在上一次迭代传递给该VN的消息的乘积，即

(18)$ \mu _{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t({{\boldsymbol{x}}_j}) = \prod\limits_i {\mu _{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1}({{\boldsymbol{x}}_j})} . $

根据（15）、（17）、（18），可得如下消息传递规则：

(19)$ \mu _{{{\boldsymbol{x}}_j} \to {y_i}}^t({{\boldsymbol{x}}_j}) \propto \frac{{\mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j})}}{{\mu _{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1}({{\boldsymbol{x}}_j})}}. $

式（16）中的外部求和涉及到对除$ {{\boldsymbol{x}}_j} $以外的所有发射符号的符号集进行全局搜索，这将带来指数级的计算复杂度. 因此使用EP算法来降低计算复杂度，将$ {{\boldsymbol{x}}_j} $视为连续随机变量，并将其真实分布近似为多元复高斯PDF，消息$ \mu _{{{\boldsymbol{x}}_j} \to {y_i}}^t({{\boldsymbol{x}}_j}) $和$\mu^t_{{y_i} \to {{\boldsymbol{x}}_j}}\;({{\boldsymbol{x}}_j}) $分别近似为$\; {\rm{CN}}\;({{\boldsymbol{x}}_j};{{\boldsymbol{m}}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t,{{\boldsymbol{V}}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t) $和${\rm{CN}}\;({{\boldsymbol{x}}_j} $；$ {{\boldsymbol{m}}}_{{y_i} \to {{\boldsymbol{x}}_j}}^t,{{\boldsymbol{V}}}_{{y_i} \to {{\boldsymbol{x}}_j}}^t) $. 基于多元复高斯PDF乘法规则^[22]，消息$ \mu _{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t({{\boldsymbol{x}}_j}) $表达式如下：

(20)$ \begin{split} \mu _{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t({{\boldsymbol{x}}_j}) \propto &\exp \Big[ - {({{\boldsymbol{x}}_j} - {{\boldsymbol{m}}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t)^{{\rm{H}}} \times }\\&{({{\boldsymbol{V}}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t)^{ - 1}}({{\boldsymbol{x}}_j} - {{\boldsymbol{m}}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t)\Big].\end{split} $

从矢量VN传递到PN的协方差矩阵和均值向量分别计算为

(21)$ {\boldsymbol{V}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t = {\left[\sum\limits_i {{{\left({\boldsymbol{V}}_{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1}\right)}^{ - 1}}} \right]^{ - 1}}, $

(22)$ {\boldsymbol{m}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t = {\boldsymbol{V}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t\sum\limits_i {\left[{{({\boldsymbol{V}}_{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1})}^{ - 1}}{\boldsymbol{m}}_{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1}\right]} . $

由于发射符号的各个元素相互独立，协方差矩阵为对角矩阵，式（21）和（22）中均值向量$ {\boldsymbol{m}}_{{y_i} \to {{\boldsymbol{x}}_j}}^0 $的元素初始化为0，协方差矩阵$ {\boldsymbol{V}}_{{y_i} \to {{\boldsymbol{x}}_j}}^0 $的对角元素初始化为$+\infty $. 为了进一步降低复杂度，将符号置信度$ \mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j}) $近似为多元复高斯$\;\;\;\;{\rm{PDF}}\; \hat \mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j}) = {\rm{CN}}({{\boldsymbol{x}}_j};{\hat {\boldsymbol{m}}}_{{{\boldsymbol{x}}_j}}^t,{\hat {\boldsymbol{V}}}_{{{\boldsymbol{x}}_j}}^t)j^t)$. 符号置信度$ \mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j}) $和近似值$ \hat \mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j}) $之间的近似可以通过矩匹配原理得到，近似之后均值向量和协方差矩阵表达式为

(23)$ {\hat {\boldsymbol{m}}}_{{{\boldsymbol{x}}_j}}^t = \sum\limits_{{{\boldsymbol{\beta}} } \in {{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}} {{\boldsymbol{\beta }}\mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j} = {{\boldsymbol{\beta}} })} , $

(24)$ {\hat {\boldsymbol{V}}}_{{{\boldsymbol{x}}_j}}^t = \sum\limits_{{{\boldsymbol{\beta}} } \in {{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}} {({{\boldsymbol{\beta}} } - {\hat {\boldsymbol{m}}}_{{{\boldsymbol{x}}_j}}^t){{({{\boldsymbol{\beta}} } - {\hat {\boldsymbol{m}}}_{{{\boldsymbol{x}}_j}}^t)}^{{\rm{H}}} }\mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j} = {{\boldsymbol{\boldsymbol{\beta }}}})} . $

由此复杂的消息传递规则可以由多元复高斯PDF表示为

(25)$ {\rm{CN}}({{\boldsymbol{x}}_j};{\boldsymbol{m}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t,{\boldsymbol{V}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t) \propto \frac{{{\rm{CN}}({{\boldsymbol{x}}_j};{\hat {\boldsymbol{m}}}_{{{\boldsymbol{x}}_j}}^t,{\hat {\boldsymbol{V}}}_{{{\boldsymbol{x}}_j}}^t)}}{{{\rm{CN}}({{\boldsymbol{x}}_j};{\boldsymbol{m}}_{{y_i} \to {{\boldsymbol{x}}_j}}^t,{\boldsymbol{V}}_{{y_i} \to {{\boldsymbol{x}}_j}}^t)}}. $

从矢量VN传递到FN的协方差矩阵和均值向量表达式分别为

(26)$ {\boldsymbol{V}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t = {\left[{\left({\hat {\boldsymbol{V}}}_{{{\boldsymbol{x}}_j}}^t\right)^{ - 1}} - {\left({\boldsymbol{V}}_{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1}\right)^{ - 1}}\right]^{ - 1}}, $

(27)$ {\boldsymbol{m}_{{\boldsymbol{x}}_j \rightarrow y_i}^t=\boldsymbol{V}_{{\boldsymbol{x}}_j \rightarrow y_i}^t\left[\left(\hat{\boldsymbol{V}}_{{\boldsymbol{x}}_j}^t\right)^{-1} \hat{\boldsymbol{m}}_{{\boldsymbol{x}}_j}^t-\left(\boldsymbol{V}_{y_i \rightarrow {\boldsymbol{x}}_j}^{t-1}\right)^{-1} \boldsymbol{m}_{y_i \rightarrow {\boldsymbol{x}}_j}^{t-1}\right].} $

由于消息$ \mu _{{{\boldsymbol{x}}_j} \to {y_i}}^t({{\boldsymbol{x}}_j}) $和$ \mu _{{y_i} \to {{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j}) $被近似为多元复高斯PDF，可得从FN传递到矢量VN的均值向量和协方差矩阵表达式如下：

(28)$ {{m}}_{{y_i} \to {{\boldsymbol{x}}_j}}^t(e) = \frac{1}{{{h_{i,j,e}}}}\left[{y_i} - \sum\limits_{(j',e') \ne (j,e)} {{h_{i,j',e'}}{{m}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t(e')} \right], $

(29)$ {{V}}_{{y_i} \to {{\boldsymbol{x}}_j}}^t(e) = \frac{1}{{{{\left| {{h_{i,j,e}}} \right|}^2}}}\left[\sum\limits_{(j',e') \ne (j,e)} {{{\left| {{h_{i,j',e'}}} \right|}^2}{{V}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t(e')+\sigma _i^2} \right]. $

式中：$e = 1, \cdots ,{N_{\text{t}}}$，$e' = 1, \cdots ,{N_{\text{t}}}$，$ {{m}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t(e) $和$ {{V}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t(e) $分别表示$ {\boldsymbol{m}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t $的第$e$个元素和$ {\boldsymbol{V}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t $的第$e$个对角元素，$ (j',e') \ne (j,e) $表示$j' \ne j$且$e' \ne e$，其中$j = 1, \cdots ,({N_{\text{t}}} - 1)MN+1, \cdots ,{N_{\text{t}}}MN$. 根据第$t$次迭代得到的符号置信度$ \mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j}) $，第$j$个发射符号的激活天线组合索引概率表达式为

(30)$ P_j^{{C}}\left( c \right) = \sum\limits_{\boldsymbol{\beta }} {\mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j} = {{\boldsymbol{\beta}} })} . $

式中：${{\boldsymbol{\beta}} } \in {{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}\left\{ {\left( {c - 1} \right){{({C_{\text{m}}})}^{{N_{\text{a}}}}}+1:c{{({C_{\text{m}}})}^{{N_{\text{a}}}}}} \right\}$，$c = 1, \cdots , {C_{\text{a}}}$，${{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}\left\{ {\left( {c - 1} \right){{({C_{\text{m}}})}^{{N_{\text{a}}}}}+1:c{{({C_{\text{m}}})}^{{N_{\text{a}}}}}} \right\}$表示发射符号集${{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}$中第$\left( {c - 1} \right){({C_{\text{m}}})^{{N_{\text{a}}}}}+1$到第$c{({C_{\text{m}}})^{{N_{\text{a}}}}}$个列向量. 激活天线组合索引的估计值表达式为${\hat c_j} = \mathop {\arg \max }\limits_{c = 1, \cdots ,{C_{\text{a}}}}\; [P_j^{{C}}(c)]$，其对应的激活天线组合为${{\boldsymbol{A}}}$中第${\hat c_j}$个列向量. 因此，第$j$个发射符号向量的估计向量为

(31)$ {\hat {\boldsymbol{x}}}_j^t = \sum\limits_{\boldsymbol{\beta }} {{{\boldsymbol{\beta}} }\mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j} = {{\boldsymbol{\beta}} })} . $

式中：${\boldsymbol{\beta }} \in {{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}\left\{ {\left( {c - 1} \right){{({C_{\text{m}}})}^{{N_{\text{a}}}}}+1:c{{({C_{\text{m}}})}^{{N_{\text{a}}}}}} \right\}$. 所有发射符号向量可以表示为${{\hat {\boldsymbol{x}}}^t} = [{({\hat {\boldsymbol{x}}}_1^t)^{\text{T}}}, \cdots , {({\hat {\boldsymbol{x}}}_j^t)^{\text{T}}}, \cdots , {({\hat {\boldsymbol{x}}}_{MN}^t)^{\text{T}}}]^{\text{T}}$.

EP-BP算法的迭代终止原则为$\left\| {{{\hat {\boldsymbol{x}}}}^t} -\right. \left.{{{\hat {\boldsymbol{x}}}}^{t - 1}} \right\|^2 < \varepsilon {\left\| {{{{\hat {\boldsymbol{x}}}}^{t - 1}}} \right\|^2}$，其中初始值的设定为$\;\;\;{{\hat {\boldsymbol{x}}}^0} = {{\mathbf{0}}_{{N_{\text{t}}}MN \times 1}}$，误差容限设为$0 < \varepsilon \ll 1.0$. 最大迭代次数为${T_{\max }}$，当$t = {T_{\max }}$，即达到最大迭代次数时，依然无法满足终止原则${\left\| {{{{\hat {\boldsymbol{x}}}}^t} - {{{\hat {\boldsymbol{x}}}}^{t - 1}}} \right\|^2} < \varepsilon {\left\| {{{{\hat {\boldsymbol{x}}}}^{t - 1}}} \right\|^2}$，此时也终止算法迭代. GSM-OTFS系统的EP-BP信号检测算法的详细过程如下.

算法1　　GSM-OTFS系统EP-BP信号检测算法

1. 输入：$ \small{{{\boldsymbol{y}}} 、{{\boldsymbol{H}}}、\sigma^2、 {\mu _{{\varphi _j} \to {{\boldsymbol{x}}_j}}}({{\boldsymbol{x}}_j}) = 1/{C_{\text{c}}} 、{T_{\max }}、\varepsilon} $

2. 初始化：$ \small{\boldsymbol{m}_{y_i \rightarrow {\boldsymbol{x}}_j}^0(e)=0 、\boldsymbol{V}_{{\boldsymbol{x}}_j \rightarrow \varphi_j}^t(e) \rightarrow \infty、 {\hat {\boldsymbol{x}}}^0=} $$ \small{{{\mathbf{0}}_{{N_{\text{t}}}MN \times 1}}} $

3. for $\small{t = 1:{T_{\max }}}$ (迭代次数)

4. 　for $\small{j = 1:MN}$

5. 　$\small{{\boldsymbol{V}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t = {\left[\sum\limits_i {{{({\boldsymbol{V}}_{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1})}^{ - 1}}} \right]^{ - 1}}}$

6. 　$\small{{\boldsymbol{m}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t = {\boldsymbol{V}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t\sum\limits_i {\left[{{({\boldsymbol{V}}_{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1})}^{ - 1}}{\boldsymbol{m}}_{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1}\right]}}$

7. 　$\small{\mu _{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t({{\boldsymbol{x}}_j}) \propto \exp \;[ - {({{\boldsymbol{x}}_j} - {\boldsymbol{m}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t)^{{\mathrm{H}}} }{({{\boldsymbol{V}}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t)^{ - 1}}}$$ \small{({{\boldsymbol{x}}_j} - {\boldsymbol{m}}_{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t)]}$

$\small{8.\;\quad \mu _{{{\boldsymbol{x}}_j}}^t\left( {{{\boldsymbol{x}}_j} = {\boldsymbol{\beta }}} \right) = \frac{{{\mu _{{\varphi _j} \to {{\boldsymbol{x}}_j}}}({{\boldsymbol{x}}_j} = {\boldsymbol{\beta }})\mu _{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t({{\boldsymbol{x}}_j} = {\boldsymbol{\beta }})}}{{\displaystyle\sum\limits_{{\mathbf{\boldsymbol{\beta }}}' \in {{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}} {{\mu _{{\varphi _j} \to {{\boldsymbol{x}}_j}}}({{\boldsymbol{x}}_j} = {\boldsymbol{\beta }}')\mu _{{{\boldsymbol{x}}_j} \to {\varphi _j}}^t({{\boldsymbol{x}}_j} = {\boldsymbol{\beta }}')} }}}$

$\small{9.\;\quad {\hat {\boldsymbol{m}}}_{{{\boldsymbol{x}}_j}}^t = \displaystyle\sum\limits_{{\mathbf{\boldsymbol{\beta }}} \in {{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}} {{\boldsymbol{\beta }}\mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j} = {\boldsymbol{\beta }})}}$

$\small{10.\;\;\; {\hat {\boldsymbol{V}}}_{{{\boldsymbol{x}}_j}}^t = \displaystyle\sum\limits_{{\boldsymbol{\beta }} \in {{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}} {({\boldsymbol{\beta }} - {\hat {\boldsymbol{m}}}_{{{\boldsymbol{x}}_j}}^t){{({\boldsymbol{\beta }} - {\hat {\boldsymbol{m}}}_{{{\boldsymbol{x}}_j}}^t)}^{{\rm{H}}} }\mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j} = {\boldsymbol{\beta }})}}$

11. for $\small{i = 1:{N_t}MN}$

12. 　　$\small{{\boldsymbol{V}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t = {[{({\hat {\boldsymbol{V}}}_{{{\boldsymbol{x}}_j}}^t)^{ - 1}} - {({\boldsymbol{V}}_{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1})^{ - 1}}]^{ - 1}}}$

13. 　　$\small{{\boldsymbol{m}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t = {\boldsymbol{V}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t\left[{({\hat {\boldsymbol{V}}}_{{{\boldsymbol{x}}_j}}^t)^{ - 1}}{\hat {\boldsymbol{m}}}_{{{\boldsymbol{x}}_j}}^t - {({\boldsymbol{V}}_{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1})^{ - 1}}{\boldsymbol{m}}_{{y_i} \to {{\boldsymbol{x}}_j}}^{t - 1}\right]}$

14. 　end for

15. end for

16. for $\small{i = 1:{N_t}MN}$

17. 　for $\small{j = 1:MN}$

$\small{18.\;\quad\quad{{m}}_{{y_i} \to {{\boldsymbol{x}}_j}}^t(e) = \dfrac{1}{{{h_{i,j,e}}}}\left[{y_i} - \displaystyle\sum\limits_{(j',e') \ne (j,e)} {{h_{i,j',e'}}{{m}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t(e')} \right]}$

$\small{19.\;\quad\quad{{V}}_{{y_i} \to {{\boldsymbol{x}}_j}}^t(e) = \dfrac{1}{{{{\left| {{h_{i,j,e}}} \right|}^2}}}\left[\displaystyle\sum\limits_{(j',e') \ne (j,e)} {{{ \left| {{h_{i,j',e'}}} \right| }^2}{{V}}_{{{\boldsymbol{x}}_j} \to {y_i}}^t(e') + \sigma _i^2} \right]}$

20. 　end for

21. end for

22. for $\small{j = 1:MN}$

$\small{23.\;\qquad P_j^{{C}}\left( c \right) = \displaystyle\sum\limits_{\boldsymbol{\beta }} {\mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j} = {\boldsymbol{\beta }})}且}$$\small{{\boldsymbol{\beta }} \in {{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}\left\{ {\left( {c - 1} \right){{({C_{\text{m}}})}^{{N_{\text{a}}}}}+1:c{{({C_{\text{m}}})}^{{N_{\text{a}}}}}} \right\}}$

$\small{24.\;\qquad{\hat c_j} = \mathop {\arg \max }\limits_{c = 1, \cdots ,{C_{\text{a}}}} \;[P_j^{{C}}(c)]}$

$\small{25.\;\qquad{{\hat{\boldsymbol{ x}}}}_j^t = \displaystyle\sum\limits_{\boldsymbol{\beta }} {{\boldsymbol{\beta }}\mu _{{{\boldsymbol{x}}_j}}^t({{\boldsymbol{x}}_j} = {\boldsymbol{\beta }})}且}$$\small{{\boldsymbol{\beta }} \in {{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}\left\{ {\left( {c - 1} \right){{({C_{\text{m}}})}^{{N_{\text{a}}}}}+1:c{{({C_{\text{m}}})}^{{N_{\text{a}}}}}} \right\}}$

26. 　end for

27. 　全部检测符号：$\small{{{\hat {\boldsymbol{x}}}^t} = {[{({\hat {\boldsymbol{x}}}_1^t)^{\text{T}}}, \cdots ,{({\hat {\boldsymbol{x}}}_j^t)^{\text{T}}}, \cdots ,{({\hat {\boldsymbol{x}}}_{MN}^t)^{\text{T}}}]^{\text{T}}}}$

28. 　终止原则：$\small{t = {T_{\max }}}$或者$\small{{\left\| {{{{\hat {\boldsymbol{x}}}}^t} - {{{\hat {\boldsymbol{x}}}}^{t - 1}}} \right\|^2} < \varepsilon {\left\| {{{{\hat {\boldsymbol{x}}}}^{t - 1}}} \right\|^2}}$

29. end for

30. 输出：$\small{{{\hat {\boldsymbol{x}}}} = {{{\hat {\boldsymbol{x}}}}^t}}$

EP-BP算法（式（17））的符号置信度计算须遍历发射符号集合${{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}$，因此计算复杂度较高. 为了降低EP-BP算法的复杂度，提出分步式的两阶段TS-EP-BP算法. 在第1阶段，执行$T_{\max }^{\text{b}}$次EP-BP算法的迭代（$ T_{\max }^{\text{b}} < {T_{\max }} $），得到激活天线组合的索引估计值. 因此，发射符号向量中的星座符号可以视为独立地分布在激活天线上，在符号置信度的计算中仅须遍历星座调制符号图，复杂度大大降低，以${{S}}_{4,2}^{{{C}_{{\text{QPSK}}}}}$为例，EP-BP算法须遍历的发射符号集合有64个符号，而TS-EP-BP的第2阶段仅须遍历QPSK星座图的4个符号. 在第2阶段每个矢量VN可以分解为${N_{\text{a}}}$个子VN（child variable nodes，cVNs），通过激活天线组合的索引估计值判断出无效的子VN（发射该GSM符号时此天线未被激活，则对应子VN无效）. 将无效的子VN从因子图中裁去，得到新的因子图. 此时可由多元复高斯PDF近似转为一元复高斯PDF近似. GSM-OTFS系统TS-EP-BP因子图如图3所示.

基于第1阶段结束时的符号置信度$ \mu _{{{\boldsymbol{x}}_j}}^{T_{\max }^{\text{b}}}({{\boldsymbol{x}}_j}) $，由式（30）可以计算出第$j$个发射符号的激活天线组合索引概率，由最大概率确定激活天线组合. 裁剪后的因子图剩余${N_{\text{a}}}MN$个子VN.

在TS-EP-BP算法的第2阶段，FN传递到子VN的均值与方差初值设定为$ m_{{y_i} \to {x_k}}^{T_{\max }^{\text{b}}} = {{m}}_{{y_i} \to {{{x}}_j}}^{T_{\max }^{\text{b}}}(e) $和$v_{{y_i} \to {x_k}}^{T_{\max }^{\text{b}}} = {v}_{{y_i} \to {{{x}}_j}}^{T_{\max }^{\text{b}}}(e)$. 后续迭代中一元复高斯消息$\mu _{{x_k} \to {\varphi _k}}^t({x_k}) = \prod\limits_i {\mu _{{y_i} \to {x_k}}^{t - 1}({x_k})} = {\rm{CN}}\;({x_k};m_{{x_k} \to {\varphi _k}}^t,v_{{x_k} \to {\varphi _k}}^t)\;$的方差和均值表达式为

(32)$ v_{{x_k} \to {\varphi _k}}^t = {\left(\sum\limits_i {\frac{1}{{v_{{y_i} \to {x_k}}^{t - 1}}}} \right)^{ - 1}}, $

(33)$ m_{{x_k} \to {\varphi _k}}^t = v_{{x_k} \to {\varphi _k}}^t\sum\limits_i {\frac{{m_{{y_i} \to {x_k}}^{t - 1}}}{{v_{{y_i} \to {x_k}}^{t - 1}}}} . $

高斯消息$\mu _{{x_k} \to {\varphi _k}}^t\left( {{x_k}} \right)$的PDF表达式为

(34)$ \mu _{{x_k} \to {\varphi _k}}^t\left( {{x_k}} \right) \propto \exp\; \left( - \frac{1}{{v_{{x_k} \to {\varphi _k}}^t}}{\left| {{x_k} - m_{{x_k} \to {\varphi _k}}^t} \right|^2}\right). $

符号置信度$\mu _{{x_k} \to {\varphi _k}}^t\left( {{x_k} = \alpha } \right)$可以归一化计算为

(35)$ \begin{split} &\mu _{{x_k}}^t\left( {{x_k} = \alpha } \right) = \\&\qquad \frac{{{\mu _{{\varphi _k} \to {x_k}}}\left( {{x_k} = \alpha } \right)\mu _{{x_k} \to {\varphi _k}}^t\left( {{x_k} = \alpha } \right)}}{{\displaystyle\sum\limits_{\alpha ' \in {C}} {{\mu _{{\varphi _k} \to {x_k}}}\left( {{x_k} = \alpha '} \right)\mu _{{x_k} \to {\varphi _k}}^t\left( {{x_k} = \alpha '} \right)} }}.\end{split} $

其中，先验概率$ {\mu _{{\varphi _k} \to {x_k}}}\left( {{x_k} = \alpha } \right) $可以定义为

(36)$ {\mu _{{\varphi _k} \to {x_k}}}\left( {{x_k}} \right) = \frac{1}{{{C_{\text{m}}}}}\sum\limits_{\alpha \in {C}} {\delta \left( {{x_k} - \alpha } \right)} . $

此时每个符号置信度的计算次数由$ {C_{\text{c}}} = \left| {{{S}}_{{N_{\text{t}}},{N_{\text{a}}}}^{{C}}} \right| $降低为${C_{\text{m}}} = \left| {{{C}}} \right|$. 将$ \mu _{{x_k}}^t\left( {{x_k}} \right) $近似为一元复高斯分布$\mu _{{x_k}}^t({x_k}) = {\rm{CN}}({x_k};m_{{x_k}}^t,v_{{x_k}}^t)$，其均值和方差表达式为

(37)$ \hat m_{{x_k}}^t = \sum\limits_{\alpha \in {C}} {\alpha \mu _{{x_k}}^t\left( {{x_k} = \alpha } \right)} , $

(38)$ \hat v_{{x_k}}^t = \sum\limits_{\alpha \in {C}} {{{\left| {\alpha - \hat m_{{x_k}}^t} \right|}^2}\mu _{{x_k}}^t\left( {{x_k} = \alpha } \right)} . $

与EP-BP算法类似，TS-EP-BP算法消息传递规则如下：

(39)$ \mu _{{x_k} \to {y_i}}^t\left( {{x_k}} \right) \propto \frac{{\mu _{{x_k}}^t\left( {{x_k}} \right)}}{{\mu _{{y_i} \to {x_k}}^{t - 1}\left( {{x_k}} \right)}}. $

将其近似为高斯PDF形式，可得

(40)$ {\rm{CN}}\;({x_k};m_{{x_k} \to {y_i}}^t,v_{{x_k} \to {y_i}}^t) \propto \frac{{{\rm{CN}}\;({x_k};m_{{x_k}}^t,v_{{x_k}}^t)}}{{{\rm{CN}}\;({x_k};m_{{y_i} \to {x_k}}^{t - 1},v_{{y_i} \to {x_k}}^{t - 1})}}. $

从子VN传递到FN的均值和方差表达式如下：

(41)$ v_{{x_k} \to {y_i}}^t = {\left(\frac{1}{{\hat v_{{x_k}}^t}} - \frac{1}{{v_{{y_i} \to {x_k}}^{t - 1}}}\right)^{ - 1}}, $

(42)$ m_{{x_k} \to {y_i}}^t = v_{{x_k} \to {y_i}}^t\left(\frac{{\hat m_{{x_k}}^t}}{{\hat v_{{x_k}}^t}} - \frac{{m_{{y_i} \to {x_k}}^{t - 1}}}{{v_{{y_i} \to {x_k}}^{t - 1}}}\right). $

从FN到子VN的均值和方差表达式如下：

(43)$ m_{{y_i} \to {x_k}}^t = \frac{1}{{{h_{i,k}}}}\left({y_i} - \sum\limits_{k' \ne k} {{h_{i,k'}}m_{{x_{k'}} \to {y_i}}^t} \right), $

(44)$ v_{{y_i} \to {x_k}}^t = \frac{1}{{{{\left| {{h_{i,k}}} \right|}^2}}}\left(\sum\limits_{k' \ne k} {{{\left| {{h_{i,k}}} \right|}^2}v_{{x_k} \to {y_i}}^t} +\sigma^2\right). $

式中：h_i,k表示发射符号里第k个子VN符号与接收信号第i个元素之间的相应系数.

根据第$t$次迭代得到的符号置信度$\mu _{{x_k}}^t\left( {{x_k}} \right)$，第$k$个子VN的符号估计值表达式如下：

(45)$ \hat x_k^t = \sum\limits_{\alpha \in C} {\alpha \mu _{{x_k}}^t\left( {{x_k} = \alpha } \right)} . $

根据式（45）以及激活天线组合索引，可以得到完整的符号估计值${{\hat x}^t}$，无效子VN的值为0. 第1阶段得到激活天线组合估计值的准确性，会对整体的性能有直接影响，将$T_{\max }^{\text{b}}$调整为更大的数可以提升估计准确性进而提高误码率，但同时会增加计算复杂度，这表明所提算法可以在复杂度与性能之间实现灵活平衡. 当$T_{\max }^{\text{b}} = {T_{\max }}$时，TE-EP-BP算法与EP-BP算法等价. TS-EP-BP信号检测算法的过程总结如下.

算法2　　GSM-OTFS系统TS-EP-BP信号检测算法

1. 输入：$\small{{{\boldsymbol{y}}} 、{{\boldsymbol{H}}}、\sigma^2、 {\mu _{{\varphi _k} \to {x_k}}}({{{x}}_k}) = 1/{C_{\text{m}}}}、{T_{\max }}、T_{\max }^{\rm{b}}$$\small{和\varepsilon}$

2. for $\small{t = 1:T_{\max }^{\rm{b}}}$

3. 　执行算法1 EP-BP信号检测算法

4. end for

5. 初始化：$\small{m_{{y_i} \to {x_k}}^{T_{\max }^{\text{b}}} = {{m}}_{{y_i} \to {{{x}}_j}}^{T_{\max }^{\text{b}}}(e)}$和$\small{v_{{y_i} \to {x_k}}^{T_{\max }^{\text{b}}} = {v}_{{y_i} \to {{{x}}_j}}^{T_{\max }^{\text{b}}}(e)}$

6. for $\small{t = T_{\max }^{\rm{b}}+1:{T_{\max }}}$

7. 　for $\small{k = 1:{N_{\text{a}}}MN}$

$\small{8.\;\;\qquad v_{{x_k} \to {\varphi _k}}^t = {\left(\sum\limits_i {\frac{1}{{v_{{y_i} \to {x_k}}^{t - 1}}}} \right)^{ - 1}},}$

$\small{9.\;\;\qquad m_{{x_k} \to {\varphi _k}}^t = v_{{x_k} \to {\varphi _k}}^t\sum\limits_i {\frac{{m_{{y_i} \to {x_k}}^{t - 1}}}{{v_{{y_i} \to {x_k}}^{t - 1}}}}}$

$\small{10.\;\qquad\mu _{{x_k} \to {\varphi _k}}^t\left( {{x_k}} \right) \propto \exp \left( - \frac{1}{{v_{{x_k} \to {\varphi _k}}^t}}{\left| {{x_k} - m_{{x_k} \to {\varphi _k}}^t} \right|^2}\right).}$

$\small{11.\;\qquad \mu _{{x_k}}^t\left( {{x_k} = \alpha } \right) = \frac{{{\mu _{{\varphi _k} \to {x_k}}}\left( {{x_k} = \alpha } \right)\mu _{{x_k} \to {\varphi _k}}^t\left( {{x_k} = \alpha } \right)}}{{\sum\limits_{\alpha ' \in {C}} {{\mu _{{\varphi _k} \to {x_k}}}\left( {{x_k} = \alpha '} \right)\mu _{{x_k} \to {\varphi _k}}^t\left( {{x_k} = \alpha '} \right)} }}}$

$\small{12.\;\qquad \hat m_{{x_k}}^t = \sum\limits_{\alpha \in {C}} {\alpha \mu _{{x_k}}^t\left( {{x_k} = \alpha } \right)}}$

$\small{13.\;\qquad \hat v_{{x_k}}^t = \sum\limits_{\alpha \in {C}} {{{\left| {\alpha - \hat m_{{x_k}}^t} \right|}^2}\mu _{{x_k}}^t\left( {{x_k} = \alpha } \right)}}$

14. 　　 for $\small{i = 1:{N_t}MN}$

$\small{15.\;\qquad\quad v_{{x_k} \to {y_i}}^t = {\left(\frac{1}{{\hat v_{{x_k}}^t}} - \frac{1}{{v_{{y_i} \to {x_k}}^{t - 1}}}\right)^{ - 1}},}$

$\small{16.\;\qquad\quad m_{{x_k} \to {y_i}}^t = v_{{x_k} \to {y_i}}^t\left(\frac{{\hat m_{{x_k}}^t}}{{\hat v_{{x_k}}^t}} - \frac{{m_{{y_i} \to {x_k}}^{t - 1}}}{{v_{{y_i} \to {x_k}}^{t - 1}}}\right).}$

17. 　　 end for

18. 　end for

19. 　for $\small{i = 1:{N_{\mathrm{t}}}MN}$

20. 　　for $\small{k = 1:{N_{\mathrm{a}}}MN}$

$\small{21.\;\qquad \;\; m_{{y_i} \to {x_k}}^t = \frac{1}{{{h_{i,k}}}}\left({y_i} - \sum\limits_{k' \ne k} {{h_{i,k'}}m_{{x_{k'}} \to {y_i}}^t} \right),}$

$\small{22.\;\qquad \;\; v_{{y_i} \to {x_k}}^t = \frac{1}{{{{\left| {{h_{i,k}}} \right|}^2}}}\left(\sum\limits_{k' \ne k} {{{\left| {{h_{i,k}}} \right|}^2}v_{{x_k} \to {y_i}}^t} +\sigma^2\right).}$

23. 　　end for

24. 　end for

25. 　for $\small{k = 1:{N_{\text{a}}}MN}$

$\small{26.\;\qquad \;\; \hat x_k^t = \sum\limits_{\alpha \in {C}} {\alpha \mu _{{x_k}}^t\left( {{x_k} = \alpha } \right)}}$

27. 　end for

28. 　全部检测符号：$\small{{{\hat {\boldsymbol{x}}}^t}}$

29. 　终止原则：$\small{t = {T_{\max }} } $或者$\small{{\left\| {{{{\hat {\boldsymbol{x}}}}^t} - {{{\hat {\boldsymbol{x}}}}^{t - 1}}} \right\|^2} < \varepsilon {\left\| {{{{\hat {\boldsymbol{x}}}}^{t - 1}}} \right\|^2}}$

30. end for

31. 输出：$\small{{{\hat {\boldsymbol{x}}}} = {{{\hat {\boldsymbol{x}}}}^t}}$

令$T$表示算法的迭代次数，${T^{\text{b}}}$表示TS-EP-BP算法第1阶段的迭代次数. 所提算法、GAMP^[14]算法以及不考虑GSM约束的MP^[21]算法计算复杂度总结如表1所示.

为了验证和对比改进后的分阶段TS-EP-BP算法在高速移动场景下的信号检测性能，对GSM-OTFS系统进行性能仿真，收发端各配备4根天线即${N_{\text{t}}} = {N_{\text{r}}} = 4$，激活天线数目$ {N_{\text{a}}} = 2 $，迭代算法设计最大迭代次数${T_{\max }} = 15$，误差容限设为$\varepsilon = {10^{ - 3}}$，其余仿真参数如表2所示.

由于理想脉冲难以物理实现，发送端与接收端均采用矩形窗进行脉冲成形^[4]. 不同路径的时延抽头采取LTE的扩展车辆A信道模型，每个时延抽头都对应一个使用Jakes公式产生的多普勒抽头^[21]. 蒙特卡洛仿真次数为1000次，每次发射一个OTFS帧，一帧包含$MN$个比特.

为了证明算法的收敛性，如图4所示给出EP-BP算法的误码率R_ber随迭代次数t_i变化的曲线. 其中，$M = N = 16$，信噪比为15 dB. 可知迭代6次后算法收敛.

如图5所示给出不同终端移动速度v下GSM-OTFS系统EP-BP检测算法的误码率性能随信噪比R_snr的变化曲线. 其中，$M = N = 16$. 可以看出，随着车速的增加，OTFS系统的性能随之提高，这是因为更高的车速提供了更高的多普勒增益.

如图6所示给出各算法的误码率性能比较. 其中，$M = N = 16$，终端移动速度为500 km/h. 可以看出，EP-BP算法具有最优的误码率性能，当$ T_{\max }^{\text{b}} = 6 $时，TS-EP-BP算法可以得到近似EP-BP算法的误码率性能，而当$ T_{\max }^{\text{b}} = 3 $时，TS-EP-BP算法则有一定的误码率性能损失，可以设定更大的$ T_{\max }^{\text{b}} $，从而在牺牲一定计算复杂度的情况下得到更近似于EP-BP算法的误码率性能. 在与图6相同的仿真条件设置下，通过2.3节可以计算出EP-BP算法迭代一轮的复杂度约为$6.9 \times {10^5}$，TS-EP-BP算法第2阶段时迭代一轮的复杂度约为$4.7 \times {10^4}$，相当于EP-BP算法的6.84%，可见TS-EP-BP算法在降低复杂度方面的优越性. GAMP^[14]算法目前是GSM-OTFS系统较优的信号检测算法，从仿真结果来看本研究所提算法的误码率性能与之相近，相同条件下该算法迭代一轮的复杂度约为$4.5 \times {10^6}$，考虑到GAMP算法仅需5次迭代，而EP-BP算法最优性能需要15次迭代，EP-BP总的算法复杂度约为GAMP的46.1%. 而文献[21]中的MP算法由于不适应GSM结构的部分天线激活特殊性，在增加了GSM约束情况下，其误码率性能仅仅略优于MMSE算法的.

如图7所示为$M = N = 32$时，不同算法的误码率性能与信噪比的关系，可知当$MN$增大时，时延多普勒域分辨率变高，整体误码率性能变好. 但相对于本研究算法，MP算法的误码率性能依然较差. 此外，TS-EP-BP算法需要更大的$ T_{\max }^{\text{b}} $来趋近于EP-BP算法的误码率性能. 例如，当$M = N = 32$时，TS-EP-BP需要设定$ T_{\max }^{\text{b}} = 9 $才能达到和EP-BP类似的误码率性能损失.

针对GSM-OTFS系统提出混合期望传播与置信传播的EP-BP联合信号检测算法，并通过减少迭代中符号置信度计算步骤需要遍历的发射符号集，进一步提出了低复杂度TS-EP-BP算法. 仿真结果表明，TS-EP-BP检测算法能以更低的复杂度取得与EP-BP相近的性能. 另外，所提出的TS-EP-BP算法可以调整第1阶段的迭代次数来实现误码率性能和复杂度之间的灵活权衡. 本研究的GSM模型是以4收4发为例，规模较小，当天线规模提升之后，激活天线组合集的复杂度会指数上升，未来应当关注这一难点，研究适应大规模天线的算法.