一种快变OFDM系统的渐进迭代时变信道估计和ICI消除方法

摘要

一种快变OFDM系统的渐进迭代时变信道估计和ICI消除方法，用于无线通信信道估计领域，用于快变信道产生ICI严重影响情况下，基于导频的信道估计。其特征在于：将载波间干扰ICI与信道噪声之和(Sum of ICI and channel noise，SIN)作为Kalman滤波器的去噪对象，消除Kalman估计时ICI的影响。此外，用于迭代的数据采用渐进增长的方式，用于测量的子载波数量，在迭代的过程中沿各个导频点两侧缓慢增长，从而抑制由ICI产生的影响。在SNR较小的情况下，该方案的性能比起现有算法有显著的提升。

说明书

技术领域

本发明涉及一种渐进的带有ICI消除的OFDM迭代信道估计方法。属于无线通信中信道估计研究的相关领域。 

技术背景

正交频分复用(Orthogonal Frequency Division Multiplexing,OFDM)是多载波调制的一种，其思想是将信道分成若干个相互正交的子载波，将高速的串行数据转换成并行的多组低速数据流，分别调至到这些子载波上传输，以提高频谱利用率。若子载波上的信号带宽小于信道的相关带宽，则子载波可以看成平坦性衰落，从而消除符号间干扰。由于其高频谱利用率和良好的抗干扰能力，OFDM技术已经被广泛应用于广播式的音频、视频领域和民用通信系统，主要的应用包括：非对称的数字用户环路（ADSL）、ETSI标准的数字音频广播（DAB）、数字视频广播（DVB）、高清晰度电视（HDTV）、无线局域网（WLAN）等。 

OFDM技术会同时受到由信道多径时延引起的频率选择特性衰落，和由信道的多普勒扩展引起的时间选择性衰落的影响，使系统性能下降。频率选择特性导致接收信号的到达时间变化从而影响到其幅度和相位。时间选择特性导致OFDM系统子载波之间的正交性受到影响，造成子载波间的干扰(intersubcarrier interference,ICI)，导致系统的性能下降。特别是在高速移动情况下，一个OFDM符号中信道也会发生很大变化，ICI的影响会更加严重。 

为此，目前有多种信道估计算法，通过插入导频和差值的方法进行估计。但是这些算法在高速移动的OFDM系统中可能不适用。近期，基扩展模型（Basis Expansion Model，BEM）算法被广泛用于模拟时频双选信道，根据所使用的基的不同，可以分为复指数基BEM算法（CE-BEM），离散余弦变换BEM（DTC-BEM），多项式基BEM（P-BEM），离散扩展球状BEM（DPS-BEM）以及离散Karhuen-LoèveBEM（KL-BEM）。在这些算法之中，P-BEM算法性能最好。为了降低ICI，自消除技术被提出。通过将信息映射到一组子载波上，产生ICI自消除，但会引起频谱效率的降低。另有算法将数据检测加入到信道估计中，数据检测用于信道估计中的迭代算法和数据恢复，从而改善估计的效果。其中的迭代算法即考虑了已知信息的卡尔曼滤波的迭代算法。数据检测算法是对信道矩阵进行QR分解，修正数据的误差以消除ICI。但是，由于ICI会影响频率域估计的准确性，在快速移动的环境下，需要进行大量的迭代，信道估计的结果在信噪比（SNR）较低的情况下也不准确。 

发明内容

1.快变OFDM系统的渐进迭代时变信道估计和ICI消除方法，其特征在于，在OFDM系统的信道估计中，将ICI与噪声之和SIN作为Kalman滤波器的去噪对象，从仅使用导频点信息开 始，渐进地增加用于迭代计算的信息，依以下步骤实现： 

步骤(1)，发射端产生发送数据，将导频数据按照梳状导频方式插入到发送数据中： 

发送端设定如下：s表示第s个OFDM符号，s＝1,2,…,s,…S，每个OFDM符号包含N个子载波，n＝1,2,…,n,…,N，其中包含N_p个导频符号和N_d个数据符号，N_d+N_p＝N，n_p＝1,2,…,N_p，导频在频域上的位置矩阵表示成：其中 且保证N_p≥L，L是信道多径数l的最大值，即l＝1，2，…,,l…,L，N_p个导频被平均的插入到N个载波之中并且在传输过程中保持不变，在N个载波中导频点符号表示为${x_p}^{\left(s\right)}=x^{\left(s\right)}\left(P_s\right)={[{x_{p_1}}^{\left(s\right)},{x_{p_2}}^{\left(s\right)},···,{x_{{p_N}_p}}^{\left(s\right)}]}^T,$

步骤(2)，数据通过OFDM系统发送至接收端，在接收端去掉循环前缀后，按以下步骤用多项式基扩展模型P-BEM对信道进行建模： 

步骤(2.1)，利用多项式基扩展模型P-BEM来描述具有时频双选特性的多径传播信道，则第s个OFDM符号的第n个子载波第l径的信道冲激响应h^(s)(n,l)表示为： 

h^(s)(n,l)＝QC_l^(s)+ξ_l^(s)(n),0≤n≤N-1， 

