本文主要来源于复现文献的部分内容，有一定的参考价值：
[1] 刘波. 基于EM的突发通信参数估计技术研究[D]. 2009.
文末有代码和参考文献网盘下载地址，有问题欢迎留言交流！

1 引言

对于多径传播，接收信号是由源信号与它的各次回波叠加而成的。对接收到的信号进行时间延迟估计广泛应用于雷达、声纳、电声测量、石油地震勘探和生物医学等领域，具有重要的实际应用意义。

对于雷达、声纳以及无线通信等应用而言，均需要对各径信号的时延进行估计，然后进行RAKE合并等处理，以避免由于多径效应带来的系统性能恶化。本文考虑发送信号已知情况下的多径时延估计问题，此类多径时延估计问题是主动式雷达与声纳系统中常见的；另外，无线通信领域中用于对信道进行估计和均衡的导频信号，非合作源定位所用的照射信号中具有某些已知时域特征的部分等也属此类问题。

经典的时间延迟估计技术采用的是匹配滤波器的方法，它利用接收信号与发射信号求相关后发生多峰的位置对时间延迟做出估计，但是这种方法的多径分辨能力依赖于发射信号相关函数的主峰宽度(近似与信号带宽成反比)。当信号带宽较窄时，基于相关运算的方法会由于各个峰值的相互重叠而带来较大的估计误差。故目前多径环境下时间延迟估计问题的重点是研究相关不可解情况下的高分辨率多径时间延迟估计算法。

最大似然估计是参数估计问题的有效方法，具有近似最佳的估计性能及稳健性。但是，鉴于多信号多参数估计问题所需多维优化的复杂性，因此在多径时延估计问题中直接使用最大似然估计方法并不可行。EM算法可以将多维优化的复杂问题分解为一系列一维优化问题的迭代，在减小运算复杂性的同时可以保证估计性能。

2 问题背景及信号模型

多径时延估计是现代通信信号处理中信号检测与参数估计问题的一个重要组成部分。发送端的发送信号由于受到传播空间中不同物体的反射和漫射，使得接收机收到的信号是多条路径上信号的合成，即接收信号r(t)r(t)r(t)可以简单表示为如下形式[1]：

r(t)=∑l=1Lhls(t−τl)+w(t)t∈[0,T)r(t) = \sum\limits_{l = 1}^L {{h_l}s(t - {\tau _l}) + w(t)} \quad t \in [0,T)r(t)=l=1∑Lhls(t−τl)+w(t)t∈[0,T) (1)

式中，s(t)s(t)s(t)为t时刻的发送信号；w(t)w(t)w(t)为t时刻的观测噪声，其符合高斯分布，平均功率为σ2{{\sigma }^{2}}σ2；参数τl{{\tau }_{l}}τl和hl{{h}_{l}}hl分别为各条多径信号的时延和复信道系数，包含有信号由发送至接收所经历的时间延迟、幅度衰减、以及相位变化等信息，反映了信号的传播距离、传播介质以及反射体的特性等；LLL为存在的总的多径个数，[0,T)[0,T)[0,T)为信号的观测区间。

3 基于EM的多径时延估计

由第2章可知，接收信号可以表示为(1)式，总的信号路径个数为L，因此带估计参数向量可以表示如下

θ=(τ1,τ2,⋯,τL,h1,h2,⋯,hL)T\theta ={{({{\tau }_{1}},{{\tau }_{2}},\cdots ,{{\tau }_{L}},{{h}_{1}},{{h}_{2}},\cdots ,{{h}_{L}})}^{T}}θ=(τ1,τ2,⋯,τL,h1,h2,⋯,hL)T (2)

本文考虑确定性信号下的[\theta ]估计问题。为不失一般性，假设

∫Ts2(t)dt=1\int_{T}{{{s}^{2}}(t)dt=1}∫Ts2(t)dt=1 (3)

则由(2)可得[r(t)]的对数似然函数为

ln⁡fR(r;θ)=C−12σ2∫[r(t)−∑l=1Lhls(t−τl)]2dt\ln {{f}_{R}}(r;\theta )=C-\frac{1}{2{{\sigma }^{2}}}{{\int{\left[ r(t)-\sum\limits_{l=1}^{L}{{{h}_{l}}s(t-{{\tau }_{l}})} \right]}}^{2}}dtlnfR(r;θ)=C−2σ21∫[r(t)−l=1∑Lhls(t−τl)]2dt (4)

