三/五/七/九点二次平滑法

三/五/七/九点二次平滑法

做时间序列处理的时候,经常需要用到平滑算法,一般的平滑法削峰严重,采用二次平滑法能够保留足够的峰值信息。这里提供一个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
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值