写在前面
空闲一周，突发奇想想测试一下cdo和python的带通滤波函数。我自己常用的滤波方法是cdo和python的butternworth，研究波动的这是跨不去的技术手段。最近看到一篇Bsiso的文章，发现信号提取的第一步基本上也离不开带通滤波。
当然，以上指的主要是时间滤波。空间滤波用到的人太少了这里不管他。
方法比较
这里，比较一下几个内容：
cdo
获取输入文件中每个网格点的时间序列，并对其进行快速傅里叶变换，转换到频域。根据特定算子及其参数，频域中的某些频率会被滤波（置零），然后频谱通过逆快速傅里叶变换转换回时域。频率的确定基于输入文件的时间轴。（由于转换需要假设时间增量恒定，因此数据应具有恒定的时间增量。但是，时间增量必须不为零。）所有作为参数给出的频率均按年解释。这是基于365天日历的假设。因此，如果要精确地进行多年滤波，则必须删除2月29日的数据。
如果输入文件采用360年日历，则频率参数fmin和fmax应乘以360/365的系数，才能获得精确的结果。设置频率滤波器时，必须将频率参数调整为数据中的某个频率。此处fmin向下取整，fmax始终向上取整。
使用cdo进行带通滤波注意点：
1、 为了获得可靠的结果，必须对时间序列进行去趋势处理（cdo detrend）
2、 不支持缺失值
cdo bandpass,x,y -del29feb infile outfile
with x<=0.5 and y >=0.5.
相对90%的人懒得去看说明文档，这里直接给出运行命令。
x，y分布表示滤波的时间，以20-80天为例，且数据分辨率为daily，计算方式为：
[365/80, 365/20] = [4.5625,18.25]
需要传入两个浮点参数，举个例子，同时确保删除了数据中02.29这一天：
signal.butterworth
defbutterworth_filter(data, low_period=80, high_period=30, order=5):
"""
Butterworth filter with higher order for sharper cutoff.
Parameters:
-----------
data : array-like
Input time series data
low_period : int
Lower period bound in days
high_period : int
Upper period bound in days
order : int
Filter order (higher = sharper, default=5)
Returns:
--------
filtered_data : array-like
Filtered time series
"""
from scipy import signal
low_freq = 2.0 / low_period
high_freq = 2.0 / high_period
b, a = signal.butter(order, [low_freq, high_freq], 'bandpass')
filtered = signal.filtfilt(b, a, data, axis=0, padtype='even')
return filtered
geocat
defgeocat_fourier_bandpass(data, low_period=80, high_period=20):
"""
Apply band-pass filter using geocat.comp.fourier_band_pass (frequency-domain).
>>> olr_filtered_geocat = geocat_fourier_bandpass(
... olr_anom_clean, low_period=80, high_period=20)
"""
from geocat.comp.fourier_filters import fourier_band_pass
import time as _time
fs = 1.0# daily sampling = 1 cycle/day
cutoff_freq_low = 1.0 / low_period # e.g. 1/80 cycles/day
cutoff_freq_high = 1.0 / high_period # e.g. 1/20 cycles/day
print(f"Applying geocat.comp.fourier_band_pass ({high_period}–{low_period} day)...")
print(f" cutoff_frequency_low = 1/{low_period} = {cutoff_freq_low:.5f} cycle/day")
print(f" cutoff_frequency_high = 1/{high_period} = {cutoff_freq_high:.5f} cycle/day")
filtered = fourier_band_pass(
data,
frequency=fs,
cutoff_frequency_low=cutoff_freq_low,
cutoff_frequency_high=cutoff_freq_high,
time_axis=0,
)
# Restore coordinates and attributes
filtered = filtered.assign_coords(data.coords)
filtered.attrs = data.attrs.copy()
filtered.attrs['long_name'] = (
f'{high_period}-{low_period} day fourier_band_pass (geocat)'
)
filtered.name = 'olr_anomaly_geocat_fourier'
return filtered
比较结果
cdo与signal
比较了一下signal不同滤波方法与cdo的结果，对去除气候态平均和季节日循环的olr进行20-80天的带通滤波
看起来还凑合，如果以cdo为真值的话。但是意外的是，cdo的速度运行非常之慢，这让我非常不理解。 anyway，结果我还是会相信cdo的更准确的/
所有比较结果
下面，对比了各种滤波函数的滤波结果，以及滤波速度，并与cdo的滤波结果进行比较。
总体来说，butterworth的参数order=5和3的时候效果都还不错，比lanczos更接近cdo的结果。其中，geocat的滤波结果最接近cdo的结果，相关系数最高，约为0.99797.当然，这里后面想了一下，也不能完全用cdo当作真值，考虑到cdo的频率参数是四舍五入的。
从计算速度上，butterworth-3最快，geocat也还不赖吧。2s的速度，后续再考虑进行dask加速，应该可以接受。
所以，后续进行带通滤波，还是优先考虑geocat和butterworth-3吧。
★https://code.mpimet.mpg.de/projects/cdo/embedded/index.html
https://geocat-comp.readthedocs.io/en/v2025.03.0/user_api/generated/geocat.comp.fourier_filters.fourier_band_pass.html
