范文一:随机信号处理
一种卡尔曼滤波与粒子滤波相结合的非线
性滤波算法
李玉莹
(河北科技大学,河北石家庄050000)
摘要:提出一种基于卡尔曼滤波与粒子滤波的非线性滤波算法。这种方法对于状态变量服从线性变化而观测方程为非线性的动态系统模型具有显著的效果。首先使用粒子滤波对状态变量进行初估计,然后对估计结果进行卡尔曼滤波,另外推导出该系统模型下状态变量估计误差的克拉美劳下界。通过计算复杂度分析及仿真实验验证,表明新方法与标准粒子滤波算法复杂度相当,但参数估计精度要高于标准粒子滤波以及扩展卡尔曼滤波算法,估计误差甚至要低于系统模型的克拉美劳下界。
关键词:非线性滤波卡尔曼滤波克拉美劳下界计算复杂度
A Nonlinear Filtering Algorithm Combining the KalmanFilter
and the Particle Filter
LI Yuying
(Hebei University of Science and Technology, Shijiazhuang, Hebei, 050000, China) Abstract:A nonlinear filtering algorithm is proposed based on the Kalman filter and the particle filter.The method can pro-vide significant performance for dynamic nonlinear system which is consist of linear state equation and nonlinear measurement equa-tion.Firstly,the particle filter is utilized for initial estimation of the state variables,and then the Kalman filter is performed.TheCramer-Rao Bound is derived for the nonlinear model.Computation complexity analysis and numerical simulations demonstrate thatthe proposed algorithm has the same complexity as the standard particle filter,but the estimation accuracy is higher than the standard article filter and the extended Kalman filter.The estimation error is even lower than the Cramer-Rao Bound of the system model. Keywords:nonlinear filteringKalman filter Cramer-Rao bound computation complexity
1 非线性滤波
近年来,非线性滤波问题越来越引起广泛的关注,其在目标跟踪、定位以及参数估计等方面都有广泛的应用。在贝叶斯框架下,给定测量样本,对一个非线性系统的状态估计可以
通过近似计算状态空间的后验概率密度函数来解决。对于线性高斯模型来说,卡尔曼滤波被公认为是最优滤波,其通过递推式更新有限维统计量来精确计算后验分布。但是大部分工程应用问题需要用非线性和/或者非高斯模型来建模。最简单的方式是利用一阶泰勒级数展开对模型进行线性化,再结合卡尔曼滤波,这种方法被称为扩展卡尔曼滤波。但是,当模型高度非线性或者后验分布呈现多峰情况时,该方法性能急剧下降。一种新的滤波方法,称为无迹卡尔曼滤波,使用无迹变换(UT)可以对非线性模型进行更加精确的近似。但是不管是扩展卡尔曼滤波还是无迹卡尔曼滤波,对于非高斯模型来说均不是最优滤波器。
粒子滤波是另外一类非线性滤波方法。该方法用一组采样点(也称为粒子)来近似表示目标状态后验概率密度函数,粒子通过重要性采样在动态系统中传递,并贯序地对后验分布进行更新。为防止粒子衰退,使用重采样技术对粒子进行更新。当粒子数量足够多时,该方法能获得最优结果。与卡尔曼滤波及改进方法相比,系统的非线性和非高斯性越强,粒子滤波的效果越明显。在标准粒子滤波基础上,又出现高斯粒子滤波、无迹粒子滤波和边缘化粒子滤波等。高斯粒子滤波通过重要性采样对未知状态变量的后验均值和方差进行递推式估计。其优点在于噪声可以是非高斯的,并且不需要粒子重采样,降低了运算复杂度。无迹粒子滤波是采用无迹卡尔曼滤波方法产生粒子的期望分布,是对系统状态后验概率密度函数的高斯近似。边缘化粒子滤波是对状态空间分为线性状态变量和非线性状态变量,线性部分由卡尔曼滤波估计而非线性部分由粒子滤波估计,其优点在于对高维状态变量进行降维处理并提高估计精度。
通过对以上粒子滤波及相关算法深入研究,我们发现虽然高斯粒子滤波算法在计算复杂度方面有较大优势,但是其估计精度与标准粒子滤波算法比较接近,而无迹粒子滤波和边缘化粒子滤波是分别对每个粒子进行卡尔曼滤波处理,计算量较大。在综合考虑计算复杂度与估计精度的基础上,本文提出一种卡尔曼滤波与粒子滤波相结合的非线性滤波算法,称为混合粒子滤波(MPF)。特别针对系统模型中状态方程为线性而观测方程为非线性的情况,在不增加计算复杂度的基础上,该方法能够提供较高的估计精度。这种模型有望在固定目标定位和参数估计等方面广泛的应用。算法首先采用粒子滤波对状态变量进行初始估计,由于状态变量服从线性变化,接下来使用卡尔曼滤波进一步滤波。本文对所提出算法的计算复杂度进行分析。另外,给出该系统模型下,状态变量估计误差下界,即克拉美劳下界(CRB)的推导过程。
MPF算法的优势在于:
(1)考虑受到噪声的影响,特别是较低信噪比条件下,标准粒子滤波性能受到限制,估
计结果会出现发散的情况。本文算法在粒子滤波的基础上采用卡尔曼滤波进行再处理,从而可以改善估计效果,提高精度。
(2)与卡尔曼滤波算法相比,本文方法可适用于较强非线性化系统模型,并且只需知道目标变量的先验分布,对初始条件的依赖性不强。而高度非线性系统以及不合适的初始条件会导致卡尔曼滤波算法发散、失效。
2 混合粒子滤波算法(MPF)
很多非线性估计问题可以通过粒子滤波方法来解决,本文所关注的状态空间模型可表示为
xn=Axn?1+un?1(1)
yn=? xn +vn(2)
其中xn∈Rm表示n时刻m维线性状态列向量,yn表示n时刻观测值,A∈Rm?m表示已知的状态转移矩阵,h(°)是任意关于xn的线性或者非线性函数,实际中根据不同的应用具有相应的解析表达,un与vn分别表示均值为0的高斯分布状态噪声和观测噪声,且满足
TE{unuk}=Qδnk,E{vnvk}=σ2δnk,?n,k其中Q为状态噪声协方差矩阵,σ2为观测噪声方差,假设状态噪声与观测噪声相互独立。
我们的目标是通过观测序列y0:n实现对状态x0:n的顺序估计。根据贝叶斯理论,后验概率p(x0:n|y0:n)包含了对状态x0:n估计的所有信息,但是我们很难直接得到其数学表达式。粒子滤波算法为解决该问题提供了有效的方法,其通过在一个所谓的重要性分布π(°)采样获
Ni得一组离散样本点{x0:n}i=1来近似计算后验概率,并得到粒子的重要性权重
iwn=ip x0:n y0:n
π x0:n y0:n =iwn?1iiip yn xn p xn xn?1
π xn yn,x0:n?1 i=1…,n 3
我们首先在已知的先验概率分布p(x0)上进行采样获得粒子,然后通过概率分布
iiiiπ xn yn,x0:n?1 =p xn xn?1 (4)
得到粒子的更新,这里由于状态噪声服从高斯分布,因此有
iiip xn xn?1 = N(Axn?1,Q) (5)
最后,由式(3)我们可以得到粒子的重要性权重为
iiiwn=wn?1p yn xn (6)
这里似然函数可以由下式得到
iip yn xn = N(h(xn),σ2) (7)
经过归一化
iw n=iwn/ N
i=1iwn (8)
及重采样,状态xn的最小均方估计(MMSE)可表示为
MMSEx n= N
i=1iiw nxn(9)
当观测噪声较大时,粒子无法准确描述后验概率,导致估计性能受到影响。考虑对状态xn的顺序估计满足
x n=Ax n?1+un?1(10)
MMSE并且基于粒子滤波的MMSE估计xn与x n之间足
MMSEx n= x n+ en(11)
这里en可以假设为最优估计与MMSE估计之间的误差,其服从均值为0,协方差矩阵为Υ的高斯分布。由于式(10)与式(11)均为线性方程,因此可以考虑卡尔曼滤波得到其最优估计 x n,有
Pn|n? 1= APn? 1|n? 1AT+ Q (12)
Kn=Pn|n? 1(?+Pn|n? 1)?1 (13)
Pn|n=(I? Kn)Pn|n? 1(14)
MMSEx n= x n?1+ Kn(x n? n?1) (15) 1? x
其中Pn|n? 1为先验估计误差协方差矩阵,Kn为卡尔曼增益,Pn|n为后验估计误差协方差矩阵。若粒子的初始分布服从均匀分布U(°),即
p(xO)= U(a,b) (16)
其中a,b∈ Rm 分别表示均匀分布的下界和上界,则卡尔曼滤波的初始值为
x 0=(a+ b)/2,P0|0=(b? a) b? a TI/2 (17)
本文所提出的混合粒子滤波算法可总结归纳如算法1。
算法1 混合粒子滤波(MPF)算法步骤
iN1. 初始化:在先验概率分布p(x0)采样得到粒子{x0}i=1
2. For n= 1,2,…
i3. 根据式(5)对粒子进行更新,得到xn
i4. 根据式(6)更新重要性权重,得到wn
i5. 根据式(8)对权重进行归一化,得到 w n
6. 对粒子进行重采样
MMSE7. 根据式(9)获得对状态xn的MMSE估计 x n
8. 根据式(12)得到先验估计误差协方差矩阵Pn|n? 1
9. 根据式(13)得到卡尔曼增益Kn
10.根据式(14)得到后验估计误差协方差矩阵Pn|n
11.根据式(15)得到状态xn的最优估计x n. 12. End n
3 算法性能分析
3.1 误差下界
本文所提出的混合粒子滤波算法对状态xn的估计收敛于最优估计。但是,当粒子数目N<>
对于固定参数的无偏估计问题,可以通过克拉美劳下界得到其方差下界,而对于贝叶斯模型来说,其状态迭代过程可以认为是一种随机过程,理论上的方差下界可通过后验克拉美劳下界(PCRB)获得。因此,我们希望推导出对xn估计的后验克拉美劳下界。
3.2 计算复杂度
在这一节中,我们对混合粒子滤波算法的计算复杂度进行分析。根据算法1,对于n时刻,我们首先分别用c1、c2、c3和c4表示产生均匀分布随机数、高斯分布随机数、计算高斯似然函数以及粒子重采样的计算量,若xk∈Rm,则算法1中第1行产生N个服从均匀分布的粒子需要计算量为mc1N。第3行对粒子更新需要计算量mc2N。第4行更新重要性权重,需要计算量(c3+1)N。第5行对重要性权重进行归一化,需要N次实数乘法。第6行粒子重采样需要计算量mc4N。第7行求最小均方误差估计需要mN次实数乘法。第8行计算先验估计误差协方差矩阵的复杂度,主要由矩阵A∈ Rm× m、Pn?1|n?1∈ Rm× m和AT∈Rm× m相乘所贡献,共需2m3次实数乘法。第9行计算卡尔曼增益,其中计算矩阵的逆需要计算量为m3,计算矩阵相乘需要的计算量为m3,则其总计算量为2m3。第10行计算后验估计误差协方差矩阵的计算量为m3。最后,第11行得到状态xk的最优估计需要m2次实数乘法。因此本问所提出的混合粒子滤波算法的计算复杂度可表示为
C m,N = c1m+ c2m+c4m+c3+m+2 N+5m3+m2(18)
在以上运算过程中忽略实数加法的计算量。由式(18)可知,对于n时刻,本文算法与标准粒子滤波算法相比,在计算复杂度上并未明显增加,均为?(N)。
4 实验仿真
本文所提出的混合粒子滤波算法是一种较为通用的参数估计及数字信号滤波的方法。我们考虑对叠加高斯分布噪声的周期信号中未知参数的估计问题,并与标准粒子滤波(SPF)以及扩展卡尔曼滤波(EKF)算法进行对比。若观测信号为
yn=αncos(2π yfn+φn)+ vn(19)
其中已知频率f= 1KHz与采样周期ΔT= 10?4s为已知参数,幅度αn及相位φn是未知参数,vn服从均值为0,方差为σ的高斯分布,取σ= 100对未知参数建模,得到状态方程有
ααn=αn?1+un?1 (20) φφn=φn?1+un?122
我们将待估计的参数表示为xn=[αn,φn]T,则由式(10)可知状态矩阵A=[1,0;0,1],设定状态噪声协方差矩阵Q=[10?4,0;0,10?4]。参数的真实值为xn=[8,2π/3]T,其中,幅度参数可通过计算观测信号能量得到其分布,这里选择均匀分布,相位参数满足0到2π的均匀分布。设定未知参数矢量的先验分布满足x0~ U([0,0]T,[15,2π]T),卡尔曼滤波的观测噪声协方差矩阵为?=[10,0;0,5],持续时间为T= 0。4s,粒子数量N= 400。如果没有特殊说明,以上参数为默认值。这里给出一次仿真实验的参数估计结果,将MPF、SPF以及EKF进行对比,如图1所示。我们可以看到,这三种算法均能够渐进地收敛到目标参数的理论值,但是SPF和EKF算法的估计结果不够稳定,在理论值上下浮动,而MPF是对SPF算法的进一步滤波,其结果更加可靠、精确。
给出在不同粒子数量情况下,MPF和SPF算法的估计误差对比。我们针对粒子的不同取值,分别进行M=100次独立实验,并对估计结果取均值。可以看到,不管是MPF还是SPF,对相位的估计精度要高于对幅度的估计精度。在粒子数量较少时,两种方法的性能比较接近,当粒子数量逐渐增多,MPF算法的估计误差要明显低于SPF算法,并且当粒子数为100时,MPF的估计性能可以达到稳定状态。综合考虑计算复杂度及估计精度等因素,选择N= 200作为本文算法的最佳粒子数。
5 结束语
非线性滤波问题是工程应用领域十分关注的一个方向。在构建系统模型时,会遇到状态方程是线性的而观测方程属于非线性的问题。本文正是针对这种情形,提出一种新的非线性滤波算法。该方法将粒子滤波与卡尔曼滤波有机地结合起来。首先采用粒子滤波对状态变量进行初估计,再通过卡尔曼滤波进行线性递推估计。同时,详细地分析了新算法的计算复杂
度。通过周期信号未知参数估计问题的计算机仿真实验,将本文算法的估计误差与标准粒子滤波的估计误差、克拉美劳下界以及扩展卡尔曼滤波算法估计误差进行对比,结果清晰地显示了新算法具有更高的估计精度,特别是在较低信噪比情况下,仍能达到比较理想的估计结果,并且对粒子数量的要求比较低,计算复杂度与标准粒子滤波算法相当。本文算法可广泛适用于信号处理中固定未知参数估计等问题。由于本文的分析与推导是建立在高斯噪声假设的基础上,未来的工作可以针对非高斯背景下的参数估计进行深入的研究。
参考文献/References:
[1]刘先省,胡振涛,金勇,杨一平.基于粒子优化的多模型粒子滤波算法[J].电子学报,2010,38(2):301- 306.
[2]刘义,赵晶,冯德军,王雪松,王国玉.高精度惯导速度信息辅助的弹目相对运动模型构建方法[ J].电子学报,2011,39(9):2207- 2211.
[3]吴孙勇,廖桂生,杨志伟.基于粒子滤波的宽带信号波达方向估计[J].电子学报,2011,39(6):1353- 1357.
[4] F Daum.Nonlinear filters beyond the Kalman filter[ J].IEEEA&E Systems Magazine,2005,20(8):57- 69.
[5] S J Julier,J K Uhlmann.Unscented filtering and nonlinear esti-mation[J].Proceedings of IEEE,2004,92(3):401- 422.
范文二:随机信号处理
学院:电子工程学院
、
马尔可夫过程概述
摘要:叙述了随机过程中的某一种--马尔可夫过程的基本定义 ,特点,以及它的应用领域;通过对离散时间马尔可夫链进行仿真分析,掌握马尔可夫的特点。 1. 随机过程发展简述
在当代科学与社会的广阔天地里,人们都可以看到一种叫作随机过程的数学模型:从银河亮度的起伏到星系空间的物质分布、从分子的布朗运动到原子的蜕变过程,从化学反应动力学到电话通讯理论、从谣言的传播到传染病的流行、从市场预测到密码破译,随机过程理论及其应用几乎无所不在。
一些特殊的随机过程早已引起注意,例如1907年前后,Α.Α.马尔可夫研究过一列有特定相依性的随机变量,后人称之为马尔可夫链(见马尔可夫过程);又如1923年N.维纳给出了布朗运动的数学定义(后人也称数学上的布朗运动为维纳过程),这种过程至今仍是重要的研究对象。虽然如此,随机过程一般理论的研究通常认为开始于30年代。1931年,Α.Η.柯尔莫哥洛夫发表了《概率论的解析方法》;三年后,Α.Я.辛钦发表了《平稳过程的相关理论》。这两篇重要论文为马尔可夫过程与平稳过程奠定了理论基础。稍后,P.莱维出版了关于布朗运动与可加过程的两本书,其中蕴含着丰富的概率思想。1953年,J.L.杜布的名著《随机过程论》问世,它系统且严格地叙述了随机过程的基本理论。1951年伊藤清建立了关于布朗运动的随机微分方程的理论(见随机积分),为研究马尔可夫过程开辟了新的道路;近年来由于鞅论的进展,人们讨论了关于半鞅的随机微分方程;而流形上的随机微分方程的理论,正方兴未艾。60年代,法国学派基于马尔可夫过程和位势理论中的一些思想与结果,在相当大的程度上发展了随机过程的一般理论,包括截口定理与过程的投影理论等,中国学者在平稳过程、马尔可夫过程、鞅论、极限定理、随机微分方程等方面也做出了较好的工作。 2. 马尔可夫过程发展
2.1 马尔可夫过程简介
马尔科夫过程(MarKov Process)是一个典型的随机过程。设X(t)是一随机过程,当过程在时刻t0所处的状态为已知时,时刻t(t>t0)所处的状态与过程在t0时刻之前的状态无关,这个特性成为无后效性。无后效的随机过程称为马尔科夫过程。马尔科夫过程中的时同和状态既可以是连续的,又可以是离散的。我们称时间离散、状态离散的马尔科夫过程为马尔科夫链。马尔科夫链中,各个时刻的状态的转变由一个状态转移的概率矩阵控制。 2.2 马尔可夫过程的发展
20世纪50年代以前,研究马尔可夫过程的主要工具是微分方程和半群理论(即分析方法);1936年前后就开始探讨马尔可夫过程的轨道性质,直到把微分方程和半群理论的分析方法同研究轨道性质的概率方法结合运用,才使这方面的研究工作进一步深化,并形成了对轨道分析必不可少的强马尔可夫性概念。1942年,伊藤清用他创立的随机积分和随机微分方程理论来研究一类特殊而重要的马尔可夫过程??扩散过程,开辟了研究马尔可夫过程的又一重要途径。
出于扩大极限定理应用范围的目的,马尔科夫在20世纪初开始考虑相依随机变量序列的规律,并从中选出了最重要的一类加以研究。1906年他在《大数定律关于相依变量的扩展》一文中,第一次提到这种如同锁链般环环相扣的随机变量序列,其中某个变量各以多大
的概率取什么值,完全由它前面的一个变量来决定,而与它更前面的那些变量无关。这就是被后人称作马尔科夫链的著名概率模型。也是在这篇论文里,马尔科夫建立了这种链的大数定律。
用一个通俗的比喻来形容,一只被切除了大脑的白鼠在若干个洞穴间的蹿动就构成一个马尔科夫链。因为这只白鼠已没有了记忆,瞬间而生的念头决定了它从一个洞穴蹿到另一个洞穴;当其所在位置确定时,它下一步蹿往何处与它以往经过的路径无关。这一模型的哲学意义是十分明显的,用前苏联数学家辛钦(1894,1959〕的话来说,就是承认客观世界中有这样一种现象,其未来由现在决定的程度,使得我们关于过去的知识丝毫不影响这种决定性。这种在已知“现在”的条件下,“未来”与“过去”彼此独立的特性就被称为马尔科夫性,具有这种性质的随机过程就叫做马尔科夫过程,其最原始的模型就是马尔科夫链。
这既是对荷兰数学家惠更斯(Ch. Huygens, 1629,1659)提出的无后效原理的概率推广,也是对法国数学家拉普拉斯(P. S. Laplace, 1749,1827)机械决定论的否定。
这里应该指出,尽管拉普拉斯对概率论的早期发展作出过重大贡献,但是他的部分哲学观点是不利于这门学科的深入发展的。十八世纪以来,随着牛顿力学的彻底胜利,一种机械唯物主义的决定论思潮开始在欧洲科学界蔓延,鼓吹最力者就是拉普拉斯。1759年他在巴黎高等师范学院发表了一篇题为《概率论的哲学探讨》的演讲,淋漓尽致地表达出了这种思想。他说:“假如有人知道了某一时刻支配自然的一切力,以及它的一切组成部分的相对位置,又假如他的智力充分发达,能把这一切数据加以充分的分析,把整个宇宙中从最巨大的天体到最微小的原子的一切运动完全包括在一个公式里面,这样对他就没有什么东西是不确定的了,未来也好,过去也好,他都能纵览无遗。”1812年,拉普拉斯又进一步提出“神圣计算者”的观念,认为这个理想的数学家只须知道世界某一时刻的初始状态,就可以从一个无所不包的微分方程中算出过去和未来的一切状态。换句话说,他认为任意系统在 t > t0时的状态 x可由其初始时刻 t0和初始状态 x0唯一决定。这可真是笔判终身、细评流年,数学家可以摆个卦摊了。马尔科夫的概率模型从根本上否定了系统中任一状态 x与其初始状态 x0之间的因果必然性,从而也否定了“神圣计算者”的神话。
还应该指出,马尔科夫所建立的概率模型不但具有深刻的哲学意义,而且具有真实的物质背景,在他的工作之前或同时,一些马尔科夫链或更复杂的随机过程的例子已出现在某些人的研究中,只不过这些人没有自觉地认识到这类模型的普遍意义或用精确的数学语言表述出来罢了。例如苏格兰植物学家布朗 ( R. Brown, 1773,1858) 于1827年发现的悬浮微粒的无规则运动、英格兰遗传学家高尔顿(F.Galton, 1822,1911) 于1889年提出的家族遗传规律、荷兰物理学家埃伦费斯特 ( P. Ehrenfest, 1880,1933) 于1907年关于容器中分子扩散的实验,以及传染病感染的人数,谣言的传播,原子核中自由电子的跃迁,人口增长的过程等等,都可用马尔科夫链或过程来描述。也正是在统计物理、量子力学、遗传学以及社会科学的若干新课题、新事实面前,决定论的方法显得百孔千疮、踵决肘见。
有趣的是,马尔科夫本人没有提到他的概率模型在物理世界的应用,但是他利用了语言文学方面的材料来说明链的性质。在《概率演算》第四版中,他统计了长诗《叶甫盖尼?奥涅金》中元音字母和辅音字母交替变化的规律:这是长诗开头的两句,意为:“我不想取悦骄狂的人生,只希望博得朋友的欣赏。”诗人那火一般的诗篇在数学家那里变成了一条冷冰冰的锁链:在这条锁链上只有两种链环,C代表辅音、 代表元音(为了使问题简化起见,不仿把两个无音字母算作辅音)。马尔科夫分别统计了在C后面出现C和 的概率p和1,p,以及在 后出现C和 的概率q和1,q,把结果与按照俄语拼音规则计算出的结果进行比较,证实了语言文字中随机的(从概率的意义上讲)字母序列符合他所建立的概率模型。
完成了关于链的大数定律的证明之后,马尔科夫又开始在一系列论文中研究链的中心极限定理。1907年他在《一种不平常的相依试验》中证明了齐次马尔科夫链的渐近正态性;
1908年在《一个链中变量和的概率计算的极限定理推广》中作了进一步的推广;1910年他发表了重要的论文《成连锁的试验》,在其中证明了两种情况的非齐次马尔科夫链的中心极限定理。与此同时他在一些假定的前提下证明了模型的各态历经性,成为在统计物理中具有重要作用的遍历理论中第一个被严格证明的结果。遍历理论亦称ergodic理论, 是奥地利物理学家玻耳兹曼(L. Boltzmann, 1844,1906) 于1781年提出来的,其大意是:一个系统必将经过或已经经过其总能量与当时状态相同的另外的任何状态。
马尔科夫链的引入,在物理、化学、天文、生物、经济、军事等科学领域都产生了连锁性的反应,很快地涌现出一系列新的课题、新的理论和新的学科,并揭开了概率论中一个重要分支,,随机过程理论蓬勃发展的序幕。
3 马尔可夫过程的应用
3.1 马尔可夫应用概述
马尔可夫随机过程的发展史说明了理论与实际之间的密切关系。许多研究方向的提出,归根到底是有其实际背景的。反过来,当这些方向被深入研究后,又可指导实践,进一步扩大和深化应用范围。下面简略介绍一下马尔可夫随机过程本身在各方面的应用情况。
在物理学方面,高能电子或核子穿过吸收体时,产生级联(或倍增)现象,在研究电了-光子级联过程的起伏问题时,要用到随机过程,常以泊松过程、弗瑞过程或波伊亚过程作为实际级联的近似,有时还要用到更新过程(见点过程)的概念。当核子穿到吸收体的某一深度时,则可用扩散方程来计算核子的概率分布。物理学中的放射性衰变,粒子计数器,原子核照相乳胶中的径迹理论和原子核反应堆中的问题等的研究,都要用到泊松过程和更新理论。湍流理论以及天文学中的星云密度起伏、辐射传递等研究要用到随机场的理论。探讨太阳黑子的规律及其预测时,时间序列方法非常有用。
化学反应动力学中,研究化学反应的时变率及影响这些时变率的因素问题,自动催化反应,单分子反应,双分子反应及一些连锁反应的动力学模型等,都要以生灭过程(见马尔可夫过程)来描述。
随机过程理论所提供的方法对于生物数学具有很大的重要性,许多研究工作者以此来构造生物现象的模型。研究群体的增长问题时,提出了生灭型随机模型,两性增长模型,群体间竞争与生尅模型,群体迁移模型,增长过程的扩散模型等等。有些生物现象还可以利用时间序列模型来进行预报。传染病流行问题要用到具有有限个状态的多变量非线性生灭过程。在遗传问题中,着重研究群体经过多少代遗传后,进入某一固定类和首次进入此固定类的时间,以及最大基因频率的分布等。
许多服务系统,如电话通信,船舶装卸,机器损修,病人候诊,红绿灯交换,存货控制,水库调度,购货排队,等等,都可用一类概率模型来描述。这类概率模型涉及的过程叫排队过程,它是点过程的特例。排队过程一般不是马尔可夫型的。当把顾客到达和服务所需时间的统计规律研究清楚后,就可以合理安排服务点。
在通信、雷达探测、地震探测等领域中,都有传递信号与接收信号的问题。传递信号时会受到噪声的干扰,为了准确地传递和接收信号,就要把干扰的性质分析清楚,然后采取办法消除干扰。这是信息论的主要目的。噪声本身是随机的,所以概率论是信息论研究中必不可少的工具。信息论中的滤波问题就是研究在接收信号时如何最大限度地消除噪声的干扰,而编码问题则是研究采取什么样的手段发射信号,能最大限度地抵抗干扰。在空间科学和工业生产的自动化技术中需要用到信息论和控制理论,而研究带随机干扰的控制问题,也要用到马尔可夫随机过程。
3.2 一种新的马尔可夫模型应用举例
隐马尔可夫模型(Hidden Markov Model,HMM)是统计模型,它用来描述一个含有隐含未知参数的马尔可夫过程。其难点是从可观察的参数中确定该过程的隐含参数。然后利用这些参数来作进一步的分析,例如模式识别。在正常的马尔可夫模型中,状态对于观察者来说是直接可见的。这样状态变迁概率便是全部的参数。而在隐马尔可夫模型中,状态并不是直接可见的,但受状态影响的某些变量则是可见的。每一个状态在可能输出的符号上都有一概率分布。因此输出符号的序列能够透露出状态序列的一些信息。
隐马尔可夫模型是马尔可夫链的一种,它的状态不能直接观察到,但能通过观测向量序列观察到,每个观测向量都是通过某些概率密度分布表现为各种状态,每一个观测向量是由一个具有响应概率密度分布的状态序列产生。所以,隐马尔可夫模型是一个双重随机过程----具有一定状态数的隐马尔可夫链和显示随机函数集。自20世纪80年代以来,HMM被应用于语音识别,取得重大成功。到了90年代,HMM还被引入计算机文字识别和移动通信核心技术“多用户的检测”。近年来,HMM在生物信息科学、故障诊断等领域也开始得到应用。 4 马尔可夫过程的仿真分析
源程序:
1. i=1; %设置当前第一个状态i为1
2. j=1; %设置初始状态
3. N=100; %设定步长N
4. n=0; %n初始化,对计数器n清零
5. s=0; %累加器初始化清零
6. p=[0.2 0.3 0.5;0.5 0.1 0.4;0.6 0.2 0.2]; %p为一个3X3维矩阵
,然后重复 7. while n
8. n=n+1; %计数器n加1
9. u=rand(); %每次循环产生一个随机数
10.i=j; %将j赋给i,接下来执行第i行
11.q(n)=i; %将x(n)放在一个一维数组中
12.j=1; %每一次外循环后再次对j初始化赋值为1
13.s=0; %依次循环后,累加器清零
14.while j<=size(p,2)&& u="">s %内重循环,找出符合要求的j值 15.s=s+p(i,j); %概率和不断累加
16.j=j+1; %计数器j加1
17.end %内重循环结束
18.j=j-1; %得到符合要求的j,也即为Y(i)
19.end
20.stairs(q); %调用画图阶梯函数
21.axis([0 100 0 5]); %标注输出的图线的最大值最小值。 22.xlabel('步长n');ylabel('x(n)');legend('样本函数1'); %标记x、y轴, 结果分析
马尔可夫过程(Markov process)是一类随机过程。它的原始模型马尔可夫链,由俄国数学家A.A.马尔可夫于1907年提出。该过程具有如下特性:在已知目前状态 (现在)的条件下,它未来的演变 (将来)不依赖于它以往的演变 ( 过去 ) 。 例如森林中动物头数的变化构成--马尔可夫过程 。在现实世界中,有很多过程都是马尔可夫过程,如液体中微粒所作的布朗
运动、传染病受感染的人数、车站的候车人数等,都可视为马尔可夫过程
马尔可夫过程的定义还可以进一步扩充。第一,所谓"过去"可以作更广泛的理解,即(2)中由,Xs所产生的σ域(见概率)可以扩大为一般的σ域Fs,只要Fs包含由{X,u?s}产生的σ域,而当 s 马尔可夫过程的许多性质可以通过转移函数来表达。转移函数P(s,x,t,A)(0?s?t,x?E,A?B)是满足某些条件的四元函数,它可以理解为过程在时刻s时处在x,在时刻t 时转移到A中的条件概率。如果P(s,x,t,A)=P(t-s,x,A)只依赖于t-s,x及A,则称转移函数及相应的马尔可夫过程为齐次的。设E是d维欧几里得空间Rd,B为Rd中的波莱尔域(见概率分布)Bd,而且齐次转移函数满足下面的登金-金尼条件:对任意 ε>0,?。式中Vε(x)={y:|y-x|?ε},那么可以选取轨道连续的齐次马尔可夫过程X,以p(t,x,A)为转移函数。一类重要的轨道连续马尔可夫过程是 d维布朗运动。 参考文献 1.《概率论与随机过程》,北京邮电大学出版社,2003 2.《随机信号分析》,科学出版社,2009 3.《随机信号分析与处理》,清华大学出版社,2006 参数功率谱估计的实现 及其与 经典谱估计的比较 参数功率谱估计的实现及其与经典谱估计的比较 一、摘要 频域分析是从频率角度对信号进行分析研究,对于确定性信号来说,通常使用傅立叶变换或者傅立叶技术展开得到频域表达,但对随机信号而言,由于其时域波形的随机性及能量无限,没有确定的时域表达式,无法用傅立叶变换直接将其变换到频域中去研究。表达随机信号通常使用概率密度函数,根据维纳-辛钦定理,广义随机过程的功率谱与自相关函数是一对傅立叶变换,所以自然想到可以用功率谱来研究随机过程的频域性质。 对随机信号的功率谱估计方法通常分为两大类:经典谱估计和现代谱估计经典谱估计是基于维纳-辛钦定理,从自相关函数出发通过傅立叶变换得到功率谱,而现代谱估计则将随机信号看成白噪声通过一个滤波器的输出。现代谱估计就是通过记录的信号序列估计滤波器参数,从而得到其频率响应,最后通过S (e jw 2 ) =σm |H (e jw ) |2得到其功率谱。 从频率分辨率来看,经典谱估计效果一般不如参数谱估计好,而且参数谱估计的不同算 法对不同的采样序列也有不同的效果,本文将对周期图法,基于L-D 快速递推算法的Y-W 法和Burg 算法使用Matlab 进行编程实现并做比较。 二、关键字 功率谱估计 周期图法 L-D 递推 Y-W 算法 Burg 算法 三、原理 1. 经典谱估计 广义平稳随机信号经典谱估计基于维纳-辛钦定理:广义平稳随机的自相关函数与其功率谱是一堆傅立叶变换。所以要求功率谱,只需由随机序列求出自相关函数然后进行傅立叶变换即可, 1N -1 R (m ) =x (n ) x (n +m ) m=0,1,2,……N-1 ∑N n =0 ∧ S (k ) =FFT (R (m )) 周期图法: 截断RN (n ) ∧ 虽然经典谱估计方法比较直观简单,但由于随机序列相当于对信号加床,所以求自相关函数后傅立叶变换的功率谱往往受到窗函数影响,不是信号真实谱,所以就产生了以下的现代谱估计。 参数功率谱估计: 周期图法谱估计运算框图 平稳随机信号模型 X(n) 根据谱分解定理,任何平稳随机信号x (n ) 都可以看成是由白噪声经过一个因果稳定可逆系统H (z ) 产生的输出,如上图。其系统函数H (z ) 全部极点都在z 平面单位圆内部,将 H (z ) 表示成H (z ) = B (z ) ,实际中,很多物理过程都可以近似表示成AR ,ARMA ,MA A (z ) 过程,其中AR 过程为H (z ) = G 1+∑a i z -i i =1,也就是可以看成白噪声通过系统后的输出。 如果根据有限的数据记录(看成系统输出)和已知输入(白噪声)确定出系统的参数,那么所要研究的随机过程就是白噪声通过已知系统之后的输出,于是可以用系统的参数来描述索要研究的物理过程,例如要分析实平稳随机信号x (n ) 的有理功率谱S x (e 机信号通过线性系统的性质可知:S x (e 率谱(为常数),|H(ejw ) |是H (e jw jw jw ) ,根据平稳随 22 ) =σw |H (e jw ) |2,其中σw 是白噪声w (n ) 的功 单位脉冲响应h (n ) ) 的能量谱对线性时不变的稳定系统, jw 是确定性能量信号,其能量谱是频谱H (e 率谱S x (e jw ) 幅度的平方。上式说明平稳随机信号x (n ) 的功 ) 可以用模型H (e jw ) 的参数来表示。所以只要建立合适的信号模型,并估计出 模型的参数,就可以分析随机信号的功率谱。上面仅提到了AR 模型,是因为柯尔莫哥洛夫 定理表明,任何ARMA 或MA 过程都可以用无限阶的AR 过程表示,即使在建模式选择了不正确的模型,只要模型阶数足够高,均可以得到一个合理的近似表示。 AR 模型的正则方程和参数计算: 正则方程,是在均方误差最小准则下建立的模型参数与相关函数的关系他是在二阶统计意义上的建模,由于高斯过程仅用二阶统计量就可以完全描述,所以这种基于相关函数的建模方法只适合于高斯过程。AR 模型的参数可以借助解线性方程获得,而且MA 和ARNA 模型都可以用阶数足够高的AR 模型表示,所以实际中通常选用AR 模型。 Yule-Walker 方程: 假定w (n ) 和x (n ) 都是实平稳随机信号,其中w (n ) 是方差为1的高斯白噪声序列, x (n ) 为AR 过程,则可将系统输出表示为: x (n ) =-∑a i x (n -i ) +Gw (n ) i =1 p 即用n 时刻之前p 个值的线性组合- ∑a x (n -i ) 来预测x (n ) ,预测误差为Gw (n ) , i i =1 p 在均方误差最小准则下,a 1, a 2, ?, a p 的选择应使Gw (n ) =e (n ) 均方值E (e 2(n ) ) 最小。 ?E [e 2(n )] 据此准则有=0,经过推导得到AR 模型的正则方程,即Y-W 方程: ?a i ?p -a i R x (m -i ) , m =1, 2, ?, p ??∑R x (m ) =?i p =1 ?-∑a i R x (i ) +G 2, m =0??i =1 直接解Y-W 方程的计算量一般较大,这里使用L-D 快速递推算法: p -1 ? ?a p (p ) =-[R x (p ) +∑a p -1(i ) R x (p -i ) ]/ρp -1 i =1 ? ?a p (i ) =a p -1(i ) +a p (p ) a p -1(p -i ), i =1, 2,..., p -1 ?ρp =ρp -1[1-a 2p (p )]? ? 其中m =1, 2,..., p ,递推初始条件:ρ0=R x (0) ,a 1(1) =- R x (1) ,自相关函数R x (m ) 求R x (0) 1N -m -1 法:R x (m ) =x (n ) x (n +m ) ,m ≥0且R x (-m ) =R x (m ) ,这样做是保证求取自∑N n =0 相关函数只用到观测数据,而不对序列补零。从观测得到随机序列x (n ) 求得R x (m ) ,在代入递推公式即可求得a p (1), a p (2),..., a p (p ) 及G 从而得到H (z ) =得到H (e jw 2 |H (e jw ) |2求得S x (e jw ) 。 ) ,再由S x (e jw ) =σw B (z ) ,调用freqz 函数后A (z ) Burg 算法: Burg 算法与L-W 算法的区别在于它没有直接使用观测随机序列的自相关函数求取AR 模型的参数,而实现估计反射系数k m (即AR 模型参数中的a m (m ), m =1, 2,..., p ) 。再利用Levinson 关系式a m (i ) =a m -1(i ) +k m a m -1(i -1) 后向预测误差递推公式为: f f b ?e m (n ) =e m -1(n ) +k m e m -1(n -1) ?b b f ?e m (n ) =e m -1(n -1) +k m e m -1(n ) f b ?e (n ) =e (n ) =x (n ) 00? 前向和, i =1, 2,..., m -1求AR 模型参数。 反 f m 2 射系数估计 b m 准 2 则是使前向预测误差功率 f b ?E {[e m (n )]2+[e m (n )]2} e (n ) 和后向预测误差功率e (n ) 的平均值最小。即=0 ?k m f b 2∑[e m -1(n ) e m -1(n -1)]n =m f m -1 b 2 (n )]2+[e m -1(n -1)]}N -1 得到k m =-N -1 n =m ∑{[e ,m =1, 2,..., p Burg 算法同样不对未知数据做人为假设。该算法的优点是求得的AR 模型保证稳定, 另外,由于Burg 算法不估算自相关函数,所以性能优于自相关法。尤其在短数据时具有较高分辨率。 四、实现程序 clear all %信号产生,两个频率的正弦信号加白噪声 t0=1; %采样间隔 tt=128; %采样时长 N=tt/t0; %采样点数-1 M=27; %模型阶数 t=0:t0:tt; f1=0.2; %信号频率 f2=0.213; y=sin(2*pi*f1*t)+sin(2*pi*f2*t)+randn(1,N+1); plot(t,y) title(' 信号加高斯噪声' ) %经典谱估计 F=fft(y,1024); F=(F.^2)/1024; m=(0:1023)/(1024); figure; plot(m,10*log(abs(F))); title(' 周期图法经典谱估计' ) xlabel(' 对采样频率归一化' ) ylabel('10lg(H)') w0=Freqval(abs(F),1/t0) %基于Y-W 方程的参数谱估计 %采用L-D 快速递推算法实现 R=zeros(1,N+1); for i=0:N %计算相关函数 for r=0:(N-i) R(i+1)=R(i+1)+y(r+1)*y(r+i+1); end R(i+1)=R(i+1)/(N+1); end a(1,1)=(-1)*R(2)/R(1); P(1)=R(1); P(2)=R(1)*(1-(a(1,1))^2); p=1; while ((abs(a(p,p)))<1 &&="">1><=m))%此处p 表示每次运算="" p="p+1;" %中最高阶数,反射系="" x="0;" %数k="" ,而直接是a(p,p)="" for="" i="">=m))%此处p> X=X+a(p-1,i)*R(p-i+1); end a(p,p)=(-1)*(R(p+1)+X)/P(p); for i=1:p-1 a(p,i)=a(p-1,i)+a(p,p)*a(p-1,p-i); end P(p+1)=P(p)*(1-(a(p,p))^2); h=p-1; %h表示模型最高阶,一般为M end for i=1:h A(i)=a(h,i); end G=P(p)^(1/2); A=[1 A]; H=freqz(G,A,1024,'whole' ); figure; plot(m,20*log(abs(H))); title(' 基于Y-W 方程的参数谱估计' ) xlabel(' 对采样频率的归一化' ) ylabel('10lg(H)') w1=Freqval(abs(H),1/t0) %Burg算法实现信号参数谱估计 p=1; fenzi=0; fenmu=0; ef=zeros(p,N+1); eb=zeros(p,N+1); ef(1,:)=y; eb(1,:)=y; for i=p:N fenzi=fenzi+ef(1,i+1)*eb(1,i); fenmu=fenmu+ef(1,i+1)^2+eb(1,i)^2; a(1,1)=(-2)*fenzi/fenmu; while (a(p,p)<1 &&="">1><=m) fenzi="0;" fenmu="0;" p="p+1;" for="" i="">=m)> ef(p,i+1)=ef(p-1,i+1)+a(p-1,p-1)*eb(p-1,i); eb(p,i+1)=eb(p-1,i)+a(p-1,p-1)*ef(p-1,i+1); end for i=p:N fenzi=fenzi+ef(p,i+1)*eb(p,i); fenmu=fenmu+ef(p,i+1)^2+eb(p,i)^2; end a(p,p)=(-2)*fenzi/fenmu; end h=p-1; for i=1:h-1 a(h,i)=a(h-1,i)+a(h,h)*a(h-1,h-i); A(i)=a(h,i); end A(h)=a(h,h); A=[1 A]; G=R(1)^(1/2); H=freqz(G,A,1024,'whole' ); figure; plot(m,20*log(abs(H))); title('Burg 算法' ); xlabel(' 对采样频率的归一化' ); ylabel('10lg(H)'); w2=Freqval(abs(H),1/t0) %FreqVal函数用于测定算法检测到的信号频率,设计为监测两个 function z=FreqVal(x,y)%x为函数序列,y 为采样频率 q=length(x)/2; w=zeros(1,2); z=zeros(1,2); for j=2:q-1 if x(j)>=x(j-1) && x(j)>=x(j+1) if x(j)>w(1) w(1)=x(j); z(1)=j; end end for j=2:q-1 if x(j)>=x(j-1) && x(j)>=x(j+1) if x(j)>w(2) if j~=z(1) w(2)=x(j); z(2)=j; end end end end z=(z.*y)./(2*q); 五、实验结果简析 程序一次执行结果如下: 信号加高斯噪声 6 4 2 -2 -4 -6 020406080100120140 周期图法经典谱估计 20 -20 10l g (H ) -40 -60 -80 -100 00.2 0.40.6对采样频率归一化 0.81 W0=0.2002 0.2158 基于Y-W 方程的参数谱估计 3530 252015 10l g (H ) 1050-5-10-150 0.2 0.40.6对采样频率的归一化 0.8 1 w1 =0.1992 0.2139 50 Burg 算法 40 30 10l g (H ) 20 10 -10 00.2 0.40.6对采样频率的归一化 0.81 w2 =0.2139 0.1973 从以上结果可以对三种方法做一比较,随即序列表达式为 y (t ) =sin(2πf 1t ) +sin(2πf 2t ) +w (n ) ,从检测到两个正弦信号的频谱幅度来看周期图法结 果约为10dB ,而Y-W 法两峰值在25~30dB之间,明显高于前者,Burg 法可以达到接近40dB 。再从谱估计图两峰值图形比较,由于加窗效应,周期图法的宽度是最宽的,而Y-W 法峰宽度比较窄,效果好于周期图法,相比之下,Burg 法峰宽度最窄,效果在三者中最好。 最后再看三种方法得到的精确估计值与真实值的比较,随机信号谐波频率真实值: f 1=0. 2, f 2=0. 213,而w0 =0.2002 0.2158,w1 =0.1992 0.2139,w2 =0.2139 0.1973, Y-W 法和Burg 法对0.213的估计好于周期图法,对0.2的估计,周期图法好于Y-W 法,而Burg 法较差。当然,这只是一次实验的结果,要想对三种方法对随机信号频率精确估值效果进行研究,可以进行大量实验,再求取各自均值和方差等统计量进行研究。 六.总结 通过本次作业实验,对使用Y-W 方法和和Burg 算法对随机信号进行谱估计有了更深的,更清晰的认识,报告撰写过程中,需要不断阅读研究有关周期图、Y-W 算法、L-D 快速递推法、预测误差滤波器及Burg 算法的基本原理和公式推导过程,加深了对这些理论原理和各种特点的认识。在实验程序编写过程中,通过对各种算法的具体实现,不仅掌握了随机信号现代谱估计Y-W 法和Burg 法的实现过程,也对matlab 的使用更加熟悉,当然,实验程序还有很多可以优化的地方,但总的来说是一个进步。 随机信号处理实验 ———— 线性FM 信号通过匹配滤波器 姓名: 黄 越 学号: 0804210204 一、 实验目的: 1、了解线性FM 信号的产生及其性质; 2、熟悉MATLAB 的基本使用方法; 3、利用MATLAB 语言编程匹配滤波器。 4、仿真实现FM 信号通过匹配滤波器实现脉压处理,观察前后带宽及增益。 5、 进一步了解雷达中距离分辨率与带宽的对应关系。 二、 实验内容: 1、 线性FM 信号的产生。 LFM 信号(也称Chirp 信号) 的数学表达式为: s (t ) =rect ( t T ) e j 2π(f c t + k 2t ) 2 rect ( ) T 为矩形信号, t 式中 f c 为载波频率, ?t 1 , ≤1t ? rect () =?T T ?0 , elsew ise ? K = B T ,是调频斜率,于是,信号的瞬时频率为 f c +K t (-T 2 ≤t ≤T ) 2 s (t ) =S (t ) e j 2πf c t 当TB>1时,LFM 信号特征表达式如下: S LFM = 2k rect ( f -f c B + (f ) ) φLFM (f ) = π(f -f c ) μt T π4 2 S (t ) =rect () e j πK t 对于一个理想的脉冲压缩系统,要求发射信号具有非线性的相位谱,并使其包络接近矩形; 其中S (t ) 就是信号s(t)的复包络。 本例中产生的线性FM 信号,时宽为25us ,带宽为 2MHz 应用了dchirp 函数: function sn=dchirp(TW,p) a=1/(2*p^2*TW); N=p*TW; n=0:N-1; sn=exp(j*2*pi*a*(n-0.5*N).^2); 生成线性FM 信号 本例中产生的线性FM 信号,时宽为25us ,带宽为 2MHz Matlab 程序如下: p=20/2; tw=50; s=dchirp(tw,p); figure(1) t=25e-6; N=p*tw; t=[-t/2:t/N:t/2];t(end)=[]; plot(t*1e6,real(s)) xlabel('t (\mus)'),ylabel('The real of s(n)') title ('The chirp Singal') figure(2) fs=p*tw/t; y=fft(y); y=fftshift(y); N1=length(y); f=[-N1/2:N1/2]/N1*fs;f(end)=[]; plot(f/1e6,abs(y)),grid on title ('The Frequence of Chirp signal') xlabel('f (Mhz)'),ylabel('Magnitude') 仿真得时域波形: 得频域波形: 2、 设计匹配滤波器并将上述信号通过此滤波器 窄脉冲具有宽频谱带宽,如果对宽脉冲进行频率、相位调制,它就可以具有和窄脉冲相同的带宽,假设LFM 信号的脉冲宽度为T ,由匹配滤波器的压缩后,带宽就变为τ,且τT =D ≥1,这个过程就是脉冲压缩。 信号s (t ) 的匹配滤波器的时域脉冲响应为: h (t ) =s (t o -t ) t 0是使滤波器物理可实现所附加的时延。理论分析时,可令t 0=0,重写, h (t ) =s (-t ) * * 得 : h (t ) =rect ( t T ) e -j πK t 2 ?e j 2πf c t 如图, s (t ) 经过系统h (t ) 得输出信号s o (t ) s o (t ) =s (t ) *h (t ) ∞ ∞ = ? -∞∞ s (u ) h (t -u ) du = 2 ? -∞ h (u ) s (t -u ) du 2 = ? -∞ e -j πK u rect ( u T ) e j 2πf c u ?e j πK (t -u ) rect ( t -u T ) e j 2πf c (t -u ) du 当0≤t ≤T 时, T 2 s 0(t ) = ? t -T 2 e j πK t 2 e -j 2πK tu du =e = j πK t 2 e -j 2πK tu T -j 2πK t t -T 2 ?e j 2πf c t sin πK (T -t ) t πK t e j 2πf c t 当-T ≤t ≤0时, t +T 2 s 0(t ) = ? -T e j πK t 2 e -j 2πK tu du =e = j πK t 2 e -j 2πK tu t +T T 2 -j 2πK t - ?e j 2πf c t sin πK (T +t ) t πK t e j 2πf c t sin πK T (1- s 0(t ) =T t ) t rect ( t 2T πK Tt ) e j 2πf c t 即为LFM 脉冲信号经匹配滤波器得输出, 它是一固定载频f c 的信号,这是因为压 缩网络的频谱特性与发射信号频谱实现了“相位共轭匹配”,消除了色散;当t ≤T 时,包络近似为辛克(sinc )函数。 S 0(t ) =TSa (πK Tt ) rect ( t 2T ) =TSa (πBt ) rect ( t 2T ) 当πBt =±π时,t =± 1B 为其第一零点坐标;当πBt =± π 2 时,t =± 12B ,习 惯上,将此时的脉冲宽度定义为压缩脉冲宽度。 τ= 12B ?2= 1B LFM 信号的压缩前脉冲宽度T 和压缩后的脉冲宽度τ之比通常称为压缩比D D = T τ =TB ≥1 压缩比也就是LFM 信号的时宽-带宽积。 Matlab 程序: tw=50; p=10; s=dchirp(tw,p); h=conj(s(end:-1:1)); t=25e-6; N=p*tw; ts=[-N/2:N/2];ts(end)=[]; th=ts; [y,ty]=conv_m(s,ts,h,th); ty=ty/(p*tw/t); figure(3) plot(ty*1e6,abs(y)) title('The output from the matched filter to Chirp signal') xlabel('t (\mus)'),ylabel('The magnitude of y(n)') figure(4) fs=p*tw/t; y=fft(y); y=fftshift(y); N1=length(y); f=[-N1/2:N1/2]/N1*fs;f(end)=[]; plot(f/1e6,abs(y)),grid on title ('The Frequence of Chirp signal') xlabel('f (Mhz)'),ylabel('Magnitude') 仿真得时域波形 频域波形: 观察4db 带宽: 三、 数值分析 从理论来说,通过匹配滤波器以后信号的增益为B τ=25us*2Mhz=50,从仿真结果来看图一和图三 在滤波器输出端产生的信号增益为50 在输出后得到的信号脉宽理论上为1/B=1/2m=0.5us,仿真得正为0.5us 4、查看4db 脉宽 4db 应为最高电平的0.631倍,即50*0.631=31.55 得4db 脉宽为0.3us+0.2us=0.5us 由此看出,雷达测距时的距离分辨率必须和匹配滤波器输出的脉冲宽度相对应时,才能进行正常的测量,不产生距离模糊 四、 实验拓展和创新 通过加入高斯白噪声更一步真实的仿真出雷达测距在噪声背景下匹配滤波器的脉压作用: 加入如下程序段 h=conj(s(end:-1:1)); noise=sqrt(10)*randn(size(s)); s=s+noise; 仿真图像如下: 同样达到仿真目的 五、 参考文献 《基于MA TLAB 的系统分析与设计》第二版楼顺天刘小东等编著 《信号检测与估计》景占荣羊彦编著 上机练习题目(1--2) 1. 熟悉下列产生随机信号的函数(可利用matlab帮助查看该函数的使用方法): rand函数:生成均匀分布的随机数,取值从0到1。 randn函数:生成正态分布的随机数,均值为0,标准差(方差)为1。 normrnd函数:生成正态分布的随机数,均值与标准差由参数指定,函数的第 一个参数是均值,第二个参数是标准差,后面的参数用于指定 生产结果的矩阵大小。具体调用格式查看matlab帮助信息。 random函数:可以生产指定分布的随机数,调用格式查看matlab帮助信息。 2. 计算长度为N=10000的高斯随机噪声信号的均值、均方值、方差和均方差(也称标准差,即对方差开根号的值)。 提示: (1)计算均值,可以调用mean函数。调用方式为: y=mean(x) 其中,x为离散随机序列,此函数返回结果为x的均值。 (2)方差调用函数为cov函数,其调用方式为: y=cov(x) 其中,x为离散随机序列,y为其方差;当x为矩阵时,则它的每一列相当于一个变量,函数返回结果为该矩阵的列与列之间的协方差矩阵。这时diag(cov(x))是该矩阵每一个列向量的方差,sqrt(diag(cov(x)))为标准差向量。元素分别为矩阵每列元素的均值。 (3)标准差调用函数std,其调用格式为: y=std(x) y=std(x,flag) 其中,x为离散随机序列,y为其标准差,flag为控制符,用来控制计算标准差的算法。 //---------------------------------------------------------------------------------------------- N=10000; %数据长度 y=randn(1,N); %产生一个均值为0,方差为1,长度为N的随机序列 disp('平均值:'); yMean=mean(y) %计算随机序列的均值 disp('均方值:'); y2p=y*y'/N %计算其均方值,这里利用了矩阵相乘的算法 disp('均方根:'); ysq=sqrt(y2p) %计算其均方根值 disp('标准差:'); ystd=std(y,1) %计算标准差,相当于ystd=sqrt(sum((y-yMean).^2)/(N-1)) disp('方差:'); yd=ystd.*ystd 3. 求一白噪声加正弦信号以及白噪声的自相关函数,并进行分析比较。(显示出信号及相关函数的波形) 提示: MATLAB提供了计算离散随机序列相关函数的函数xcorr,其调用格式为: r=xcorr(x,y) r=xcorr(x,y,’option’) r=xcorr(x,y,maxlags,’option’) 其中,x,y为长度是N的两个独立离散随机序列,r为互相关函数的估计,’option’为控制项。该函数也可用于求随机序列x(n)的自相关函数,调用格式为:c=xcorr(x,maxlags)。具体可查看matlab帮助信息。 【部分参考答案】 clf; N=1000; Fs=500; %数据长度和采样频率 n=0:N-1; t=n/Fs; %时间序列 Lag=100; %延迟样点数 x=sin(2*pi*20*t)+0.6*randn(1,length(t)); %白噪声加正弦信号 [c,lags]=xcorr(x,Lag,'unbiased'); %估计原始信号x的无偏自相关 subplot(2,2,1), plot(t,x); xlabel('时间/s'); ylabel('x(t)'); title('带噪声周期信号'); grid on; subplot(2,2,2), plot(lags/Fs,c); %绘x信号的自相关,lags/Fs为时间序列 xlabel('时间/s'); ylabel('Rx(t)'); title('带噪声周期信号的自相关'); grid on; x1=randn(1,length(x)); %产生一与x长度一致的随机信号x1 [c,lags]=xcorr(x1,Lag,'unbiased'); %求随机信号x1的无偏自相关 subplot(2,2,3), plot(t,x1); %绘制随机信号x1 xlabel('时间/s'); ylabel('x1(t)'); title('噪声信号'); grid on; subplot(2,2,4); plot(lags/Fs,c); %绘制随机信号x1的无偏自相关 xlabel('时间/s'); ylabel('Rx1(t)'); title('噪声信号的自相关'); grid on 4. 已知两个周期信号x(t)=sin(2πft),y(t)=0.2sin(2πft+600),其中f=20Hz,求互相关函数Rxy(τ),并将这2个周期信号以及互相关的图形显示出来。 【部分参考答案】 %Samp9_4 clf; N=1000; Fs=500; %数据长度和采样频率 n=0:N-1;t=n/Fs; %时间序列 Lag=200; %最大延迟单位 x=sin(2*pi*20*t); %周期信号x y=0.2*sin(2*pi*20*t+60*pi/180); %与x有90o相移的信号y [c,lags]=xcorr(x,y,Lag,'unbiased'); %求无偏互相关 subplot(2,1,1);plot(t,x,'r'); %绘制x信号 hold on;plot(t,y,':'); %在同一幅图中绘y信号 legend('x信号', 'y信号') xlabel('时间/s');ylabel('x(t) y(t)'); title('原始信号');grid on; hold off subplot(2,1,2),plot(lags/Fs,c,'r'); %绘制x,y的互相关 xlabel('时间/s');ylabel('Rxy(t)'); title('信号x和y的相关 ');grid on 5. 在某音乐厅内,原始音频信号的回音由于墙壁和天花板等的反射而产生,听众所感受到的音频信号是x(n)和它的回音的合成。令 y(n)=x(n)+ax(n-k) 其中,x(n)=cos(0.3πn)+0.5cos(0.6πn),a=0.1,k=50。产生100个样本,求出其自相关,从中观测确定a和k。 //-------------------------------------------------------------------------------------------------- 第5题clear all clc figure(1) a=0.1; k=40; n=100:200; x(n)=cos(0.3*pi*n)+0.5*cos(0.6*pi*n); y(n)=x(n)+a*x(n-k); stem(n,y(100:200)); figure(2) t=-199:199; R=xcorr(y); stem(t,R); figure(3) %k变化(k=40:5:60),a不变(a=0.1) clear y(n) hold on k=20:20:80; d=['r' 'g' 'b' 'c' ]; for e=1:length(k) y(n)=x(n)+a*x(n-k(e)); R=xcorr(y); stem(t,R,d(e)); end figure(4) %k变化(k=0:1:60),a不变(a=0.1)的R(0) hold on clear y(n) k=0:1:60; for e=1:length(k) y(n)=x(n)+a*x(n-k(e)); R=xcorr(y); stem(k(e),max(R)); end figure(5) %a变化(a=0:0.05:0.2) k不变(k=50) clear y(n) hold on k=50; a=0:0.05:0.2; for e=1:length(a) y(n)=x(n)+a(e)*x(n-k); R=xcorr(y); stem(t,R,d(e)); end figure(6) %a变化(a=0:0.01:0.2) k不变(k=50)的R(0) hold on a=0:0.01:0.2; for e=1:length(a) y(n)=x(n)+a(e)*x(n-k); R=xcorr(y); stem(a(e),max(R)); end 6. 两个Sinc信号有0.2秒的时移(用MATLAB程序产生),用互相关函数计算 时移的大小。 【部分参考答案】 %Samp9_6 clf N=1000;n=0:N-1;Fs=500;t=n/Fs; %数据个数采样频率和时间序列 Lag=200; %最大延迟单位数 x1=90*sinc(pi*(n-0.1*Fs)); %第一个原始信号,延迟0.1s y1=50*sinc(pi*(n-0.3*Fs)); %第二个原始信号,延迟0.3s [c,lags]=xcorr(x1,y1,Lag,'unbiased'); %计算两个函数的互相关 subplot(2,1,1),plot(t,x1,'r'); %绘第一个信号 hold on;plot(t,y1,'b:'); %在同一幅图中绘第二个信号 legend('信号x', '信号y'); %绘制图例 xlabel('时间/s');ylabel('x(t) y(t)'); title('信号x和y');hold off subplot(2,1,2),plot(lags/Fs,c,'r'); %绘制互相关信号 xlabel('时间/s');ylabel('Rxy(t)'); title('信号x和y的相关 '); 范文三:随机信号处理
范文四:随机信号处理
范文五:随机信号处理