其中，ξ_l^(s)表示建模时每一个OFDM符号的第l径的模型误差，其值小于10^-3，在计算时忽略，即认为h^(s)(n,l)＝QC_l^(s)，Q是一个N×B的正交基函数矩阵，C_l^(s)则是由基函数对应的B个系数组成的向量${C_l}^{\left(s\right)}={[c_{1,l}^{\left(s\right)},c_{2,l}^{\left(s\right)},···,c_{b,l}^{\left(s\right)},···,c_{B,l}^{\left(s\right)}]}^T,$f_max是信道的最高频率，T_s是采样时间， 

步骤(2.2)，将在接收端的接收信号表示成以下形式： 

y^(s)＝H^(s)x^(s)+W^(s)， 

其中，x^(s)＝[x₁^(s),x₂^(s)…x_N^(s)]^T、y^(s)＝[y₁^(s),y₂^(s),…,y_N^(s)]^T分别表示频域上第s个去掉循环前缀后的发送信号和接收信号，W^(s)是其频域上的白噪声，H^(s)是N×N的信道矩阵： 

其中，矩阵的每一个元素为多径信道的信道冲击响应的和，计算方式如下： 

$H^{\left(s\right)}(m,k)=\underset{l=0}{\overset{L-1}{\Sigma }}{G_l}^{\left(s\right)}(M,K)e^{-j2\pi (\frac{k-1}{N}-\frac{1}{2}){\tau }_l},$

m，k表示上述矩阵H^(s)的第m行k列的值，m＝1,2,…,m,…,M，M≤N，k＝1,2,…,k,…,K，K≤N，τ_l是第l径的时延，G_l^(s)(M,K)为信道冲击相应的频域表达矩阵，其每一个元素计算如下： 

${G_l}^{\left(s\right)}(m,k)=\frac{1}{N}\underset{n=0}{\overset{N-1}{\Sigma }}{h}^{\left(s\right)}(n,l)e^{-j2\pi (m-k)n/N},$

步骤(2.3)，根据P-BEM模型将接收信号进行重新建模，表示成带有P-BEM系数的表达式如下： 

y^(s)＝Φ^(s)g^(s)+W^(s)， 

其中， 

g^(s)=[C₁^(s)T,C₂^(s)T…C_L^(s)T]^T，表示P-BEM算法中的系数矩阵， 

表示重新建模后，与发送数据相关的系数矩阵，其计算方法如下： 

$Z_l^{\left(s\right)}=\frac{1}{N}[D_1\mathrm{diag}\left(x^{\left(s\right)}\right){\Gamma }_l,···,D_b\mathrm{diag}\left(x^{\left(s\right)}\right){\Gamma }_l,···,D_B\mathrm{diag}\left(x^{\left(s\right)}\right){\Gamma }_l],$其中， 

${\Gamma }_l={\left(\begin{array}{cccc}{e}^{-j2\pi (\frac{p_1}{N}-\frac{1}{2}){\tau }_l}& e^{-j2\pi (\frac{p_2-1}{N}-\frac{1}{2}){\tau }_l}& ···& e^{-j2\pi (\frac{{p_N}_p-1}{N}-\frac{1}{2}){\tau }_l}\end{array}\right)}^T,$是第l径的傅里叶变换， 

Γ＝[Γ₁,Γ₂,…,Γ_L]，表示L径的总傅里叶变换矩阵， 

diag(x^(s))表示以向量x^(s)为对角元素的矩阵， 

步骤(3)，利用AR模型对信道BEM系数进行建模： 

步骤(3.1)，按下式计算C_l^(s)相关矩阵： 

$R_{C_l}^{\left(j\right)}={\left({\left(Q\right)}^{H}Q\right)}^{-1}{\left(Q\right)}^H{R}_{h(n,l)}^{\left(j\right)}Q{\left({\left(Q\right)}^{H}Q\right)}^{-1},$

其中，j表示相关的阶数，即进行相关运算的OFDM符号的符号间隔，j的取值为[-1,0,1]，分别表示当前OFDM信号的C_l^(s)和前一个符号的C_l^(s-1)的相关矩阵，当前OFDM信号的C_l^(s)的自相关矩阵，当前OFDM信号的C_l^(s)和后一个符号的C_l^(s+1)的相关矩阵。(·)^H表示 Hermitian运算，$R_{h(n,l)}^{\left(j\right)}=E\left[h(n,l)h^{*}(n+j,l)\right]={\sigma }_{h(n,l)}^2{J}_0\left(2\pi f_d{T}_{s}j\right),$其中E[·]表示均值,J₀(·)表示第一类的零阶贝塞尔函数，f_d=vf_c/c是终端的移动速度为v时的最大多普勒频移，f_c是载波频率，c是光速，代表第l径的信道冲击响应的方差，并假设

步骤(3.2)，根据Yule-Walker方程得到信道P-BEM参数的状态转移方程： 

g_(s)＝Ag_(s-1)+U_(s)， 

将OFDM系统发送符号的时间顺序g^(s)看做控制系统中状态转移过程g_(s)，即g^(s)＝g_(s)，状态转移方程系数A＝diag(a₁,a₂,…,a_l,…a_L)，diag(x)表示以向量x为对角元素的矩阵，U_(s)代表第s个OFDM符号的AR模型的建模误差； 

步骤(4)，对Kalman滤波器进行初始化并计算初始更新方程： 

步骤(4.1)，按照下式对Kalman滤波器进行初始化： 

$g_{\left(0\right|0)}=0_{\mathrm{LB},1},$$P_{\left(0\right|0)}=\mathrm{diag}(R_{C_1}^{\left(0\right)},R_{C_2}^{\left(0\right)}···R_{C_L}^{\left(0\right)}),$