其中，C为常数。为获得hl{{h}_{l}}hl和τl{{\tau }_{l}}τl的最大似然估计量，必须解决式(5)的最优化问题。

min⁡θ∫[r(t)−∑l=1Lhls(t−τl)]2\underset{\theta }{\mathop{\min }}\,{{\int{\left[ r(t)-\sum\limits_{l=1}^{L}{{{h}_{l}}s(t-{{\tau }_{l}})} \right]}}^{2}}θmin∫[r(t)−l=1∑Lhls(t−τl)]2 (5)

这是一个2N元的复杂优化问题。当然，使用暴力也可以解决问题，在粗略的网格上通过搜索目标函数求值，粗略定位全局最小值，然后应用高斯法、牛顿－拉夫森法或其他迭代梯度搜索算法进行求解。然而，这些方法往往是非常复杂的，且需要大量的计算时间。

根据文献[1]，对进行分解，以参与叠加的各个单径信号分量为基础构造完备数据如下所示

x(t)=(x1(t),x2(t),⋯,xL(t))Tx(t)={{\left( {{x}_{1}}(t),{{x}_{2}}(t),\cdots ,{{x}_{L}}(t) \right)}^{T}}x(t)=(x1(t),x2(t),⋯,xL(t))T (6)

其中

xl(t)=hls(t−τl)+wl(t){{x}_{l}}(t)={{h}_{l}}s(t-{{\tau }_{l}})+{{w}_{l}}(t)xl(t)=hls(t−τl)+wl(t)

∑l=1Lwl(t)=w(t)\sum\limits_{l=1}^{L}{{{w}_{l}}(t)}=w(t)l=1∑Lwl(t)=w(t)

r(t)=1Tx(t)1T=(1,1,⋯,1)r(t)={{\mathbf{1}}^{T}}x(t)\begin{matrix} {} & {} \\ \end{matrix}{{\mathbf{1}}^{T}}=(1,1,\cdots ,1)r(t)=1Tx(t)1T=(1,1,⋯,1)

则有

r(t)=∑l=1Lxl(t)r(t)=\sum\limits_{l=1}^{L}{{{x}_{l}}(t)}r(t)=l=1∑Lxl(t) (7)

我们希望通过使ln⁡fR(r;θ)\ln {{f}_{R}}(r;\theta )lnfR(r;θ)最大来求解的最大似然估计量（MLE），但这是困难的，以求ln⁡fX(x;θ)\ln {{f}_{X}}(x;\theta )lnfX(x;θ)的最大值来代替，因x无法求得，故使用对数似然函数的条件数学期望来代替对数似然函数[2]，即

E{ln⁡fX(x;θ)∣r;θ}=∫ln⁡fX(x;θ)f(x∣r;θ)dxE\left\{ \ln {{f}_{X}}(x;\theta )|r;\theta \right\}=\int{\ln }{{f}_{X}}(x;\theta )f(x|r;\theta )dxE{lnfX(x;θ)∣r;θ}=∫lnfX(x;θ)f(x∣r;θ)dx (8)

上式中，必须知道θ\thetaθ才能确定ln⁡fX(x;θ)\ln {{f}_{X}}(x;\theta )lnfX(x;θ)，故对数似然函数的期望将作为当前的猜测，令θi{{\theta }^{i}}θi表示θ\thetaθ的MLE的第iii次猜测，则根据

ln⁡fX(x;θ)=∑l=1Lln⁡fXl(xl;θ)=C1−∑l=1L12σl2∫[xl(t)−hls(t−τl)]2dt\ln {{f}_{X}}(x;\theta )=\sum\limits_{l=1}^{L}{\ln {{f}_{{{X}_{l}}}}({{x}_{l}};\theta )}={{C}_{1}}-\sum\limits_{l=1}^{L}{\frac{1}{2\sigma _{l}^{2}}}{{\int{\left[ {{x}_{l}}(t)-{{h}_{l}}s(t-{{\tau }_{l}}) \right]}}^{2}}dtlnfX(x;θ)=l=1∑LlnfXl(xl;θ)=C1−l=1∑L2σl21∫[xl(t)−hls(t−τl)]2dt (9)

式中，C1为与待估参数无关的常数。确定完备数据的平均对数似然函数为

