【不连续有噪音的时间序列：如何获得光滑连续的上包络线？(1)】

不连续有噪音的时间序列：如何获得光滑连续的上包络线？(1)
思路：
1.上包络：滑动窗口取最大值（movmax）；
2.噪音过滤：slope,Z-score,IQR等;
3.连续光滑：多项式拟合，SG/低通滤波，样条插值等。

函数代码

1. 该主函数直接粘贴复制存好，可根据自己需求删选。
2. 光有主函数还不够，主函数中还包含10个子函数代码：包括3+1=4个过滤器子函数，2个设定nan子函数，4个光滑器子函数是matlab自带的，调用即可；都是很实用的，虽然简单，但全放这里太多了，分成3个帖子分别存放，无脑存就可
3. 说明部分比较全，不同方法的filter_info,smooth_info需要输入的信息都在里面了，按照每个子函数参数说明粘贴删改就可以
4. 不同方法的效果不同，调参方法也不同，说明里讲了一部分，其他未说明的有待使用者自己摸索，欢迎在评论区指正补充
5. 已知bug：多项式拟合法在拟合阶数特别大时会报黄色警告，比如阶数取100，但是我的结果正常我就没理，但是有的时候阶数大于80，结果全部为nan，也不知道为啥，这时候可以试试别的方法，效果都挺好的，尚未发现其他bug