形式如和P_(s|s)中下标前一个s均表示当前的状态为g_(s)，后一个s表示第s个OFDM符号， P_(0|0)是用于计算的初始值，表示OFDM符号的g_(s)的初始值，P_(0|0)表示对应的误差相关矩阵，O_LB,1是LB×1的零矩阵， 

步骤(4.2)，按下式计算Kalman的初始时间更新方程： 

i＝1，s＝1， 

${\hat{g}}_{\left(s\right)}=A{\hat{g}}_{\left(0\right|0)},$

P_(s)＝AP_(0|0)(A)^H+V[U_(s)]， 

i表示迭代次数，表示Kalman方程中状态估计g_(s)的中间变量，P_(s)表示中间变量对应的误差相关矩阵；用V[·]表示协方差矩阵，V[U_(s)]=diag(u₁,u₂…u_L)，

步骤(5)，进行第一次信道估计迭代运算，此时迭代次数i＝1，本次迭代中仅使用导频点所在的子载波处接收到数据做信道估计，除了导频点所在子载波外其他子载波上的数据视为ICI，用SIN方法消除未知数据对导频处信道估计的影响，实现无ICI干扰的导频辅助Kalman信道估计，具体步骤如下： 

步骤(5.1)，仅将接收信号中各个导频点对应的载波位置的接收信号用于计算，并将接收信号分成各个导频点所在子载波上的数据，除导频点子载波外的其它子载波上的数据对各个导频点所在子载波的干扰和噪声三部分，如下式所示： 

${y_p}^{\left(s\right)}=y^{\left(s\right)}\left(P_s\right)={H^{\left(s\right)}}_{[P_s,P_s]}{x_{p_{n_p}}}^{\left(s\right)}+{H^{\left(s\right)}}_{[P_s,d_{n_p{n}^{\prime }}]}{x_{d_{n_p{n}^{\prime }}}}^{\left(s\right)}+W^{\left(s\right)}\left(P_s\right)$

其中，$P_s=[p_1,p_2,···,p_{N_p}],$${x_{d_{n_p{n}^{\prime }}}}^{\left(s\right)}={[x^{\left(s\right)}\left(d_{11}\right),x^{\left(s\right)}\left(d_{12}\right),···{,x}^{\left(s\right)}\left(d_{{1N}^{\prime }}\right),x^{\left(s\right)}\left(d_{21}\right),···,x^{\left(s\right)}\left(d_{N_p{N}^{\prime }}\right)]}^T,$n'＝1,2,…,N'表示相邻的各个导频点之间的距离，是一个N_p×N_p的单位矩阵，σ²高斯白噪声W^(s)的方差，上式中第二项为除导频点子载波外的其它子载波上的数据对导频点所在子载波的干扰ICI， 

步骤(5.2)，将数据ICI干扰项考虑成信道噪声W^(s)(P_s)的一部分作为滤波器的去噪对象，将步骤（2）中的算法按照SIN估计的方法改写，令则SIN估计的Kalman观测方程表示为： 

$y_{p\left(s\right)}={{\Phi }^{\mathrm{SIN}}}_{\left(s\right)}{g}_{\left(s\right)}+W_{\left(s\right)}^{\mathrm{SIN}},$

其中： 

${{\Phi }^{\mathrm{SIN}}}_{\left(s\right)}=\frac{1}{N}[{z_1^{\left(s\right)}}_{\mathrm{SIN}},{Z_2^{\left(s\right)}}_{\mathrm{SIN}}···{Z_L^{\left(s\right)}}_{\mathrm{SIN}}],$

${Z_l^{\left(s\right)}}_{\mathrm{SIN}}=\frac{1}{N}[D_1^{\mathrm{SIN}}\mathrm{diag}\left({x_p}^{\left(s\right)}\right){\Gamma }_l^{\mathrm{SIN}}···D_B^{\mathrm{SIN}}\mathrm{diag}\left({x_p}^{\left(s\right)}{\Gamma }_l^{\mathrm{SIN}}\right)],$

${\Gamma }_l^{\mathrm{SIN}}={\left(\begin{array}{cccc}{e}^{-j2\pi (\frac{p_1}{N}-\frac{1}{2}){\tau }_1}& e^{-j2\pi (\frac{p_2-1}{N}-\frac{1}{2}){\tau }_l}& ···& e^{-j2\pi (\frac{{p_N}_p-1}{N}-\frac{1}{2}){\tau }_l}\end{array}\right)}^T,$

${\Gamma }^{\mathrm{SIN}}=\left(\begin{array}{cccc}{{\Gamma }_1}^{\mathrm{SIN}},& {{\Gamma }_2,}^{\mathrm{SIN}}& ···,& {{\Gamma }_l}^{\mathrm{SIN}}\end{array}\right)$

步骤(5.3)，计算的协方差矩阵

计算中假设ICI为高斯白噪声，令因为噪声与ICI二者相互独立，所以$V\left[W_{\left(s\right)}^{\mathrm{SIN}}\right]=U_{\mathrm{ICI}}+V\left[W^{\left(s\right)}\left(P_s\right)\right];$

U_ICI矩阵中的每一个元素的计算式为： 