Q(θ,θi)=E{ln⁡fX(x;θ)∣r;θ}=C2−∑l=1L12σl2∫[xli(t)−hls(t−τl)]2dtQ(\theta ,{{\theta }^{i}})=E\left\{ \ln {{f}_{X}}(x;\theta )|r;\theta \right\}={{C}_{2}}-\sum\limits_{l=1}^{L}{\frac{1}{2\sigma _{l}^{2}}}{{\int{\left[ x_{l}^{i}(t)-{{h}_{l}}s(t-{{\tau }_{l}}) \right]}}^{2}}dtQ(θ,θi)=E{lnfX(x;θ)∣r;θ}=C2−l=1∑L2σl21∫[xli(t)−hls(t−τl)]2dt (10)

式中，C2为与待估参数无关的常数。其中，xli(t)x_{l}^{i}(t)xli(t)由观测数据r(t)r(t)r(t)以及第iii次迭代中获得的待估参数θi{{\theta }^{i}}θi给出，利用联合高斯随机变量的条件期望的标准结果，可以得到如下表示

xli(t)=E{xl(t)/∑l=1Lxl(t)=r(t);θi}=E(xl)+CxrCrr−1(r−E(r))x_{l}^{i}(t)=E\left\{ {{{x}_{l}}(t)}/{\sum\limits_{l=1}^{L}{{{x}_{l}}(t)=r(t);}}\;{{\theta }^{i}} \right\}=E\left( {{x}_{l}} \right)+{{C}_{xr}}C_{rr}^{-1}(r-E(r))xli(t)=E{xl(t)/l=1∑Lxl(t)=r(t);θi}=E(xl)+CxrCrr−1(r−E(r)) (11)

由此可以看出，EM算法的M步骤中，对条件期望函数Q(θ,θi)Q(\theta ,{{\theta }^{i}})Q(θ,θi)的最大化问题可以通过对求和项中的各分量分别最大化来实现，即

θ^i+1=arg⁡max⁡θQ(θ,θ^i){{\hat{\theta }}^{i+1}}=\arg \underset{\theta }{\mathop{\max }}\,Q(\theta ,{{\hat{\theta }}^{i}})θ^i+1=argθmaxQ(θ,θ^i) (12)

则EM算法步骤可以描述为：

E-step：根据所获得的条件期望函数Q(θ,θi)Q(\theta ,{{\theta }^{i}})Q(θ,θi)的表达式，利用第iii次迭代得到的参数估计θ^i{{\hat{\theta }}^{i}}θ^i。对 l=1,2,…,L，构造数据

xli(t)=hlis(t−τli)+βl[r(t)−∑l=1Lhlis(t−τli)]x_{l}^{i}(t)=h_{l}^{i}s(t-\tau _{l}^{i})+{{\beta }_{l}}\left[ r(t)-\sum\limits_{l=1}^{L}{h_{l}^{i}s(t-\tau _{l}^{i})} \right]xli(t)=hlis(t−τli)+βl[r(t)−l=1∑Lhlis(t−τli)] (13)

M-step：根据上一步得到的数据，对l=1,2,…,L，计算使(14)式最优化的参数hli+1,τli+1h_{l}^{i+1},\tau _{l}^{i+1}hli+1,τli+1

min⁡hl,τl∫[xli(t)−hls(t−τl)]2dt\underset{{{h}_{l}},{{\tau }_{l}}}{\mathop{\min }}\,{{\int{\left[ x_{l}^{i}(t)-{{h}_{l}}s(t-{{\tau }_{l}}) \right]}}^{2}}dthl,τlmin∫[xli(t)−hls(t−τl)]2dt (14)

其中，

∑l=1Lβl=1\sum\limits_{l=1}^{L}{{{\beta }_{l}}=1}l=1∑Lβl=1

经过化简计算，得对应于同一信号的参数估计(h^l,τ^l)({{\hat{h}}_{l}},{{\hat{\tau }}_{l}})(h^l,τ^l)为

τ^l=arg⁡max⁡τl∫xli(t)s(t−τl)dt{{\hat{\tau }}_{l}}=\arg \underset{{{\tau }_{l}}}{\mathop{\max }}\,\int{x_{l}^{i}(t)s(t-{{\tau }_{l}})}dtτ^l=argτlmax∫xli(t)s(t−τl)dt (15)

