菜鸟进阶系列·MATLAB数学建模·数据预处理(一）剔除异常值及平滑处理
最重点的事情在本篇的最后，你也可以直接看末尾——medfilt优先级最高
数据预处理，顾名思义，预处理操作有很多，总的来说，它的一个形象的起搏就是“胃肠前端的咀嚼”，较为复杂的预处理和正式处理的界限会模糊——对于消化系统，胃其实也算肠的预处理前端，但它比咀嚼要复杂的多。
我们在咀嚼时，会把硌牙的异物吐出来已免其进入肠胃——剔除异常值，为此需要判定硌牙的程度——置信
拉依达方法就是高中数学所学的3西格玛原则，大于3西格玛原则的我们认为它不可信，也就是认为它是假的，就把它剔除。话说代码实例略是几个意思？考题不会出现？或者说，这事干脆不用MATLAB，用Excel/WPS表格就行?——但我用后者发现不会弄……还是用MATLAB自己写一个吧~
我们现在来把异常值踢了：
行列两重循环遍历所有元素，将异常值置0
这个wn叫肖伟勒系数，它随着测量次数n的增加而缓慢增加，
Sx是数据xi的标准差，拉依达方法的3西格玛原则相当于让肖伟勒系数wn=3。
x=load('error.dat');%导入数据
n=length(x);%求得测量次数n
subplot(2,1,1);
plot(x,'o');
title('原始数据')
axis([0,n+1,min(x)-1,max(x)+1]);%axis函数设置图的横纵视野范围比数据范围大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]);
与拉依达和肖维勒方法相比，下边这种方法风格差异较大：一阶差分法
看起来这个方法似乎用不太上的样子……得实时采集，还得有较强的线性，讲义上提了一嘴，然后就话锋一转：
对于平滑处理，形象的起搏是“给食物剥皮”——噪声就好比皮壳，它并不与“肠胃”亲和。
如图，我们倾向于认为一段数据是趋于稳定的，而突出的毛刺被置信为需要滤去噪声。相比于剔除异常数据，平滑处理不仅是剔除了毛刺，还填充了一个我们的期望数据：
先建立一个函数文件Y，
% 建立“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点平滑后');
这个方法对于2n+1的各个点是平权的，当然按邻域影响观，越近的应该有越高的权重，于是我们还想不平权——加权。
权重系数可以采用最小二乘原理，使平滑后的数据以最小均方差逼近原始数据。即令：
话说，这操作求了个导，比赛时能用得上么……也没给个代码示例……或者这些原理干脆就是屁话，直接调用库函数拉倒了~
比赛优先级最高的应该是smooths函数：
因为它给了个示例是读取240条股市数据，并进行平滑处理：
以下内容copy讲义：
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平滑');
运行结果：
我个人感觉高斯窗是最为我们常用的，如图它的平滑处理最贴近原始数据，这让我想起了超分辨率荧光显微技术中的STORM，用高斯拟合重建物像。
5、 用medfilt1函数（一维中值滤波）
medfilt顾名思义，med（medium）表示“中”，filt（filter）表示“过滤”，该篇文档首先介绍的是“均值滤波”的平滑处理，用局部平均数代替噪声值，现在的中值滤波就是用局部中位数代替噪声值。
调用格式：
y = medfilt1(x,n)%由此，该函数的两个参量是输入信号和选定的取中值的窗口宽度
y = medfilt1(x,n,blksz)
y = medfilt1(x,n,blksz,dim)
例： 产生一列正弦波信号，加入噪声信号，然后调用medfilt1函数对加入噪声的正弦波进行滤波（平滑处理）。
代码：
t = linspace(0,4*pi,500)'; % 产生一个从0到4*pi的向量，长度为500个点，当t轴用
y = 100*sin(t); % 产生正弦波信号
noise = normrnd(0,15,500,1); % 产生500行1列的服从N(0,15)分布的随机数，作为噪声信号
%%话说讲义上也不知咋回事，写了个“服从N(0,152)分布”，盲猜是九键键盘5和2相邻误打了
y = y + noise; % 将正弦波信号叠加一个噪声信号
subplot(2,1,1);
plot(t,y);
xlabel('时间'); ylabel('加噪声的正弦波');
yy = medfilt1(y,30); % 调用medfilt1对加噪正弦波信号y进行中值滤波、绘制窗宽为30的波形图%讲义上这个地方句子冗余，我给它改成这样了
subplot(2,1,2);
plot(t,y,'b:'); % b:表示蓝色虚线
hold on
plot(t,yy,'k','LineWidth',2); % 绘制平滑后曲线，黑色实线，线宽2
xlabel('时间'); ylabel('中值滤波'); legend('加噪波形','平滑后波形');
运行结果：
yy还是500个点，没有缺损
yy的每一个点都是局部取中位数，这个平滑作用是很强大的。
这个图看起来效果上比smooth更好，这是因为——
“中值滤波是图像处理中的一个常用步骤，它对于斑点噪声（en:speckle noise）和椒盐噪声（en:salt-and-pepper noise）来说尤其有用。保存边缘的特性使它在不希望出现边缘模糊的场合也很有用。当要求在降低噪声的同时要求保持边缘，中值滤波较卷积有更好的效果”——引用自https://blog.csdn.net/u010849228/article/details/48216919
我觉得，medfilt的优先级最高