$U_{\mathrm{ICI}}(m,k)=R_{\mathrm{ICI}}^{\left(j\right)}\approx 4{\pi }^2{T_s}^2{E}_s\left(\underset{l=0}{\overset{L-1}{\Sigma }}{\sigma }_{h(n,l)}^2{\sigma }_{D_l}\right)\rho (\alpha ,\mathrm{rag},N),$

其中，m，k表示矩阵的第m行k列，E_s是发送数据的功率，是功率为P_v时的多普勒功率普函数，f是传输频率，α表示用于迭代的最边缘数据距对应的各个导频点的距离，第一次计算时为0，rag为计算的精度，rag=[0,1,2,3]，并且: 

ρ(α,rag，N)=ρ(0,rag,N)-ρ₁(α,rag,N)， 

步骤(5.4)，分别按下述三式计算kalman增益K_(s)，第s个OFDM符号转移到状态的状态估计矩阵和与对应的协防差矩阵P_(s|s)，构成观测更新方程组，其中，Φ=Φ^SIN,$W_{\left(s\right)}=W_{\left(s\right)}^{\mathrm{SIN}},$

K_(s)＝P_(s)(Φ_(s))^H(Φ_(s)P_(s)(Φ_(s))^H+V[W_(s)])^-1， 

${\hat{g}}_{\left(s\right|s)}={\hat{g}}_{\left(s\right)}+K_{\left(s\right)}(y_{\left(s\right)}-{\Phi }_{\left(s\right)}{\hat{g}}_{\left(s\right)}),$

P_(s|s)＝P_(s)-K_(s)Φ_(s)P_(s)， 

步骤(5.5)，根据下式计算出信道矩阵H^(s)的估计值： 

$H^{\left(s\right)}=\frac{1}{N}\underset{b=1}{\overset{B}{\Sigma }}{D}_b\mathrm{diag}\left(\Gamma {\hat{g}}_{\left(s\right|s)}\right),$

步骤(5.6)，利用下式对信道矩阵进行QR分解，得到矩阵R^(s)： 

H^(s)＝IR^(s)

其中I是一个单位矩阵，R^(s)是一个上三角矩阵， 

步骤(5.7)，按下式对数据进行QR数据检测： 

其中y′^(s)＝(I)^Hy^(s)，和分别是数据的检测值和检测值星座图量化后的结果，[·]_m，k代表矩阵的第m行k列，[·]_m是向量的第m个元素，[·]_k是向量的第k个元素，O(·)表示解调运算，m，k表示矩阵H^(s)的第m行k列的值，m＝1,2,…,m,…，M，M≤N，k＝1,2,…,k,…,K，K≤N， 

步骤(6)，迭代次数i=i+1，迭代计算次数i>1时，第二次迭代用导频点所在处的子载波和其两侧相邻的一个子载波处接收到的数据做信道估计，其余子载波上的数据视为ICI，其后迭代每次增加的用于计算的数据，均为上次迭代用于计算的数据所在的子载波两侧的子载波处接收到的数据，其余子载波上的数据视为ICI，SIN方法进行信道估计迭代运算如下所述： 

步骤(6.1)，将接收信号分成导频点所在子载波数据，用于计算的数据所在的子载波数据，未用于计算的数据对导频点和用于计算的数据所在子载波的干扰和噪声四部分，如下式所示： 

上式中第三项为更新后的ICI干扰， 

步骤(6.2)，步骤(5.2)中所述的方法计算Φ^SIN_(s)， 

步骤(6.3)，随迭代次数增加，用于迭代的数据增加，此时更新的的协方差矩阵U_ICI计算如下： 

$U_{\mathrm{ICI}}\approx 4{\pi }^2{T_s}^2{E}_s\left(\underset{l=0}{\overset{L-1}{\Sigma }}{\sigma }_{h(n,l)}^2{\sigma }_{D_l}\right)[\frac{1}{2}\rho (A+m,\mathrm{rag},N)+\frac{1}{2}\rho (A-m,\mathrm{rag},N)]$

其中，A+m表示导频点右侧的用于计算的数据到对应的导频点的距离，A-m表示导频点左侧的用于计算的数据到对应的导频点的距离， 

步骤(6.4)，按步骤（5.5）所述方法计算kalman增益K_(s)，第s个OFDM符号转移到状态的状态估计矩阵和与对应的协防差矩阵Ｐ_(s|s)，构成观测更新方程组，其中，Φ=Φ^SIN，

K_(s)＝P_(s)(Φ_(s))^H(Φ_(s)P_(s)(Φ_(s))^H+V[W_(s)])^-1， 

${\hat{g}}_{\left(s\right|s)}={\hat{g}}_{\left(s\right)}+K_{\left(s\right)}(y_{\left(s\right)}-{\Phi }_{\left(s\right)}{\hat{g}}_{\left(s\right)}),$

P_(s|s)＝P_(s)-K_(s)Φ_(s)P_(s)， 

步骤(6.5)，按步骤（5.5）所述方法计算出信道矩阵的估计值H^(s)： 

$H^{\left(s\right)}=\frac{1}{N}\underset{b=1}{\overset{B}{\Sigma }}{D}_b\mathrm{diag}\left(\Gamma {\hat{g}}_{\left(s\right|s)}\right),$

步骤(6.6)，按步骤（5.6）所述方法对信道矩阵进行QR分解得到R^(s)： 

H^(s)＝IR^(s)， 

步骤(6.7)，按步骤（5.7）所述方法对数据进行QR数据检测： 

