三/五/七/九点二次平滑法
做时间序列处理的时候,经常需要用到平滑算法,一般的平滑法削峰严重,采用二次平滑法能够保留足够的峰值信息。这里提供一个MATLAB多点[3/5/7/9]二次平滑的代码。
下面是复制的代码,也可以通过这里下载源代码。
% Name: runningmean.m
% Function: 计算滑动平均(气象上的三/五/七/九..点二次平滑)的MATLAB程序(函数)
% Author: dwcai
% Email: 1510336267@qq.com
% Createtime: 20190724
% Usage: [rm] = runningmean(X, n, [endpoint])
% Version: V0.1
% Parameters: rm: 滑动平均序列
% X: 数据序列
% n: 平滑窗口长度(3, 5, 7, 9)
% endpoint: 处理端点['ave','drop'],ave算术平均[默认],drop去掉端点值
% Reference:
% https://wenku.baidu.com/view/556d418f3086bceb19e8b8f67c1cfad6185fe953.html
% Test code:
%--------------以下为测试代码(复制到此函数相同路径的 m 脚本中)------------------
% ts = randi(10, 50, 1);
% r3 = runningmean(ts, 3); % 默认端点值处理为算数平均法
% r5 = runningmean(ts, 5, 'drop'); % 设置为抛除端点值
% r7 = runningmean(ts, 7, 'ave');
% r9 = runningmean(ts, 9);
% plot(ts)
% hold on
% plot(r3)
% plot(r5)
% plot(r7)
% plot(r9)
% hold off
%------------------------------测试结束-----------------------------------
function [rm] = runningmean(X, n, endpoint)
if nargin < 2
error('输入参数不足!');
elseif nargin == 2
endpoint = 'ave'; % 默认首尾采用算术平均法
end
if min(size(X)) > 1
error('不能处理维度>=2的数据序列!');
end
if max(size(X)) <= n
error('输入数据个数必须大于滑动点数!');
end
if isodd(n) == false
error('滑动点数n必须为奇数[3,5,7,9].');
else
half = (n - 1) / 2; % 求半区间
end
rm = ones(size(X))*nan;
ln = length(X);
% 计算1+half:ln-half区间的 n 点二次滑动平均值
startp = 1 + half;
endp = ln - half;
switch n
case 3 % 三点二次平滑
for i = startp:endp
rm(i) = (X(i-1)+2*X(i)+X(i+1))/4;
end
case 5 % 五点二次平滑
for i = startp:endp
rm(i) = (-3*X(i-2)+12*X(i-1)+17*X(i)+12*X(i+1)-3*X(i+2))/35;
end
case 7 % 七点二次平滑
for i = startp:endp
rm(i) = (-2*X(i-3)+3*X(i-2)+6*X(i-1)+7*X(i)+6*X(i+1)+...
3*X(i+2)-2*X(i+3))/21;
end
case 9 % 九点二次平滑
for i = startp:endp
rm(i) = (-21*X(i-4)+14*X(i-3)+39*X(i-2)+54*X(i-1)+...
59*X(i)+54*X(i+1)+39*X(i+2)+14*X(i+3)-21*X(i+4))/231;
end
%case 11... 待扩展
end
% 处理边界值
switch endpoint
case 'ave' % 算术均值平滑
% 前段端点
for i = 1:startp
rm(i) = mean(X(1:i+1));
end
% 后端端点
for i = endp:ln
rm(i) = mean(X(i-1:end));
end
case 'drop' % 抛弃端点值, 则返回值的首尾出现 NaN 值
warning('端点值被抛弃, 返回值的首尾以NaN填充.');
disp('端点值被抛弃, 返回值的首尾以NaN填充.');
end
end
% 奇数判断函数
function [odd] = isodd(n)
if rem(n,2) == 1
odd = true;
else
odd = false;
end
end