h^l=∫xli(t)s(t−τli+1)dt∫s2(t−τli+1)dt{{\hat{h}}_{l}}=\frac{\int{x_{l}^{i}(t)s(t-\tau _{l}^{i+1})}dt}{\int{{{s}^{2}}(t-\tau _{l}^{i+1})dt}}h^l=∫s2(t−τli+1)dt∫xli(t)s(t−τli+1)dt (16)

重复以上两个步骤直到算法收敛。

4 实验结果及分析

4.1 仿真参数设置

为了验证上述所提出的对多径时延EM估计方法的有效性，进行蒙特卡罗仿真实验，实验次数N=1000。仿真中，采用加窗的线性调频信号s(t)s(t)s(t)作为已知的导频信号，如下所示

s(t)=g(t)ej2π(t−T2)2,0<t≤Ts(t)=g(t){{e}^{j2\pi {{\left( t-\frac{T}{2} \right)}^{2}}}},\ 0<t\le Ts(t)=g(t)ej2π(t−2T)2, 0<t≤T (17)

其中，g(t)g(t)g(t)为钟形的窗函数，如下所示：

g(t)={0.5−0.5cos⁡(πt/Tw)0<t<Tw1Tw⩽t⩽T−Tw,Tw=T100.5−0.5cos⁡[π(t−T)/Tw]T−Tw<t⩽Tg(t)=\left\{\begin{array}{ll}0.5-0.5 \cos \left(\pi t / T_{w}\right) & 0<t<T_{w} \\ 1 & T_{w} \leqslant t \leqslant T-T_{w}, T_{w}=\frac{T}{10} \\ 0.5-0.5 \cos \left[\pi(t-T) / T_{w}\right] & T-T_{w}<t \leqslant T\end{array}\right.g(t)=⎩⎨⎧0.5−0.5cos(πt/Tw)10.5−0.5cos[π(t−T)/Tw]0<t<TwTw⩽t⩽T−Tw,Tw=10TT−Tw<t⩽T

信号的调频斜率K=100000，持续时间T=20.6ms，采样速率fs{{f}_{s}}fs=10080，采样周期Ts=1/fs{{T}_{s}}=1/{{f}_{s}}Ts=1/fs，信号带宽B=2KT=4.12kHzB=2KT=4.12kHzB=2KT=4.12kHz。设接收信号r(t)r(t)r(t)由2路多径信号构成，其参数分别为

τl=16Ts,h1=ejπ/4;τ2=20Ts,h1=ejπ/8{{\tau }_{l}}=16{{T}_{s}},{{h}_{1}}={{e}^{j\pi /4}};{{\tau }_{2}}=20{{T}_{s}},{{h}_{1}}={{e}^{j\pi /8}}τl=16Ts,h1=ejπ/4;τ2=20Ts,h1=ejπ/8 (19)

在EM算法的迭代过程中，设置迭代收敛的判定条件为对数似然函数的增量小于其取值的0.1%，满足此条件时的参数取值即为最终估计值，即定义如下的收敛判定条件：

Λ(r;θi+1)−Λ(r;θi)≤Λ(r;θi)×0.1\Lambda \left( \mathbf{r};{{\theta }^{i+1}} \right)-\Lambda \left( \mathbf{r};{{\theta }^{i}} \right)\le \Lambda \left( \mathbf{r};{{\theta }^{i}} \right)\times 0.1%Λ(r;θi+1)−Λ(r;θi)≤Λ(r;θi)×0.1 (20)

4.2 仿真结果及分析

钟形函数g(t)g(t)g(t)的时域波形如图1所示，线性调频信号s(t)s(t)s(t)的实部、虚部的时域波形如图2，线性调频信号的时频图如图3所示。

图1 钟形窗函数g(t)时域波形图 图2 线性调频信号s(t)的实部、虚部波形图 图3 线性调频信号s(t)的时频图

图4～图６为对时延参数和传播系数参数估计的均方根误差性能随信噪比变化的关系曲线。

图４时延估计均方根误差随信噪比变化的关系曲线 图5 幅度衰减因子实部均方根误差随信噪比变化的关系曲线 图6 幅度衰减因子虚部均方根误差随信噪比变化的关系曲线

从图中所示的仿真结果可以看出，在信噪比较低的情况下，会存在多径时延估计不准的偶然情况，除信噪比SNR=2.5dB时，时延估计存在误差外，其他信噪比下误差均为0。此外，在信噪比SNR≥5.5dB时，幅度衰减因子的实部、虚部均方根误差可达到10-3左右。仿真结果表明，高斯白噪声下基于EM的多径时延估计算法能够较好地估计多径时延和幅度衰减因子。

5 总结与展望

当多径信号时延不是采样周期的整数倍时，要提高时延估计的精度需要通过内插处理才能得到。在EM算法的每次迭代中都使用内插将会大大增加算法的运算量，因此可以首先对接收数据进行预处理，通过DFT将其变换到频率域，再进行EM算法迭代处理，即可避免了原先在迭代过程中频繁的内插处理。迭代初始值的设置既影响EM算法的收敛结果，同时也影响算法的收敛速度。初值的选择和迭代更新可参考轮换投影（AP）方法。

6 参考文献

[1] 刘波. 基于EM的突发通信参数估计技术研究[D]. 2009.
[2] Steven M. Kay. 统计信号处理基础——检测与估计理论（卷I、卷II合集）[M], 2014. 154-156.
[3] Meir Feder, Ehud Weinstein. Multipath Time-Delay Estimation via the EM Algorithm[D]. Woods Hole Oceanographic Institution, 1987.

7 代码

下载:https://download.csdn.net/download/wlwdecs_dn/12620599

8 获取更多内容请关注：

高斯白噪声下基于EM的多径时延估计相关推荐

仿射投影算法 matlab,基于仿射投影算法的时延估计
Signal Process & System 号与系基于仿射投影算法的时延估计摘要:在时间延迟计中,通常利用互相关箅法对8寸延进行估计.然而,互相关算法受操声影较人,在低信比时无法准确对 ...
matlab删失数据威布尔,随机删失数据下基于EM算法的Weibull分布参数估计
1引言生存数据分析的统计方法在生物医学和可靠性研究中,甚至在犯罪学.社会学.市场学以及健康保险领域的研究中,有重要的应用.但是,在生存数据的观测中常常由于动物的意外死亡.被调查者拒绝回答有关调查项目. ...
基于线性调频（LFM）信号分数阶傅里叶变换的多径时延估计
仿真信号仿真信号为线性调频信号,信号采样率为2kHz. 线性调频(LFM)信号参数为:初始频率为300Hz,带宽400Hz,脉宽1s,调频斜率为 k ^ = − F s Δ t c o t ( p ...
如何在信号中添加指定信噪比的高斯白噪声，为何深度学习去噪研究采用高斯白噪声？
在信号或者图像的降噪研究中,很多学者采用高斯白噪声添加到干净的样本中,来模拟含有噪声的样本,并以此来验证提出模型的降噪效果(比如降噪自编码器--Denoising Autoencoder).有一次投稿 ...
Matlab高斯白噪声通过互相关函数求功率谱密度
参考文章:百度文库文章但是该文章有很多部分,不知是matlab更新了还是咋地,有错误,踩了很多坑,仅做参考. 产生高斯白噪声的方法在"信号处理"专栏中有写,其时域和频域图如下: ...
基于EM算法的高斯混合模型参数估计
一算法思想对于样本数据集所含变量都是可观测的,我们一般应用极大似然估计法或者贝叶斯估计法估计模型参数.但若样本数据集中存在不可观测的变量(隐变量),那么单纯的极大似然估计法是不可用的:EM算法(e ...
论文R语言复现 | 基于 EM 算法的高斯混合模型参数估计
高斯混合概率在众多领域都有重要应用,依据已知观测数据估计高斯模型中未知参数就显得尤为重要,由于观测值具体来自于高斯分布的哪个分模型是未知的,那么利用传统的极大似然( MLE) 方法进行参数估计就变得十 ...
【信道模型分析1】通过伪码相关峰分析信道模型对信号的影响——高斯白噪声信道/瑞丽多径信道
目录一.理论基础二.核心程序三.仿真测试结果 3.1高斯白噪声信道模型
高斯噪声、高斯白噪声解析
如何用matlab 产生均值为0,方差为5的高斯噪声 2011-07-15 19:36 y=randn(1,2500); y=y/std(y); y=y-mean(y); a=0; b=5; y=a ...

高斯白噪声下基于EM的多径时延估计