步骤(7)，判断是否所有的数据是否已经都用于迭代，如果是，则算法结束，如果不是，则继续， 

步骤(8)，通过判断迭代的次数，决定是否需要增加迭代算法的输入数据，判断算法如下： 

设定步长Δ，比较迭代次数i和Δ·μ+2，其中，$\mu =\mathrm{1,2},···,\Delta ·\mu +2<\frac{N}{N_p},$如果i＝Δ·μ+2，则选择已用于计算的数据所在的子载波两侧的子载波处的数据作为下次迭代新加入的用于计算的数据，与已用于计算的数据一起带入SIN算法，带入到步骤(6)中重新计； 

如果i≠Δ·μ+2，则不增加用于计算的数据，返回步骤(6)进行迭代； 

结束。 

附图说明

图1为本发明适用的OFDM系统； 

图2为本发明的基本原理流程图； 

图3为本发明中涉及的渐进步长示意图。其中图3（a）表示第一次计算时的接收信号，其中①即黑色的部分表示到频点所在子载波③表示ICI所在的子载波${x_{d_{n_p{n}^{\prime }}}}^{\left(s\right)}={[x^{\left(s\right)}\left(d_{11}\right),x^{\left(s\right)}\left(d_{12}\right),···,x^{\left(s\right)}\left(d_{{1N}^{\prime }}\right),x^{\left(s\right)}\left(d_{21}\right),···,x^{\left(s\right)}\left(d_{N_p{N}^{\prime }}\right)]}^T;$图3（b） 表示第二次计算时的接收信号，其中①表示导频点所在子载波②表示本次迭代新加入计算的数据所在的子载波和③表示表示ICI所在的子载波，④表示α，即用于迭代的最边缘的数据与导频点之间的距离，此时α＝1，⑤表示导频点之间的距离图3（c）表示迭代次数i＝Δ·μ+2，且i＞2时的接收信号，其中①表示上次用于计算的数据所在子载波 

${x^{\left(s\right)}\left(p_{N_p}\right)+[x^{\left(s\right)}\left(d_{11}\right),x^{\left(s\right)}\left(d_{12}\right),···,x^{\left(s\right)}\left(d_{1(i-1)}\right),···x^{\left(s\right)}\left(d_{n_{p}1}\right),···,x^{\left(s\right)}\left(d_{n_{p(i-1)}}\right),···,x^{\left(s\right)}\left(d_{N_p(i-1)}\right)]}^T,$②表示本次迭代新加入计算的数据所在的子载波和④表示α，此时，α＝Δ+1； 

图4为本发明与基于Kalman但未处理ICI的信道估计算法的性能对比。其中

和分别表示传统方法中1次、3次和10次迭代后的结果，

和表示本发明传中1次、3次和10次迭代后的结果，表示数据全部已知时， 

此种算法的理论值上限。 

具体实施方式

快变OFDM系统的渐进迭代时变信道估计和ICI消除方法，其特征在于，在OFDM系统的信道估计中，将ICI与噪声之和SIN作为Kalman滤波器的去噪对象，从仅使用导频点信息开始，渐进地增加用于迭代计算的信息，依以下步骤实现： 

步骤(1)，发射端产生发送数据，将导频数据按照梳状导频方式插入到发送数据中： 

发送端设定如下：s表示第s个OFDM符号，s＝1,2,…,s,…S，每个OFDM符号包含N个子载波，n＝1,2,…,n,…,N，其中包含N_p个导频符号和N_d个数据符号，N_d+N_p＝N，n_p＝1,2,…,N_p，导频在频域上的位置矩阵表示成：其中 且保证N_p≥L，L是信道多径数l的最大值，即l＝1,2,…,,l…,L，N_p个导频被平均的插入到N个载波之中并且在传输过程中保持不变，在N个载波中导频点符号表示为${x_p}^{\left(s\right)}=x^{\left(s\right)}\left(P_s\right)={[{x_{p_1}}^{\left(s\right)},{x_{p_2}}^{\left(s\right)},···,{x_{{p_N}_p}}^{\left(s\right)}]}^T,$

步骤(2)，数据通过OFDM系统发送至接收端，在接收端去掉循环前缀后，按以下步骤用多项式基扩展模型P-BEM对信道进行建模： 

步骤(2.1)，利用多项式基扩展模型P-BEM来描述具有时频双选特性的多径传播信道，则 第S个OFDM符号的第n个子载波第l径的信道冲激响应h^(s)(n,l)表示为： 

h^(s)(n,l)＝QC_l^(s)+ξ_l^(s)(n),0≤n≤N-1， 

其中，ξ₁^(s)表示建模时每一个OFDM符号的第l径的模型误差，其值小于10^-3，在计算时忽略，即认为h^(s)(n,l)＝QC_l^(s)，Q是一个N×B的正交基函数矩阵，C_l^(s)则是由基函数对应的B个系数组成的向量${{C_l}^{\left(s\right)}=[c_{1,l}^{\left(s\right)},c_{2,l}^{\left(s\right)},···,c_{b,l}^{\left(s\right)},···,c_{B,l}^{\left(s\right)}]}^T,$f_max是信道的最高频率，T_s是采样时间， 

步骤(2.2)，将在接收端的接收信号表示成以下形式： 

y^(s)＝H^(s)x^(s)+W^(s)， 