function out_put = smooth_upper_envelope(date,y, window_size,...
                           filter_info,smooth_info)
    %% 函数说明
    % 基本输入：
    % date-时间信息
    % y - 时间序列数据 (1D double)，包含 NaN 值
    % window_size=14; %- 包络线窗口大小

    %- 过滤器
    % filter_info.no=1; %-Slope斜率,itr_thresh: 筛选次数，thresh_idx:标准差倍数阈值 (通常为 3)
    % filter_info.no=2; %-Z-score,Z-score 的阈值 (通常为 3)
    % filter_info.no=3; %-IQR,四分位数范围（Interquartile Range, IQR），标记超出1.5倍IQR的点为离群点
   
    %- 不同过滤器所需参数及调参方法
    % filter_info=table;
    % filter_info.missing_thresh=30; % 允许最大连续缺失数据个数
    % filter_info.no=2; %- 过滤器选择
    % filter_info.itr_thresh=5; % 筛选次数，最好10次以内
    % filter_info.thresh_idx=3; % 标准差倍数阈值 (通常为 3)
                                % Z-score 的阈值 (通常为 3)
                                % IQR 的阈值 (通常为 1.5)

    %- 平滑方法
    % smooth_info.no=1; %-多项式拟合,n-拟合阶数
    % smooth_info.no=2; %-SG滤波,k为多项式阶数, f为滤波窗口长度
    % smooth_info.no=3; %-样条插值,p为控制平滑程度的参数，范围在[0,1]之间，数值越小平滑程度越高
    % smooth_info.no=4; %-低通滤波,n为滤波器阶数, cutoff_freq为截止频率    
    %- 不同方法所需参数及调参方法
    % smooth_info=table;
    % smooth_info.no=2; %- 平滑方法
    % smooth_info.n=100; % 多项式拟合阶数，越大越平滑

    % smooth_info.k=3;%k为多项式阶数，越大越贴近原数据
    % smooth_info.f=101;%f为滤波窗口长度,必须为奇数，越大越平滑

    % smooth_info.p=0.000001;%p为平滑参数, [0,1], 越小越平滑

    % smooth_info.n_low=3;% n为滤波器阶数，越大越贴近原数据
    % smooth_info.cutoff=0.015;% cutoff_freq为截止频率，越小越平滑
    
    % 输出：
    % y_smooth - 光滑连续的上包络线 (1D double)

    %% 1 前期处理
    x=(1:length(y))';
    % 计算上包络线，采用滑动窗口最大值
    y_max = movmax(y, window_size);

    %% 2 不同异常值筛选方法——关键
    filter=filter_info.no;
    missing_thresh=filter_info.missing_thresh;
    % 法一：slope过滤    
    if filter==1 
        itr_thresh=filter_info.itr_thresh;
        thresh_idx=filter_info.thresh_idx;
        y_cleaned = removeOutliersSLOPE(y_max, itr_thresh,thresh_idx);
    end
    % 法二：ZScore过滤    
    if filter==2 
        thresh_idx=filter_info.thresh_idx;
        y_cleaned = removeOutliersZScore(y_max, thresh_idx);
    end
    % 法三：IQR过滤    
    if filter==3 
        thresh_idx=filter_info.thresh_idx;
        y_cleaned = removeOutliersIQR(y_max, thresh_idx);
    end

    y_cleaned=setNaN_edge(y,y_cleaned);% 两端恢复nan
    y_cleaned= setNaN_center(y,y_cleaned,missing_thresh);% 两端恢复nan

    %% 3 不同平滑方法
    % Step 4: 平滑时间序列
    y_filled = fillmissing(y_cleaned, 'linear'); % 先插值填充NaN
    smooth=smooth_info.no;
    % 法一：多项式拟合    
    if smooth==1  
       n=smooth_info.n;    
       p = polyfit(x, y_filled, n); 
       % 选择适当的n阶多项式，阶数100能达到较好效果
       y_fit = polyval(p, x);
       y_smooth=y_fit;
    end

    % 法二：Savitzky-Golay滤波
    if smooth==2
       k=smooth_info.k;
       f=smooth_info.f;
       y_smooth = sgolayfilt(y_filled,k,f); 
                            % k为多项式阶数≥1, f为滤波窗口长度≥60，奇数     
    end

    % 法三：样条插值 
    % csaps（Cubic Smoothing Splines, 三次平滑样条）:
    % csaps 函数允许你通过一个平滑参数 p 来控制插值的平滑度。
    % 参数的范围是 [0, 1]，其中 p=0 代表线性回归，p=1 代表普通三次样条插值。
    % 适用于需要平滑曲线的场景，可以忽略数据中的一些局部变化，只保留大趋势。
    
    % spline 函数实现的是普通的三次样条插值，它会通过所有的数据点，不会有任何平滑处理。
    % 更适合那些你想要精确通过所有数据点的场景。
    % 对于含有噪声或局部异常值的数据，spline 可能会在这些位置产生不平滑的曲线。   
    if smooth==3
       p=smooth_info.p;
       spline_struct=csaps(x, y_filled, p); % 控制平滑程度的参数，范围在[0,1]之间
       % 使用 fnval 函数评估样条，生成平滑后的数据
       xx = linspace(min(x), max(x), length(y)); % 插值的点
       y_smooth = fnval(spline_struct, xx)';
    end

    % 法四：低通滤波    
    if smooth==4
        n=smooth_info.n_low;
        cutoff_freq=smooth_info.cutoff;
        [b, a] = butter(n, cutoff_freq, 'low'); % n为滤波器阶数, cutoff_freq为截止频率
        y_smooth = filtfilt(b, a, y_filled);
    end

    y_filled= setNaN_edge(y,y_filled);% 两端恢复nan
    y_filled= setNaN_center(y,y_filled,missing_thresh);% 两端恢复nan

    y_smooth= setNaN_edge(y,y_smooth);% 两端恢复nan
    y_smooth= setNaN_center(y,y_smooth,missing_thresh);% 两端恢复nan


    %% 输出
    out_put=table;
    out_put.y=y;
    out_put.y_max=y_max;
    out_put.y_smooth=y_smooth;

    %% 粗略绘制结果比较
    figure;
    subplot(2, 1, 1);
    scatter(date,y, 'b','SizeData',1);
    hold on;
    scatter(date,y_max, 'r+','SizeData',20);
    scatter(date,y_cleaned, 'go','SizeData',5);
    title('原始数据');

    subplot(2, 1, 2);
    scatter(date,y, 'b','SizeData',1);
    hold on;      
    plot(y_filled, 'g', 'LineWidth', 1.5);
    plot(y_smooth, 'k', 'LineWidth', 1.5);
    title('拟合效果');
end