Matlab笔记——数据预处理——剔除异常值及平滑处理012

012.数据预处理（1）剔除异常值及平滑处理测量数据在其采集与传输过程中，由于环境干扰或人为因素有可能造成个别数据不切合实际或丢失，这种数据称为异常值。为了恢复数据的客观真实性以便将来得到更好的分析结果，有必要先对原始数据（1）剔除异常值；另外，无论是人工观测的数据还是由数据采集系统获取的数据，都不可避免叠加上“噪声”干扰（反映在曲线图形上就是一些“毛刺和尖峰”）。为了提高数据的质量，必须对数据进行（2）平滑处理（去噪声干扰）；（一）剔除异常值。注：若是有空缺值，或导入Matlab数据显示为“NaN”（非数）,需要忽略整条空缺值数据，或者填上空缺值。填空缺值的方法，通常有两种：A.使用样本平均值填充；B.使用判定树或贝叶斯分类等方法推导最可能的值填充（略）。一、基本思想：规定一个置信水平，确定一个置信限度，凡是超过该限度的误差，就认为它是异常值，从而予以剔除。二、常用方法：拉依达方法、肖维勒方法、一阶差分法。注意：这些方法都是假设数据依正态分布为前提的1. 拉依达方法（非等置信概率）如果某测量值与平均值之差大于标准偏差的三倍，则予以剔除。x一x3Six为（xx）21为样本的标准偏其中，x，1邑x为样本均值，S，nixi，i差。注：适合大样本数据，建议测量次数$50次。代码实例（略）。2. 肖维勒方法（等置信概率）在n次测量结果中，如果某误差可能出现的次数小于半次时，就予以剔除。这实质上是规定了置信概率为1-1/2n,根据这一置信概率，可计算出肖维勒系数，也可从表中查出，当要求不很严格时，还可按下列近似公式计算：，1<0.41n（n）nTab1.肖维勒系数表n34567891011123n1.381.531.651.731.801.861.921.962.002.03n131415203040501002005003n2.072.102.132.242.392.492.582.813.023.20如果某测量值与平均值之差的绝对值大于标准偏差与肖维勒系数之积，则该测量值被剔除。x一x,Sinx例1.利用肖维勒方法对下列数据的异常值(2.5000)进行剔除：1.50341.50621.50341.50241.49852.50001.50071.50671.49931.4969上述数据保存于文件erro.dat代码：x=load('error.dat');n=length(x);subplot(2,1,1);plot(x,'o');title('原始数据')axis(0,n+1,min(x)-1,max(x)+1);w=1+0.4*log(n);yichang=abs(x-mean(x)>w*std(x);%若用拉依达方法，把w改成3即可，但本组数据将不能成功剔除异常值。x(yichang)=;saveerrornew.datx-ASCIIsubplot(2,1,2);plot(x,'rs');title('异常值剔除后数据');axis(0,n+1,min(x)-1,max(x)+1);运彳丁结果：x=1.50341.50621.50671.49931.49691.50341.50241.49852.50001.5007y=1.50341.50621.50341.50241.49851.50071.50671.49931.49693. 一阶差分法（预估比较法）用前两个测量值来预估新的测量值，然后用预估值与实际测量值比较，若大于事先给定的允许差限值，则剔除该测量值。预估值X，x+（xx）nn1n1n2比较判别：xx<Wnn注：该方法的特点是（1）适合于实时数据采集与处理过程；（2）精度除了与允许误差限的大小有关外，还与前两点测量值的精确度有关；（3）若被测物理量的变化规律不是单调递增或单调递减函数，这一方法将在函数的拐点处产生较大的误差，严重时将无法使用。（二）数据的平滑处理对于一组测量数据（xi,yi）i=1,.,n，不要直接就想着求出的拟合多项式的线性参数，而是要先平滑处理去掉“噪声”。平滑处理在科学研究中广泛使用，它可以减少测量中统计误差带来的影响，尤其被用于无法利用多次重复测量来得到其平均值的情况和当yi随xi有徒然变化的那些测量段。1.“（2n+1点）单纯移动平均”平滑滤波取出以yi为中心的前后各n个数据（yi-n,.,yi-1,yi,.yi+n）求平均值代替yi，即叮2ni工爲k=一n优点：方法简单，计算方便。缺点：方法产生误差会造成信号失真；前后各n个数据无法平滑。适用性：适用于变化缓慢的数据。注：n越大平滑效果越好，但失真也越大。例2.“9点单纯移动平均”平滑滤波代码：%建立“n点单纯移动平均”的滤波函数%注意函数要单独保存为与函数名同名的.m文件functionY=smooth_data(y,n)m=length(y);j=1;fori=(n-1)/2+1:(m-(n-1)/2)p=i-(n-1)/2;q=i+(n-1)/2;Y(j)=sum(y(p:q)/n;j=j+1;endend%主程序clccleart=-15:0.5:15;n=length(t);Y=5./(1+t.A2);%原始测试数据y=Y+(0.5-rand(1,n);%给测试数据加上噪声干扰y1=smooth_data(y,9);%调用函数作9点滤波处理plot(1:n,Y,1:n,y,'-o',5:n-4,y1,'-*');legend('无噪声','含噪声','9点平滑后');运行结果：IIIII0102030405060702.“加权移动平均”平滑滤波加权的基本思想：作平均的区间内中心处数据的权值最大，愈远离中心处的数据权值越小小。这样就减小了对真实信号本身的平滑作用。权重系数可以采用最小二乘原理，使平滑后的数据以最小均方差逼近原始数据。即令min,(y'-y)2ikikk通常采用“五点二次平滑”(n=5,k=-2,-1,0,1,2)工(y-A-Ak-Ak2)二0ik012k=2工(yAAk-Ak2)k=0ik012k=，2工(yAAkAk2)k2=0ik0127=，2五点二次平滑权重系数表：归一系数y-2y-1y。y】丫2y-235319-3-53y-1'35913126-5y。35-3121712-3y】35-5612139y;353-5-3931y0=35("3y-2+12y，i17y0+12T3y2)3.用“smooth函数”平滑滤波调用格式：Z=smooth(Y,span,method)说明：Z：平滑后的数据向量Y：被平滑的数据向量span：平滑点数，缺省为5点method：平滑方法，缺省为移动平滑，其它还有moving'Movingaverage(default)单纯移动平均lowess'Lowess(linearfit)线性加权平滑loess'Loess(quadraticfit)二次加权平滑'sgolay'Savitzky-Golay'rlowess'RobustLowess(linearfit)'rloess'RobustLoess(quadraticfit)例3.用matlab自带的平滑函数作平滑滤波实例。代码：t=-10:0.5:10;n=length(t);y=5./(1+t.A2);%原始测试数据y1=y+0.5*(0.5-rand(1,n);%给测试数据加上噪声干扰%调用多个滤波函数作滤波处理y2=smooth(y1,3);y3=smooth(y1,9);y4=smooth(y1,3,'lowess');y5=smooth(y1,9,'lowess');y6=smooth(y1,3,'loess');y7=smooth(y1,9,'loess');y8=smooth(y1,3,'rloess');y9=smooth(y1,9,'rloess');figure(1);%第一张图subplot(3,2,1);plot(t,y);axis(-1010-16);gridontitle('无噪声信号');subplot(3,2,2);plot(t,y1,'-*');axis(-1010-16);gridontitle('含噪声信号');subplot(3,2,3);plot(t,y2,'-*');axis(-1010-16);gridontitle('3点单纯移动平均');subplot(3,2,4);plot(t,y3,'-*');axis(-1010-16);gridontitle('9点单纯移动平均');subplot(3,2,5);plot(t,y4,'-*');axis(-1010-16);gridontitle('3点线性加权平滑');subplot(3,2,6);plot(t,y5,'-*');axis(-1010-16);gridontitle('9点线性加权平滑');figure(2);%第二张图subplot(3,2,1);plot(t,y);axis(-1010-16);gridontitle('无噪声信号');subplot(3,2,2);plot(t,y1,'-*');axis(-1010-16);gridontitle('含噪声信号');subplot(3,2,3);plot(t,y6,'-*');axis(-1010-16);gridontitle('3点二次加权平滑');subplot(3,2,4);plot(t,y7,'-*');axis(-1010-16);gridontitle('9点二次加权平滑');subplot(3,2,5);plot(t,y8,'-*');axis(-1010-16);gridontitle('3点rloess平滑');subplot(3,2,6);plot(t,y9,'-*');axis(-1010-16);gridontitle('9点rloess平滑');运行结果：Figure14.用“smoothts函数”（盒子法、高斯窗法、指数法）平滑滤波调用格式：output=smoothts(input)output=smoothts(input,b',wsize)%盒子法output=smoothts(input,g',wsize,stdev)%高斯窗方法output=smoothts(input,e',n)%指数法例4.读取股市数据，对开盘价的240条数据，调用smoothts函数进行平滑处理。代码：x=xlsread('D:ProgramFilesMATLABMyWorksgupiaoshuju01.xls');%读取数据文件p0=x(1:240,1)'%用开盘价所在列的前240条数据%注意若不转置可能导致后面处理结果异常subplot(2,2,1);plot(p0,'k','LineWidth',1.5);%绘制平滑后曲线图，黑色实线，线宽1.5xlabel('观测序号');ylabel('股市日开盘价');axis(025010001400);p1=smoothts(p0,'b',30);%用盒子法平滑数据，窗宽为30subplot(2,2,2);plot(p0,'.');%绘制日开盘价散点图plot(p0,'.','markersize',3);可以改变点的大小holdonplot(p1,'k','LineWidth',1.5);xlabel('观测序号');ylabel('盒子法');legend('原始散点','平滑曲线'，'location','northwest');axis(025010001400);p2=smoothts(p0,'g',30);%高斯窗方法，窗宽为30，标准差为默认值0.65subplot(2,2,3);plot(p0,'.');holdonplot(p2,'k','LineWidth',1.5);xlabel('观测序号');ylabel('高斯窗方法');legend('原始散点','平滑曲线','location','northwest');axis(025010001400);p3=smoothts(p0,'e',30);%用指数法平滑数据，窗宽为30subplot(2,2,4);plot(p0,'.');holdonplot(p3,'k','LineWidth',1.5);xlabel('观测序号');ylabel('指数方法');legend('原始散点','平滑曲线','location','northwest');axis(025010001400);gridontitle('9点rloess平滑');运行结果：ooo543r)03Joooo5D5D2211JIJIJITJIOOH05054332ooooID5D211JIJIJIXo5DHI1DD0111105010015D20025D观捌序号05010015020D250观测序号匣姑釈点1DD01)04HIU04di原始散点平滑曲线200250100150观测序号也-Kss衆05010015020D250观测序号*原始骸点平滑曲线503oOO53200250JI-1100W501D005.用medfilt1函数(一维中值滤波)调用格式：y=medfilt1(x,n)y=medfilt1(x,n,blksz)y=medfilt1(x,n,blksz,dim)例5.产生一列正弦波信号，加入噪声信号，然后调用medfilt1函数对加入噪声的正弦波进行滤波(平滑处理)。代码：t=linspace(0,4*pi,500)'%产生一个从0到4*pi的向量，长度为500y=100*sin(t);%产生正弦波信号noise=normrnd(0,15,500,1);%产生500行1列的服从N(0,152)分布的随机数，作为噪声信号y=y+noise;%将正弦波信号加入噪声信号subplot(2,1,1);plot(t,y);xlabel('时间');ylabel('加噪声的正弦波');%调用medfiltl对加噪正弦波信号y进行中值滤波，并绘制波形图yy=medfiltl(y,30);%指定窗宽为30，对y进行中值滤波subplot(2,1,2);plot(t,y,'b:');%b:表示蓝色虚线holdonplot(t,yy,'k','LineWidth',2);%绘制平滑后曲线，黑色实线，线宽2xlabel('时间');ylabel('中值滤波');legend('加噪波形','平滑后波形');运行结果：15041M加噪液形平滑后皱形nr101214时间