其中，x^(s)＝[x₁^(s),x₂^(s)…x_N^(s)]^T、y^(s)＝[y₁^(s),y₂^(s),…,y_N^(s)]^T分别表示频域上第s个去掉循环前缀后的发送信号和接收信号，W^(s)是其频域上的白噪声，H^(s)是N×N的信道矩阵： 

其中，矩阵的每一个元素为多径信道的信道冲击响应的和，计算方式如下： 

$H^{\left(s\right)}(m,k)=\underset{l=0}{\overset{L-1}{\Sigma }}{G_l}^{\left(s\right)}(M,K)e^{-j2\pi (\frac{k-1}{N}-\frac{1}{2}){\tau }_l},$

m，k表示上述矩阵H^(s)的第m行k列的值，m＝1,2,…,m,…,M，M≤N，k＝1,2,…,k,…,K，K≤N，τ_l是第l径的时延，G_l^(s)(M,K)为信道冲击相应的频域表达矩阵，其每一个元素计算如下： 

$G_l^{\left(s\right)}(m,k)=\frac{1}{N}\underset{n=0}{\overset{N-1}{\Sigma }}{h}^{\left(s\right)}(n,l)e^{-j2\pi (m-k)n/N},$

步骤(2.3)，根据P-BEM模型将接收信号进行重新建模，表示成带有P-BEM系数的表达式如下： 

y^(s)＝Φ^(s)g^(s)+W^(s)， 

其中， 

g^(s)=[C₁^(s)T,C₂^(s)T…C_L^(s)T]^T，表示PBEM算法中的系数矩阵， 

表示重新建模后，与发送数据相关的系数矩阵，其计算方法如下： 

$Z_l^{\left(s\right)}=\frac{1}{N}[D_1\mathrm{diag}\left(x^{\left(s\right)}\right){\Gamma }_l,···,D_b\mathrm{diag}\left(x^{\left(s\right)}\right){\Gamma }_l,···,D_B\mathrm{diag}\left(x^{\left(s\right)}\right){\Gamma }_l],$其中， 

${\Gamma }_l={\left(\begin{array}{cccc}{e}^{-j2\pi (\frac{p_1}{N}-\frac{1}{2}){\tau }_l}& e^{-j2\pi (\frac{p_2-1}{N}-\frac{1}{2}){\tau }_l}& ···& e^{-j2\pi (\frac{{p_N}_p-1}{N}-\frac{1}{2}){\tau }_l}\end{array}\right)}^T,$是第l径的傅里叶变换， 

Γ＝[Γ₁,Γ₂,…,Γ_L]，表示L径的总傅里叶变换矩阵， 

diag(x^(s))表示以向量x^(s)为对角元素的矩阵， 

步骤(3)，利用AR模型对信道BEM系数进行建模： 

步骤(3.1)，按下式计算C_l^(s)相关矩阵： 

$R_{C_l}^{\left(j\right)}={\left({\left(Q\right)}^{H}Q\right)}^{-1}{\left(Q\right)}^H{R}_{h(n,l)}^{\left(j\right)}Q{\left({\left(Q\right)}^{H}Q\right)}^{-1},$

其中，j表示相关的阶数，即进行相关运算的OFDM符号的符号间隔，j的取值为[-1,0,1]，分别表示当前OFDM信号的C_l^(s)和前一个符号的C_l^(s-1)的相关矩阵，当前OFDM信号的C_l^(s)的自相关矩阵，当前OFDM信号的C_l^(s)和后一个符号的C_l^(s+1)的相关矩阵。(·)^H表示Hermitian运算，$R_{h(n,l)}^{\left(j\right)}=E\left[h(n,l)h^{*}(n+j,l)\right]={\sigma }_{h(n,l)}^2{J}_0\left(2\pi f_d{T}_{s}j\right),$其中E[·]表示均值,J₀(·)表示第一类的零阶贝塞尔函数，f_d=vf_c/c是终端的移动速度为v时的最大多普勒频移，f_c是载波频率，c是光速，代表第l径的信道冲击响应的方差，并假设

步骤(3.2)，根据YuleWalker方程得到信道P-BEM参数的状态转移方程： 

g_(s)＝Ag_(s-1)+U_(s)， 

将OFDM系统发送符号的时间顺序g^(s)看做控制系统中状态转移过程g_(s)，即g^(s)＝g_(s)，状态转移方程系数A＝diag(a₁,a₂,…,a_l,…a_L)，diag(x)表示以向量x为对角元素的矩阵，U_(s)代表第s个OFDM符号的AR模型的建模误差； 

步骤(4)，对Kalman滤波器进行初始化并计算初始更新方程： 

步骤(4.1)，按照下式对Kalman滤波器进行初始化： 

${\hat{g}}_{\left(0\right|0)}=0_{LB,1},P_{\left(0\right|0)}=\mathrm{diag}(R_{C_1}^{\left(0\right)},R_{C_2}^{\left(0\right)}···R_{C_L}^{\left(0\right)}),$

形式如和P_(s|s)中下标前一个s均表示当前的状态为g_(s)，后一个s表示第s个OFDM符号， P_(0|0)是用于计算的初始值，表示OFDM符号的g_(s)的初始值，P_(0|0)表示对应的误差相关矩阵，O_LB,1是LB×1的零矩阵， 

步骤(4.2)，按下式计算Kalman的初始时间更新方程： 

i＝1，s＝1， 

