1、数据预处理 剔除异常值 及平滑处理 测量数据在其 采集 与 传输 过程中， 由于环境干扰 或人为因素 有可能造成 个别 数据不切合实际 或丢失 ，这种数据称为 异常 值 。 为了 恢复数据的客观真实性 以便将来得到更好的分析结果，有必要先对原始数据 剔除异常值。 另外， 无论是 人工观测的数据还是由数据采集系统获取的数据，都不可避免叠加上 “ 噪声 ”干扰 （反映在曲线图形上就是一些“毛刺和尖峰 ”） 。为了提高数据的质量，必须 对数据进行 平滑处理（去噪声干扰）。 （一） 剔除异常值。 注：若是有空缺 值， 或导入 Matlab 数据显示为“ NaN”（非数），需要 忽略整条 空缺值数据，或
2、者 填上空缺值。 填空缺值的方法，通常有两种： A. 使用 样本 平均值填充； B. 使用判定树或贝叶斯分类等方法推导最可能的值填充（略）。 一、 基本思想： 规定一个置信水平，确定一个置信限度，凡是超过 该 限度的误差，就认为它是 异常 值，从而予以剔除。 二、 常用方法：拉依达方法、肖维勒方法、一阶差分法。 注意：这些方法都是假设数据依正态分布为前提的。 1. 拉依达方法（非等置信概率 ） 如果某测量值与平均值之差大于标准偏差的三倍，则予以剔除。 3 xix x S 其中，11 n iixxn 为样本均值 ， 12211 ()1 nx iiS x xn 为样本的标准偏差。 注： 适合大样本
3、数据，建议测量次数 50 次 。 代码实例（略）。 2. 肖维勒方法（等置信概率） 在 n 次测量结果中，如果某误差可能出现的次数小于半次时，就予以剔除。 这实质上是规定了置信概率为 1-1/2n，根据这一置信概率，可计算出肖维勒系数，也可从表中查出，当要求不很严 格时，还可按下列近似公式计算： 1 0.4 ln( )n n Tab1. 肖维勒系数 表 n 3 4 5 6 7 8 9 10 11 12 n 1.38 1.53 1.65 1.73 1.80 1.86 1.92 1.96 2.00 2.03 n 13 14 15 20 30 40 50 100 200 500 n 2.07 2.1
4、0 2.13 2.24 2.39 2.49 2.58 2.81 3.02 3.20 如果某测量值与平均值之差的绝对值大于标准偏差与肖维勒系数之积，则该测量值被剔除。 nxix x S 例 1. 利用肖维勒方法对下列数据的 异常 值 （ 2.5000） 进行剔除 ： 1.5034 1.5062 1.5034 1.5024 1.4985 2.5000 1.5007 1.5067 1.4993 1.4969 上 述 数据保存于文件 erro.dat 代码： x=load(error.dat); n=length(x); subplot(2,1,1); plot(x,o); title(原始数据 )
5、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)=; save errornew.dat x -ASCII subplot(2,1,2); plot(x,rs); title(异常值剔除后数据 ); axis(0,n+1,min(x)-1,max(x)+1); 运行结果： x = 1.5034 1.5062 1.5034 1.5024 1.4985 2.5000 1.5007 1.50
6、67 1.4993 1.4969 y = 1.5034 1.5062 1.5034 1.5024 1.4985 1.5007 1.5067 1.4993 1.4969 3. 一阶差分法（预估比较法） 用前两个测量值来预估新的测量值，然后用预估值与实际测量值比较， 若大于 事先给定 的允许 差限值， 则剔除该测量值 。 预估值 1 1 2()n n n nx x x x 比较判别： nnx x W 注： 该方法的特点是 （ 1） 适合于实时数据采集与处理过程 ； （ 2） 精度除了与 允许 误差 限 的大小有关外，还与前两点测量值的精确度有关 ； （ 3）若 被测物理量的变化规律不是单调递增或单
7、调递减函数，这一方法将在函数的拐点处产生较大的误差，严重时将无法使用 。 （二）数据的平滑处理 对于一组测量数据（ xi， yi） i=1,n，不要直接就想着 求出的拟合多项式的线性参数， 而是要先 平滑处理 去掉“噪声” 。平滑处理在科学研究中广泛使用，它可以减少测量中统计误差带来的影响，尤其被用于无法利用多次重复测量来得到其平均值的情况和当 i 随 i有徒然变化的那些测量段。 1. “ （ 2n+1 点） 单纯移动平均 ”平滑滤波 取出以 yi 为中心的前后各 n 个数据（ yi-n, ,yi-1,yi,yi+n） 求平均值代替 yi，即 1121 ni iknyyn 优点：方法简单，计算
8、方便。 缺点：方法 产生 误差 会造成 信号失真 ； 前后各 n 个数据无法平滑。 适用性：适用于变化缓慢的数据。 注： n 越大 平滑效果越好，但失真也越大 。 例 2. “ 9 点 单纯移动平均 ”平滑滤波 代码： % 建立“ n点 单纯移动平均 ”的滤波函数 % 注意函数要单独保存为与函数名同名的 .m文件 function Y=smooth_data(y,n) m=length(y); j=1; for i=(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; end end % 主程序 cl
9、c clear t=-15:0.5:15; n=length(t); Y=5./(1+t.2); % 原始测试数据 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点平滑后 ); 运行结果： 2. “加权移动平均”平滑滤波 加权的基本思想： 作 平均 的 区间内中心处数据的权值最大，愈 远 离中心处的数据权值 越小 小。这样就减小了对真实信号本身的平滑作用。 权重系数可以采用最 小二乘原理， 使平滑后
10、的数据以最小均方差逼近原始数据。 即 令 2m in ( )i k i kk yy通常采用“五点二次平滑” （ n=5， k=-2,-1,0,1,2） 220 1 22220 1 222220 1 22( ) 0( ) 0( ) 0ikkikkikky A A k A ky A A k A k ky A A k A k k 五点二次平滑权重系数表 ： 归一系数 y-2 y-1 y0 y1 y2 y-2 35 31 9 -3 -5 3 y-1 35 9 13 12 6 -5 y0 35 -3 12 17 12 -3 y1 35 -5 6 12 13 9 y2 35 3 -5 -3 9 31 0
11、2 1 0 1 21 3 1 2 1 7 1 2 335y y y y y y 3. 用 “ smooth 函数 ” 平滑滤波 调用格式： Z = smooth(Y, span, method) 说明： Z： 平滑后的数据向量 Y： 被平滑的数据向量 span： 平滑点数，缺省为 5 点 method ：平滑方法，缺省为移动平滑 ，其它还有 moving Moving average (default)单纯移动平均 lowess Lowess (linear fit)线性加权平滑 loess Loess (quadratic fit)二次加权平滑 sgolay Savitzky-Golay r
12、lowess Robust Lowess (linear fit) rloess Robust Loess (quadratic fit) 例 3. 用 matlab 自带的平滑函数作平滑滤波 实例。 代码： t=-10:0.5:10; n=length(t); y=5./(1+t.2); % 原始测试数据 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=smoot
13、h(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(-10 10 -1 6); grid on title(无噪声信号 ); subplot(3,2,2); plot(t,y1,-*); axis(-10 10 -1 6); grid on title(含噪声信号 ); subplot(3,2,3); plot(t,y2,-*); axis(-10 10 -1 6); grid
14、 on title(3点单纯移动平均 ); subplot(3,2,4); plot(t,y3,-*); axis(-10 10 -1 6); grid on title(9点单纯移动平均 ); subplot(3,2,5); plot(t,y4,-*); axis(-10 10 -1 6); grid on title(3点线性加权平滑 ); subplot(3,2,6); plot(t,y5,-*); axis(-10 10 -1 6); grid on title(9点线性加权平滑 ); figure(2); % 第二张图 subplot(3,2,1); plot(t,y); axis(
15、-10 10 -1 6); grid on title(无噪声信号 ); subplot(3,2,2); plot(t,y1,-*); axis(-10 10 -1 6); grid on title(含噪声信号 ); subplot(3,2,3); plot(t,y6,-*); axis(-10 10 -1 6); grid on title(3点二次加权平滑 ); subplot(3,2,4); plot(t,y7,-*); axis(-10 10 -1 6); grid on title(9点二次加权平滑 ); subplot(3,2,5); plot(t,y8,-*); axis(-1
16、0 10 -1 6); grid on title(3点 rloess平滑 ); subplot(3,2,6); plot(t,y9,-*); axis(-10 10 -1 6); grid on title(9点 rloess平滑 ); 运行结果： Figure 1 Figure 2 4. 用 “ smoothts 函数 ” （盒子法、高斯窗法、指数法） 平滑滤波 调用格式： output = smoothts(input) output = smoothts(input, b , wsize) % 盒子法 output = smoothts(input, g , wsize, stdev)
17、 % 高斯窗方法 output = smoothts(input, e , n) % 指数法 例 4. 读取 股市数据， 对开盘价的 240 条数据， 调用 smoothts 函数进行平滑处理 。 代码： x=xlsread(D:Program FilesMATLABMyWorksgupiaoshuju01.xls); % 读取数据文件 p0=x(1:240,1); % 用开盘价所在列的前 240条数据 % 注意若不转置可能导致后面处理结果异常 subplot(2,2,1); plot(p0,k,LineWidth,1.5); % 绘制平滑后曲线图，黑色实线，线宽 1.5 xlabel(观测序
18、号 ); ylabel(股市日开盘价 ); axis(0 250 1000 1400); p1 = smoothts(p0,b,30); % 用盒子法平滑数据，窗宽为 30 subplot(2,2,2); plot(p0,.); % 绘制日 开 盘价散点图 plot(p0,.,markersize,3); 可以改变点的大小 hold on plot(p1,k,LineWidth,1.5); xlabel(观测序号 ); ylabel(盒子法 ); legend(原始散点 ,平滑曲线 ,location,northwest); axis(0 250 1000 1400); p2 = smooth
19、ts(p0,g,30); % 高斯窗方法，窗宽为 30，标准差为默认值 0.65 subplot(2,2,3); plot(p0,.); hold on plot(p2,k,LineWidth,1.5); xlabel(观测序号 ); ylabel(高斯窗方法 ); legend(原始散点 ,平滑曲线 ,location,northwest); axis(0 250 1000 1400); p3 = smoothts(p0,e,30); % 用指数法平滑数据，窗宽为 30 subplot(2,2,4); plot(p0,.); hold on plot(p3,k,LineWidth,1.5);
20、 xlabel(观测序号 ); ylabel(指数方法 ); legend(原始散点 ,平滑曲线 ,location,northwest); axis(0 250 1000 1400);grid on title(9点 rloess平滑 ); 运行结果： 5. 用 medfilt1 函数 （ 一维中值滤波） 调用格式： y = medfilt1(x,n) y = medfilt1(x,n,blksz) y = medfilt1(x,n,blksz,dim) 例 5. 产生一列正弦波信号，加入噪声信号，然后调用 medfilt1 函数对加入噪声的正弦波进行滤波（平滑处理） 。 代码： t = l
21、inspace(0,4*pi,500); % 产生一个从 0到 4*pi的向量，长度为 500 y = 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(加噪声的正弦波 ); % 调用 medfilt1对加噪正弦波信号 y进行中值滤波，并绘制波形图 yy = medfilt1(y,30); % 指定窗宽为 30，对 y进行中值滤波 subplot(2,1,2); plot(t,y,b:); % b:表示蓝色虚线 hold on plot(t,yy,k,LineWidth,2); % 绘制平滑后曲线，黑色实线，线宽 2 xlabel(时间 ); ylabel(中值滤波 ); legend(加噪波形 ,平滑后波形 ); 运行结果：
