Matlab学习系列012. 数据预处理1剔除异常值及平滑处理

012. 数据预处理（1）——剔除异常值及平滑处理测量数据在其采集与传输过程中，由于环境干扰或人为因素有可能造成个别数据不切合实际或丢失，这种数据称为异常值。为了恢复数据的客观真实性以便将来得到更好的分析结果，有必要先对原始数据（1）剔除异常值；

另外，无论是人工观测的数据还是由数据采集系统获取的数据，都不可避免叠加上“噪声”干扰（反映在曲线图形上就是一些“毛刺和尖峰”）。为了提高数据的质量，必须对数据进行（2）平滑处理（去噪声干扰）；

（一）剔除异常值。

注：若是有空缺值，或导入Matlab数据显示为“NaN”（非数），需要①忽略整条空缺值数据，或者②填上空缺值。

填空缺值的方法，通常有两种：A. 使用样本平均值填充；B. 使用判定树或贝叶斯分类等方法推导最可能的值填充（略）。

一、基本思想：

规定一个置信水平，确定一个置信限度，凡是超过该限度的误差，就认为它是异常值，从而予以剔除。

二、常用方法：拉依达方法、肖维勒方法、一阶差分法。

注意：这些方法都是假设数据依正态分布为前提的。

1. 拉依达方法（非等置信概率）

如果某测量值与平均值之差大于标准偏差的三倍，则予以剔除。

3x i x x S ->

其中，11n i

i x x n ==∑为样本均值，1

2

2

11()1n

x i i S x x n =?? ???

=--∑为样本的标准偏差。

注：适合大样本数据，建议测量次数≥50次。

代码实例（略）。

2. 肖维勒方法（等置信概率）

在 n 次测量结果中，如果某误差可能出现的次数小于半次时，就予以剔除。

这实质上是规定了置信概率为1-1/2n ，根据这一置信概率，可计算出肖维勒系数，也可从表中查出，当要求不很严格时，还可按下列近似公式计算：

10.4ln()n n ω=+

Tab1. 肖维勒系数表

如果某测量值与平均值之差的绝对值大于标准偏差与肖维勒系数之积，则该测量值被剔除。

n x i x 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('原始数据')

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.5067 1.4993 1.4969

y =

1.5034 1.5062 1.5034 1.5024 1.4985 1.5007 1.5067 1.4993 1.4969

3.一阶差分法（预估比较法）

用前两个测量值来预估新的测量值，然后用预估值与实际测量值比较，若大于事先给定的允许差限值，则剔除该测量值。

预估值

112()?n n n n x

x x x ---+-= 比较判别：

?n n x x

W -< 注：该方法的特点是

（1）适合于实时数据采集与处理过程；

（2）精度除了与允许误差限的大小有关外，还与前两点测量值的精确度有关；

（3）若被测物理量的变化规律不是单调递增或单调递减函数，这一方法将在函数的拐点处产生较大的误差，严重时将无法使用。

（二）数据的平滑处理

对于一组测量数据（xi ，yi ） i=1,…,n ，不要直接就想着求出的拟合多项式的线性参数，而是要先平滑处理去掉“噪声”。平滑处理在科学研究中广泛使用，它可以减少测量中统计误差带来的影响，尤其被用于无法利用多次重复测量来得到其平均值的情况和当ｙi 随ｘi 有徒然变化的那些测量段。

1. “（2n+1点）单纯移动平均”平滑滤波

取出以yi 为中心的前后各n 个数据（yi-n, …,yi-1,yi,…yi+n ）求

平均值代替yi ，即

'

1121n

i i k n

y y n +=-=+∑ 优点：方法简单，计算方便。

缺点：方法产生误差会造成信号失真；前后各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 % 主程序 clc 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. “加权移动平均”平滑滤波

加权的基本思想：

作平均的区间内中心处数据的权值最大，愈远离中心处的数据权值越小小。这样就减小了对真实信号本身的平滑作用。

权重系数可以采用最小二乘原理，使平滑后的数据以最小均方差逼近原始数据。即令

'2min ()i k i k k

y y ++-∑

通常采用“五点二次平滑” （n=5， k=-2,-1,0,1,2）

22

012222

01222

22

0122

()0()0()0i k k i k k i k k y A A k A k y A A k A k k y A A k A k k +=-+=-+=-?---=???---=???---=??∑∑∑ 五点二次平滑权重系数表：

()'0210121

3121712335y 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

'rlowess' ——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=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([-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 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([-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([-10 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) % 高斯窗方法

output = smoothts(input, ‘e’, n) % 指数法

例4.读取股市数据，对开盘价的240条数据，调用smoothts函数进行平滑处理。

代码：

x=xlsread('D:\Program

Files\MATLAB\MyWorks\gupiaoshuju01.xls'); % 读取数据文件

p0=x(1:240,1)'; % 用开盘价所在列的前240条数据

% 注意若不转置可能导致后面处理结果异常

subplot(2,2,1);

plot(p0,'k','LineWidth',1.5);

% 绘制平滑后曲线图，黑色实线，线宽1.5

xlabel('观测序号');

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 = smoothts(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);

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 = linspace(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('加噪波形','平滑后波形');

运行结果：