${\hat{g}}_{\left(s\right)}=A{\hat{g}}_{\left(0\right|0)},$

P_(s)＝AP_(0|0)(A)^H+V[U_(s)]， 

i表示迭代次数，表示Kalman方程中状态估计g_(s)的中间变量，P_(s)表示中间变量对应的误差相关矩阵；用V[·]表示协方差矩阵，V[U_(s)]=diag(u₁,u₂…u_L)，

步骤(5)，进行第一次信道估计迭代运算，此时迭代次数i＝1，本次迭代中仅使用导频点所在的子载波处接收到数据做信道估计，除了导频点所在子载波外其他子载波上的数据视为ICI，用SIN方法消除未知数据对导频处信道估计的影响，实现无ICI干扰的导频辅助Kalman信道估计，具体步骤如下： 

步骤(5.1)，仅将接收信号中各个导频点对应的载波位置的接收信号用于计算，并将接收信号分成各个导频点所在子载波上的数据，除导频点子载波外的其它子载波上的数据对各个导频点所在子载波的干扰和噪声三部分，如下式所示： 

${y_p}^{\left(s\right)}=y^{\left(s\right)}\left(P_s\right)={H^{\left(s\right)}}_{[P_s,P_s]}{x_{{p_n}_p}}^{\left(s\right)}+{H^{\left(s\right)}}_{[P_s,d_{n_p{n}^{\prime }]}}{x_{d_{n_p{n}^{\prime }}}}^{\left(s\right)}+W^{\left(s\right)}\left(P_s\right)$

其中，$P_s=[p_1,p_2,···,{p_N}_p],$${x_{d_{n_p{n}^{\prime }}}}^{\left(s\right)}={[x^{\left(s\right)}\left(d_{11}\right),x^{\left(s\right)}\left(d_{12}\right),···{,x}^{\left(s\right)}\left(d_{{1N}^{\prime }}\right),x^{\left(s\right)}\left(d_{21}\right),···,x^{\left(s\right)}\left(d_{N_p{N}^{\prime }}\right)]}^T,$n'＝1,2,…,N'表示相邻的各个导频点之间的距离，是一个N_p×N_p的单位矩阵，σ²高斯白噪声W^(s)的方差，上式中第二项为除导频点子载波外的其它子载波上的数据对导频点所在子载波的干扰ICI， 

步骤(5.2)，将数据ICI干扰项考虑成信道噪声W^(s)(P_s)的一部分作为滤波器的去噪对象，将步骤（2）中的算法按照SIN估计的方法改写，令则SIN估计的Kalman观测方程表示为： 

$y_{p\left(s\right)}={{\Phi }^{\mathrm{SIN}}}_{\left(s\right)}{g}_{\left(s\right)}+W_{\left(s\right)}^{\mathrm{SIN}},$

其中： 

${{\Phi }^{\mathrm{SIN}}}_{\left(s\right)}=\frac{1}{N}[{z_1^{\left(s\right)}}_{\mathrm{SIN}},{Z_2^{\left(s\right)}}_{\mathrm{SIN}}···{Z_L^{\left(s\right)}}_{\mathrm{SIN}}],$

${Z_l^{\left(s\right)}}_{\mathrm{SIN}}=\frac{1}{N}[D_1^{\mathrm{SIN}}\mathrm{diag}\left({x_p}^{\left(s\right)}\right){\Gamma }_l^{\mathrm{SIN}}···D_B^{\mathrm{SIN}}\mathrm{diag}\left({x_p}^{\left(s\right)}\right){\Gamma }_l^{\mathrm{SIN}}],$

${\Gamma }_l^{\mathrm{SIN}}={\left(\begin{array}{cccc}{e}^{-j2\pi (\frac{p_1}{N}-\frac{1}{2}){\tau }_l}& e^{-j2\pi (\frac{p_2-1}{N}-\frac{1}{2}){\tau }_l}& ···& e^{-j2\pi (\frac{{p_N}_p-1}{N}-\frac{1}{2}){\tau }_l}\end{array}\right)}^T,$

Γ^SIN＝[Γ₁^SIN,Γ₂,^SIN…,Γ_l^SIN] 

步骤(5.3)，计算的协方差矩阵

计算中假设ICI为高斯白噪声，令因为噪声与ICI二者相互独立，所以$V\left[W_{\left(s\right)}^{\mathrm{SIN}}\right]=U_{\mathrm{ICI}}+V\left[W^{\left(s\right)}\left(P_S\right)\right];$

U_ICI矩阵中的每一个元素的计算式为： 

$U_{\mathrm{ICI}}(m,k)=R_{\mathrm{ICI}}^{\left(j\right)}\approx 4{\pi }^2{T_s}^2{E}_s\left(\underset{l=0}{\overset{L-1}{\Sigma }}{\sigma }_{h(n,l)}^2{\sigma }_{D_l}\right)\rho (\alpha ,\mathrm{rag},N),$

其中，m，k表示矩阵的第m行k列，E_s是发送数据的功率，是功率为P_v时的多普勒功率普函数，f是传输频率，α表示用于迭代的最边缘数据距对应的各个导频点的距离，第一次计算时为0，rag为计算的精度，rag=[0,1,2,3]，并且: 

ρ(α,rag，N)=ρ(0,rag,N)-ρ₁(α,rag,N) 

