1、第七届华中地区大学生数学建模邀请赛承 诺 书我们仔细阅读了第七届华中地区大学生数学建模邀请赛的竞赛细则。我们完全明白,在竞赛开始后参赛队员不能以任何方式(包括电话、电子邮件、网上咨询等)与队外的任何人(包括指导教师)研究、讨论与赛题有关的问题。我们知道,抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料(包括网上查到的资料),必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺,严格遵守竞赛规则,以保证竞赛的公正、公平性。如有违反竞赛规则的行为,我们将受到严肃处理。我们的参赛报名号为: 10487069参赛队员 (签名) :队员1:队员2:队员3:武
2、汉工业与应用数学学会第七届华中地区大学生数学建模邀请赛组委会第七届华中地区大学生数学建模邀请赛编 号 专 用 页选择的题号: A参赛的编号: 10487069(以下内容参赛队伍不需要填写) 竞赛评阅编号:1第七届华中地区大学生数学建模邀请赛题目: 基于复合滤波和高斯拟合的声屏障加速度校正【摘 要】本文主要解决声屏障的加速度校正问题。由于声屏障检测仪内部的加速度传感器测得的数据通常会存在误差,在使用数值积分方法计算振动位移的过程中,也会累积更多的干扰,为了得到本文将从从系统误差校正、随机误差数据滤波等方面对声屏障的加速度数据进行校正。首先,基于加速度-速度和加速度-位移物理公式,利用梯形数值积分
3、公式,建立速度和位移的数值积分计算模型,将高频采样下测得的加速度数据代入仿真得到速度和位移的离散点,详见图2。并由此分别画出加速度、速度和位移的曲线图,详见图3-图11。其次,对AB段、CD段和EF段的位移曲线图做误差分析,详见图12、图13、图14。分别是从随机误差、系统误差2个角度对数据进行误差分析;发现加速度存在严重的突变和短区间内的双向变化,结合物理实际情形,运用对称性原理和物体做弹性碰撞时的平滑运动知识可知,突变和短区间双向变化均不是物体本身应具有的运动特性,属于随机误差的范畴。此外,由于随机误差的随机性,采样足够的情形下不会影响最终积分得到的最终速度,然而速度时间图像显示末速度不为
4、零,不符合物理情景。因此必然存在系统误差导致整体加速度测量的偏移。再次,基于速度和位移的数值积分计算模型和误差分析结果,建立加速度校正模型来对加速度数据进行校正,尽量消除系统误差与随机误差,使得速度和位移的计算结果基本符合物体运动事实。首先,结合ARIMA模型的思想,运用了基于拉依达准则的奇异数据滤波法、中值滤波公式和滑动平均公式相结合的复合滤波方式来尽量的消除随机误差;然后,用五点三次算法和经验筛选滤波消除系统误差;最后,运用Matlab拟合工具箱(CurveFittingTool)依次对加速度、速度和位移进行高斯拟合,并找出其与时间的函数关系。最终即可实现对测得的加速度的校正。最后,对所建
5、立的数据处理方法和模型进行评价和推广。首先分析了所建模型的优缺点,然后从马航时间联系到法航事件,本文中的数据处理方法和模型可以运用于导航测速中的累积误差的消除。关键词:数值积分 复合滤波 高斯拟合 误差分析 校正 matlab2目 录第一部分 问题重述3第二部分 问题分析3第三部分 模型假设3第四部分 符号说明3第五部分 模型的建立与算法实现41模型一:速度和位移的数值积分计算模型42判断是否有明显误差并进行误差分析93模型二:加速度校正模型11第六部分 模型评价与推广171.模型评价172.数据处理方法和模型的推广183.改进的加速度检测仪的其他应用场景18第七部分 参考文献19第八部分 附
6、录 2031、问题重述声屏障是一种控制铁路、公路、高速铁路等各种道路行车对周围环境的噪声污染有效措施之一,随着列车的大幅度加速,脉动风交替出现在列车两侧,从而引起对声屏障的拉压作用,声屏障发生摆动。正常状态下,声屏障的摆动应当在一定的范围内,当超过正常范围则需要对其进行加固维修。由于声屏障维修或重建费用高昂,故需声屏障检测仪对声屏障的工作状态进行检测,有针对性的对声屏障进行维修。声屏障检测仪的工作原理是:通过内部的加速度传感器来记录车辆经过时声屏障振动而产生的加速度数值(密集采样)。将加速度数据通过数值积分,按照加速度-位移的物理公式将加速度数据转化为震动的位移,并通过震动位移对声屏障状态进行
7、判断。在试验中,传感器测得的数据通常会存在误差,误差包括系统误差、随机误差。其中系统误差,又称为固有误差,一般其存在是具有一定的规律性,是可以被分析掌握的;随机误差,又称为测量误差,一般它的出现是不具有规律并且不可避免的。由于误差的存在,在使用数值积分方法计算振动位移的过程中,就会累积较多的干扰,故而在测得数据后,需要经过系统误差校正、随机误差数据滤波等对数据进行校正。本题中,给出了3组正常状态下的声屏障实验采样数据(见附表),所给出的加速度数据是模拟声屏障震动而用加速度检测仪所测量出的加速度数据,加速度传感器采集频率1000Hz,加速度单位为 2/g sm (g为重力加速度),数据为三组:单
8、方向从A点运动至B点;从C到D后再返回到C;从E点到F点,再由F到E,并再重复一次,其初速度皆为0。如下图1所示。图1 模拟振动示意图问题1:建立适当的数学模型,基于加速度-速度和加速度-位移物理公式,通过数值积分的方法计算声屏障的速度、位移,并基于给定数据对模型进行仿真计算,判断声屏障检测仪是否存在明显误差,从随机误差、系统误差2个角度对数据进行误差分析;问题2:基于速度和位移的数值积分计算模型和误差分析结果,建立数学模型来对加速度数据进行校正,要求能尽量消除系统误差与随机误差,使得速度和位移的计算结果基本符合物体运动事实;问题3:对你所建立的数据处理方法和模型进行推广,所改进过的加速度检测
9、仪除了可以用于声屏障监测以外,还可以应用于哪些场景,请结合改进方案阐述理由。 2、问题分析4问题一考察数值积分、随机和系统的误差分析。声屏障检测仪通过内部的加速度传感器等距采样,测得加速度数值,需要基于加速度-速度和加速度-位移物理公式,对这些离散点做数值积分,从而得到速度和位移的数值,这就是建立的数值积分计算模型。此外,在试验中,传感器测得的数据通常会存在误差,在使用数值积分方法计算振动位移的过程中,就会累积较多的干扰,故而在测得数据后,需要进行误差分析,分别为随机误差分析和系统误差分析。问题二考察系统误差校正、随机误差数据滤波,需要根据问题一中的数值积分计算模型、误差分析,对实验测得的加速
10、度进行校正,尽量消除系统误差与随机误差,使得速度和位移的计算结果基本符合物体运动事实。3、模型假设1假设不考虑重力加速度g的地区效应;2假设加速度为零的点测量值是准确的,即速度的稳定域不受误差影响;3假设所有的随机误差都属于高斯白噪声,即幅度分布服从高斯分布,功率谱密度服从均匀分布;4假设产生的系统误差对脉动风的敏感度是确定的,即系统误差是线性的。4、符号说明符号 含义 单位a 测量加速度 2/g smv 速度 /g m sx 位移 g miv残余误差 /g m s 标准差 /g m s5、模型的建立与算法实现5.1模型一:速度和位移的数值积分计算模型5.1.1模型一的分析本题要涉及到数值积分
11、。要由传感器测得的加速度数值,计算得到声屏障的速度和位移。实现的思路:根据加速度-速度和加速度-位移的物理公式,分别对加速度测试仪测得的加速度信号做数值积分,即可得到原始速度信号和原始位移信号。由于加速度传感器采集频率1000Hz,即加速度数据离散点是定步长的,而且鉴于步长较小,所以可以选择梯形公式建立数值积分模型。5.1.2模型一的建立5已知加速度-速度和加速度-位移的物理公式,如下:210 ttv v adt 210 ttx x atdt 由题意可知:0 0v 0 0x 1 0t 2 0t t所以此题中,v-a、x-a的表达式简化为:00tv adt 00tx atdt 由于a是关于时间t
12、的等步长离散点,所以不能直接做定积分,需要利用数值积分公式估算:在数学里,梯形公式是用作估算定积分 ( )ba f x dx梯形公式会把函数图像 当作成梯形并估算它的面积。以下就是估算所用的公式 ( ) ( )( ) ( ) 2ba f a f bf x dx b a 为了计算出更加准确的定积分,可以把积分的区间 分成 份,当中 趋向无限,分割出的每一个区间长度必定要是一样的,然后就可以应用梯形公式:11( ) ( )( ) ( )2 2 Nba kb a f a f b b af x dx f a kN N 亦可以写成: 0 1 2 1( ) ( ( ) 2 ( ) 2 ( ) . 2 (
13、) ( )2b N Na b af x dx f x f x f x f x f xN 当中 ,k b ax a k N for 0,1,.k N5.1.3模型一的仿真计算代入各段的加速度数据,仿真得到各段的速度和位移的离散数值。通过Matlab中应用范围广的trapz积分函数来实现。由于要逐步算出每个T的速度与位移,因此需要用到FOR循环嵌套语句。先将需要的加速度的值赋给一个临时矩阵(从第一个数据开始再逐步累加),6且建立一个储存结果的矩阵(velocity和displacement)。通过trapz函数进行积分,结果储存在对应矩阵中。经过循环加一,同样将需要的加速度的值赋给同一个临时矩阵,
14、比之上一次,将多一个加速度数据,再得到这个周期内的速度和位移值。依次循环,最终得到速度和位移的处理数据。对于此数值积分计算模型,采用matlab进行仿真计算,代码详见附录一。分别对AB段、CD段、EF段的加速度数据进行仿真得如下结果:图2.仿真结果的部分数据图3.AB段加速度7图4.AB段速度图5.AB段位移 图6.CD段加速度8图7.CD段速度图8.CD段位移图9.EF段加速度9图10.EF段速度图11.EF段位移5.2判断是否有明显误差并进行误差分析:首先判断是有明显误差的!从如下图像分析可看出:图12.AB段的位移误差分析图10在方框中的曲线应该是直接逐渐趋于平缓,而不应该有一个位移减小
15、的过程。图13.CD段的位移误差分析图(从左到右依次为方框一,方框二,方框三)方框一中的曲线近乎是斜率为恒定负值的直线,这是由于第一重积分(速度)在此时保持为0而是保持为负值导致的。正常的曲线应是斜率为0的直线,即在曲线达到第一个极大值时,不再下降。方框二中的曲线是斜率在变化的下降曲线,这就是速度变为负值,位移逐渐下降的阶段。方框三中直线原理与方框一类似。皆是由于速度未在降为0是趋于平缓而是滞后到负值时才平缓导致的。图14.EF段的位移误差分析图(从左到右依次为方框一、方框二)方框一中的图像应该在横轴上逐渐趋于平缓,方框二中的曲线应该在曲线与横轴的相交后保持斜率为0.115.2.1分析数据的随
16、机误差随机误差,又称为测量误差,一般它的出现是不具有规律并且不可避免的,在采样点足够多的情形下,为了实验模型的简便,随机误差一般视为中心值为零的高斯分布:观察AB段、CD段和EF段的加速度图像可看出:加速度检测仪所测得的加速度数据是存在很多无规律的噪声和波动的,测得的数据本身就带有误差,经过积分后,还有二重积分的累积误差。这些就是数据的随机误差。为了尽量消除这些随机误差,采用复合滤波:包括中值滤波、滑动平均滤波、拉依达准则的奇异数据滤波法。5.2.2分析数据的系统误差其中系统误差,又称为固有误差,一般其存在是具有一定的规律性,是可以被分析掌握的。传感器由于种种原因,必然存在系统误差。从AB段、
17、CD段和EF段的速度和位移曲线可以看出,在曲线的一些极点、零点和最后的那段曲线都不是正常的情况,这是由于系统的误差而导致的。系统误差的发现一般比较困难,原始数据本身一般不能揭示系统误差的存在,故只能通过一些分析途径来发现系统误差。5.3模型二:加速度校正模型5.3.1模型二的分析由5.2中的误差分析可知,所记录的加速度是既有随机误差,也有系统误差的。因此,要建立加速度校正模型,就需要消除随机误差和系统误差。1、消除随机误差:采用复合滤波:有时为了提高滤波的效果,尽量减少噪声数据对结果的影响,常将两种或两种以上的滤波算法结合在一起。对于这个问题,我们采用基于拉依达准则的奇异数据滤波法、中值滤波公
18、式和滑动平均公式相结合的复合滤波方式。2、消除系统误差:认识系统误差产生原因,重点和难点是系统非线性校正。通过一组反映被测值的离散数据,此处即为加速度检测仪所测得的加速度的采样数据,利用这些离散数据建立起一个反应被测量值变化的近似数学模型(即校正模型)。有时即使有了数学模型,例如n次多项式,但其次数过高,计算太复杂、太费时,常常要从系统的实际精度要求出发,用逼近法来降低一个已知非线性特性函数的次数,以简化数学模型,便于计算和处理。5.3.1模型二的建立1、用复合滤波法消除随机误差:复合滤波法:基于拉依达准则的奇异数据滤波法、中值滤波公式和滑动平均 2 21 e x p 22 xf x 12公式
19、相结合的复合滤波方式。步骤一:基于拉依达准则的奇异数据滤波法进行滤波拉依达准则:当测量次数N足够多且测量服从正态分布时,在各次测量值中,若某次测量值 ix所对应的剩余误差 3iv ,则认为该 ix为坏值,予以剔除。其中, iv为残余误差,为标准差拉依达准则法实施步骤:(1)求N次测量值 1x至 Nx 的算术平均值以AB为例:(2)求各项的残余误差 ivi iv x x (3)计算标准偏差(4)判断并剔除奇异项 3iv ,则认为该 ix为坏值,予以剔除。步骤二:中值滤波公式滤波中值滤波法是一种非线性平滑技术,它将每一像素点的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值。相当于y=new(
20、x),new的操作是,从在以x为中心,长度为2k的原信号中(区间为x-k+1,x+k),提取出这段区间内中间的那个值,作为y=new(x)的结果。中值滤波对孤立的噪声像素即椒盐噪声、脉冲噪声具有良好的滤波效果。由于其并不是简单的取均值,所以,它产生的模糊也就相对比较少。211 1 n ii vn 211 1 n ii vn 1n ii xx n 13图15.中值滤波的流程图步骤三:滑动平均公式滤波原理:如果滑动窗长为n的话,滑动平均就是让数据通过一个n点的FIR滤波器,滤波器抽头系数都是1,这样取滑动平均就是起到序列平滑的作用。现介绍基于filter函数的计算滑动平均:设原始数据为x,平均窗口
21、设为a(a为正整数),那么无权重滑动平均后的数据y为:2、用五点三次算法和经验筛选滤波消除系统误差(1)五点三次算法:建立了基于五点三次算法的mean5_3函数对数据进行进一步平滑处理,去除信号上的毛刺:以AB为例:共1395个数据:n=1395;y为经过处理后的数据;y=mean5_3(x,.m);(注:m为循环次数)1 11 1 1.a a ay x x xa a a 1 11y xa 2 2 11 1y x xa a 1 11 1 1.i i i i ay x x xa a a 141 2 3 1394 1395 , , ,., , ;x x x x x x a x ;(1 )for k
22、 m 1 2 4 3 51 69 4( ) 6a ;70a a a ab 1 5 2 3 42 2( ) 27a 12 8 ;35a a a ab (3 2)for n 2 2 1 13( ) 12( ) 1735j j j j jj a a a a ab 4 1 2 31 2( ) 27 12 8 ;35n n n n nn a a a a ab 1 3 2 469 4( ) 670n n n n nn a a a a ab ;;a b ;y a(2)经验筛选滤波:观察剔除了随机误差后的数据,尝试不同经验值低通滤波,来对数据进行筛选。同时进行了模糊对称化处理,再对数据进行非线性曲线拟合。3、
23、各数据拟合修正高斯拟合:运用Matlab拟合工具箱(CurveFittingTool)依次对加速度、速度和位移进行高斯拟合,并找出其与时间的函数关系。(1)对AB段速度做1阶高斯拟合:置信区间为95%时的系数为:a1= 0.1387b1= 0.5682c1= 0.06483(2)对CD段位移做8阶高斯拟合:211( )1( ) x bcf x a e 15置信区间为95%时的系数为:a1= 0.008159b1= 2.332c1= 0.06863a2= 0.01133b2= 1.721c2= 0.2714a3= 0.009678b3= 1.254c3= 0.1789a4= 0.006937b4
24、= 0.9761c4= 0.0802a5= 0.006415b5= 1.082c5= 0.1144a6= 0.004734b6= 1.465c6= 0.1525a7= 0.008172b7= 2.036c7= 0.1754a8= 0.008143b8= 2.217c8= 0.1167(3)对EF段位移做2阶高斯拟合:2 21 21 2( ) ( )1 2( ) x b x bc cf x a e a e 置信区间为95%时的系数为:a1=0.02542b1=1.54c1=0.205322 2 231 2 431 2 4( )( ) ( ) ( )1 2 3 4( ) x bx b x b x
25、bcc c cf x a e a e a e a e 2 2 2 25 6 7 85 6 7 8( ) ( ) ( ) ( )5 6 7 8x b x b x b x bc c c ca e a e a e a e 16a2=0.0195b2=0.9443c2=0.18395.3.2模型的求解与检验(检验的是CD和EF)用matlab编程,对数据进行滤波等处理后仿真得到图形如下(详细代码见附录):图16.CD速度校正前后对比图(红色为校正前,蓝色为校正后)图17.CD位移校正前后对比图(红色为校正前,蓝色为校正后)17图18.EF速度校正前后对比图(红色为校正前,蓝色为校正后)图19.EF位移
26、校正前后对比图(红色为校正前,蓝色为校正后)6、模型评价与推广6.1模型评价:6.1.1优点:(1)在数据处理上,由于对离散数据拟合来求积分得速度和位移的精确度较低,我们采用了等步长梯形积分的方法,直接对离散数据进行积分处理,得到较为精确速度和位移数据。由于数据的可信性极高,为下一步的误差分析亦提供了有力18依据。(2)在对加速度数据的随机误差进行滤波处理时,我们首先采用的是中值滤波。中值滤波对孤立的噪声像素即椒盐噪声和脉冲噪声具有良好的滤波效果。由于其并不是简单的取均值,所以,它产生的模糊也就相对比较少。数据的可信度高。其次,我们采用了滑动平均滤波,通过Matlab自带的函数,简介有力的处理
27、数据。低通滤波,将噪声剔除。其次,我们还基于五点三次算法,建立了mean5_3函数,反复检验最合适的参数值,去除毛刺。(3)在模型的建立上,我们尽量简化了无用的公式,贴近“用最简单的方法来解决问题”的数模思想。(4)在模型的图像处理和显示上,我们采用了高斯拟合,拟合效果好,置信概率大。将校正前后曲线放置在同一张图上显示,能够清晰的判别校正情况的良好程度。6.1.2缺点:(1)数据的随机误差由滤波器处理,未比较其他方法,没有全面的分析各种方法的孰优孰劣。(2)数据在校正后的图像大致符合科学运动规律,但仍存在积分误差带来的影响。(3)没有考虑声屏障加速度传感器的频率响应函数,是对加速度数据来进行分
28、析。(4)系统误差的修正是根据需求的结果进行反向推演,以追求最高可信度为原则。然而需求的结果并不是唯一确定的,只能进行预估。6.2数据处理方法和模型的推广本文应用了数值积分、复合滤波、拟合等数据处理方法,建立了数值积分计算模型和加速度误差校正模型。由于是直接对离散数据进行的积分处理,而实际中测得的数据往往都是采样之后的离散数据,所以此处建立的数值积分计算模型在实际中适用范围很广,适用于所有实际测得的定步长且步长较小的离散数据的数值积分。由于在进行加速度误差校正时,不是单纯的使用一种滤波方法,而是采用了多种滤波组合使用的复合滤波,滤波效果很好,而实际中的随机误差是不确定的,所以这种滤波方式在实际
29、中的运用也很广泛。而此处对于系统误差的分析和消除的方式则推广型不强,因为是根据加速度监测仪分析的系统误差,得出的特定的系统误差函数,进而消除的系统误差,所以在实际中其他测试仪器的推广并不适用,但是分析误差和消除误差的数据处理方法却是可以在其他的测试仪中进行推广的。6.3改进的加速度检测仪可以应用的其他场景:6.3.1场景一:由于机械加工过程中的最终精度是由机床的刀具以及相对的位移所决定的。对于刀具和工件间的相对位移则同样是通过运动学模型来计算的。由于位移最终由运动学公式积分而得,所以同样会产生累积误差。这需要运用误差剔除模型。19我们建立的误差分析模型原理与机械加工制造中产生的误差类似,可以在
30、修改相关参数后有所运用。6.3.2场景二:在天体物理学研究中,同样会在对天体测速的过程中产生引力噪音干扰,造成误差。我们的模型可以用于初步的校正这些的数据,让天体运动学计算常数的过程中不产生较大偏差。6.3.3场景三:本小组一直持续关注马航事件,而根据这次的加速度校正模型,受马航事件启发,自然联想到2009年6月1日的法航空难事件,空难事件发生的重要原因之一正是接收飞行速度的皮托管(空速管)被冻住,使得飞行速度读取错误,导致AFS(自动飞行系统)停止运行最后高空失速造成的。因此,我们认为皮托管测速并不足以保障飞行的安全。在恶劣的环境下,它亦有概率发生故障。且皮托管实测的是空速,受风速影响,而并
31、非是地速。我们知道,飞机还能够导航,可以通过一组陀螺仪精确测量自己的一切运动来推算位置、姿态、速度,优点是启动后长时间、不依赖外部来源地计算导航参数,可是缺点是长时间后误差会累积。因此校正误差至关重要。而误差累积的校正方法则与我们所建立的模型相关。在通过误差剔除校正后,让飞机得到精确的数据,可以进一步确保飞机安全,愿空难事件不再发生!7、参考文献1周小祥,陈尔奎,吕桂庆,等.基于数字积分和LMS的振动加速度信号处理J.自动化仪表,2006,27(9):51253.2姜启源,数学模型.北京:高等教育出版社,20033薛冬梅,ARIMA模型及其在时间序列分析中的应用A.吉林化工学院学报2010年6
32、月4MartinMeyer,信号处理模拟与数字信号、系统及滤波器,机械工业出版社,2011.1.15DuaneHanselman,BruceLittlefield著,朱仁峰译,精通MATLAB7,清华大学出版社,2006.1.6金祖东,机械制造中积累误差的存在与解决研究,TG506,科学之友2010年11月208、附录附录一:(问题一的源代码)%solvingdatadata=load(mydataAB.txt);%换成mydataCD.txtmydataEF.txt即为CDEF速度,位移数据。t_size=numel(data);velocity=zeros(t_size,1);displa
33、c=zeros(t_size,1);fori=1:t_sizetempdata=zeros(i,1);t=1:i.*(1/1000);forj=1:itempdata(j,1)=data(j,1);endvelocity(i,1)=trapz(t,tempdata,1);%梯形积分endfori=1:t_sizetempdata=zeros(i,1);t=1:i.*(1/1000);forj=1:itempdata(j,1)=velocity(j,1);enddisplac(i,1)=trapz(t,tempdata,1);endt=1:t_size.*(1/1000);附录二:(问题二的源代
34、码一)functiony=mean5_3(x,m)%x为被处理的数据%m为循环次数n=length(x);a=x;fork=1:mb(1)=(69*a(1)+4*(a(2)+a(4)-6*a(3)-a(5)/70;b(2)=(2*(a(1)+a(5)+27*a(2)+12*a(3)-8*a(4)/35;21forj=3:n-2b(j)=(-3*(a(j-2)+a(j+2)+12*(a(j-1)+a(j+1)+17*a(j)/35;endb(n-1)=(2*(a(n)+a(n-4)+27*a(n-1)+12*a(n-2)-8*a(n-3)/35;b(n)=(69*a(n)+4*(a(n-1)+a
35、(n-3)-6*a(n-2)-a(n-4)/70;a=b;endy=a;附录三:(问题二的源代码二)%solvingdata%去掉噪声信号data=load(mydataCD.txt);%同理,可换data_impro=mean5_3(data,5);data_impro=data_impro;data_rem=data_impro;t_size=numel(data);datanew=medfilt1(data_impro,5);%5为控制滤波区间的参数win_size=5;data_impro=filter(ones(1,win_size)/win_size,1,datanew);fori
36、=1:t_sizeifdata_impro(i)0.12data_impro(i)=0;endendifabs(data_impro(i)thr_sigma);22num_w=numel(M);fori=1:t_sizetempdata=zeros(i,1);t=1:i.*(1/1000);forj=1:itempdata(j,1)=data_impro(j,1);endvelocitynew(i,1)=trapz(t,tempdata,1);tempdata=;t=;endfori=2488:t_sizeifvelocitynew(i)0ifvelocitynew(i)0.05velocitynew(i)=0;endendendfori=1:t_sizetempdata=zeros(i,1);t=1:i.*(1/1000);forj=1:itempdata(j,1)=velocitynew(j,1);enddisplacnew(i,1)=trapz(t,tempdata,1);tempdata=;t=;endt=1:t_size.*(1/1000);