三维网格化气候数组的25小时滤波/移动平均实现问询

针对你的n×m×p三维数组（纬度×经度×30分钟间隔时间序列），以下是几种实用的Matlab解决方案，替代downsample_ts和解决movmean的使用问题：

方法1：直接使用movmean指定时间维度 #

movmean完全支持多维数组，只需明确指定操作的**第三维（时间维度）**即可，之前的问题大概率是未指定维度导致默认处理了非时间维度：

% 计算25小时对应的时间步长：30分钟间隔 → 25×2=50个时间点
window_size = 25 * 2;

% 对第三维（时间）做滑动平均，端点处理可选'discard'/'fill'/'wrap'等
filtered_data = movmean(your_data, window_size, 3, 'Endpoints', 'discard');

• 'Endpoints'参数根据需求选择：'discard'会丢弃窗口不完整的边缘数据；'fill'用边缘值填充；'wrap'循环填充。

方法2：分格点循环处理（低内存场景） #

如果数组过大导致直接用movmean内存不足，可以逐经纬度格点单独处理时间序列：

[n, m, p] = size(your_data);
window_size = 50;
% 初始化结果数组（这里用discard端点的情况，长度为p - window_size + 1）
filtered_data = zeros(n, m, p - window_size + 1);

for lat_idx = 1:n
    for lon_idx = 1:m
        % 提取单个格点的时间序列
        ts = squeeze(your_data(lat_idx, lon_idx, :));
        % 对时间序列做滑动平均
        filtered_data(lat_idx, lon_idx, :) = movmean(ts, window_size, 'Endpoints', 'discard');
    end
end

方法3：Reshape降维优化效率 #

通过reshape将三维数组转为二维（经纬度合并为一维），处理后再恢复形状，比双循环效率更高：

window_size = 50;
% 将n×m×p转为(n*m)×p的二维数组，时间维度保留为第二维
data_2d = reshape(your_data, [], p);
% 对第二维（时间）做滑动平均
filtered_2d = movmean(data_2d, window_size, 2, 'Endpoints', 'discard');
% 恢复为三维数组
filtered_data = reshape(filtered_2d, n, m, size(filtered_2d, 2));

进阶：精准潮汐去除（滤波法） #

如果滑动平均的低通效果不够精准，可使用零相位低通滤波直接滤除潮汐频段（半日潮12h、全日潮24h均属于25小时周期以内的信号）：

fs = 2/60; % 采样频率：2次/小时（30分钟间隔）
cutoff_period = 25; % 截止周期25小时
cutoff_freq = 1/cutoff_period; % 转换为截止频率（Hz）

% 设计4阶Butterworth低通滤波器
[b, a] = butter(4, cutoff_freq/(fs/2), 'low');

[n, m, p] = size(your_data);
filtered_data = zeros(n, m, p);

for lat_idx = 1:n
    for lon_idx = 1:m
        ts = squeeze(your_data(lat_idx, lon_idx, :));
        % filtfilt实现零相位滤波，避免时间序列相位偏移
        filtered_data(lat_idx, lon_idx, :) = filtfilt(b, a, ts);
    end
end

内容的提问来源于stack exchange，提问作者lisse