步骤(5.4)，分别按下述三式计算kalman增益K_(s)，第s个OFDM符号转移到状态的状态估计矩阵和与对应的协防差矩阵P_(s|s)，构成观测更新方程组，其中，Φ=Φ^SIN$W_{\left(s\right)}=W_{\left(s\right)}^{\mathrm{SIN}},$

K_(s)＝P_(s)(Φ_(s))^H(Φ_(s)P_(s)(Φ_(s))^H+V[W_(s)])^-1， 

${\hat{g}}_{\left(s\right|s)}={\hat{g}}_{\left(s\right)}+K_{\left(s\right)}(y_{\left(s\right)}-{\Phi }_{\left(s\right)}{\hat{g}}_{\left(s\right)}),$

P_(s|s)＝P_(s)-K_(s)Φ_(s)P_(s)， 

步骤(5.5)，根据下式计算出信道矩阵H^(s)的估计值： 

$H^{\left(s\right)}=\frac{1}{N}\underset{b=1}{\overset{B}{\Sigma }}{D}_b\mathrm{diag}\left(\Gamma {\hat{g}}_{\left(s\right|s)}\right),$

步骤(5.6)，利用下式对信道矩阵进行QR分解，得到矩阵R^(s)： 

H^(s)＝IR^(s)

其中I是一个单位矩阵，R^(s)是一个上三角矩阵， 

步骤(5.7)，按下式对数据进行QR数据检测： 

其中y′^(s)＝(I)^Hy^(s)，和分别是数据的检测值和检测值星座图量化后的结果，[·]_m，k代表矩阵的第m行k列，[·]_m是向量的第m个元素，[·]_k是向量的第k个元素，O(·)表示解调运算，m，k表示矩阵H^(s)的第m行k列的值，m＝1,2,…,m,…，M，M≤N，k＝1,2,…,k,…,K，K≤N， 

步骤(6)，迭代次数i=i+1，迭代计算次数i>1时，第二次迭代用导频点所在处的子载波和其两侧相邻的一个子载波处接收到的数据做信道估计，其余子载波上的数据视为ICI，其后迭代每次增加的用于计算的数据，均为上次迭代用于计算的数据所在的子载波两侧的子载波处接收到的数据，其余子载波上的数据视为ICI，SIN方法进行信道估计迭代运算如下所述： 

步骤(6.1)，将接收信号分成导频点所在子载波数据，用于计算的数据所在的子载波数据，未用于计算的数据对导频点和用于计算的数据所在子载波的干扰和噪声四部分，如下式所示： 

上式中第三项为更新后的ICI干扰， 

步骤(6.2)，步骤(5.2)中所述的方法计算Φ^SIN_(s)， 

步骤(6.3)，随迭代次数增加，用于迭代的数据增加，此时更新的的协方差矩阵U_ICI计算如下： 

$U_{\mathrm{ICI}}\approx 4{\pi }^2{{T_s}^{2}E}_s\left(\underset{l=0}{\overset{L-1}{\Sigma }}{\sigma }_{h(n,l)}^2{\sigma }_{D_l}\right)[\frac{1}{2}\rho (A+m,\mathrm{rag},N)+\frac{1}{2}\rho (A-m,\mathrm{rag},N)]$

其中，A+m表示导频点右侧的用于计算的数据到对应的导频点的距离，A-m表示导频点左侧的用于计算的数据到对应的导频点的距离， 

步骤(6.4)，按步骤（5.5）所述方法计算kalman增益K_(s)，第s个OFDM符号转移到状态的状态估计矩阵和与对应的协防差矩阵P_(s|s)，构成观测更新方程组，其中，Φ=Φ^SIN，

K_(s)＝P_(s)(Φ_(s))^H(Φ_(s)P_(s)(Φ_(s))H+V[W_(s)])^-1， 

${\hat{g}}_{\left(s\right|s)}={\hat{g}}_{\left(s\right)}+K_{\left(s\right)}(y_{\left(s\right)}-{\Phi }_{\left(s\right)}{\hat{g}}_{\left(s\right)}),$

P_(s|s)＝P_(s)-K_(s)Φ_(s)P_(s)， 

步骤(6.5)，按步骤（5.5）所述方法计算出信道矩阵的估计值H^(s)： 

$H^{\left(s\right)}=\frac{1}{N}\underset{b=1}{\overset{B}{\Sigma }}{D}_b\mathrm{diag}\left(\Gamma {\hat{g}}_{\left(s\right|s)}\right),$

步骤(6.6)，按步骤（5.6）所述方法对信道矩阵进行QR分解得到R^(s)： 

H^(s)＝IR^(s)， 

步骤(6.7)，按步骤（5.7）所述方法对数据进行QR数据检测： 

步骤(7)，判断是否所有的数据是否已经都用于迭代，如果是，则算法结束，如果不是，则继续， 

步骤(8)，通过判断迭代的次数，决定是否需要增加迭代算法的输入数据，判断算法如下： 

设定步长Δ，比较迭代次数i和Δ·μ+2，其中，$\mu =\mathrm{1,2}···,\Delta ·\mu +2<\frac{N}{N_p},$如果i＝Δ·μ+2，则选择已用于计算的数据所在的子载波两侧的子载波处的数据作为下次迭代新加入的用于计算的数据，与已用于计算的数据一起带入SIN算法，带入到步骤(6)中重新计； 

如果i≠Δ·μ+2，则不增加用于计算的数据，返回步骤(6)进行迭代； 

结束。