时间滤波

time_filter.png

本示例滤波频率范围分布使用了频率与周期两种参数表示，两者是完全等价的。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
"""
一维滤波演示
"""

import numpy as np
from matplotlib import pyplot as plt
import moisten_ew as mew
import xarray as xr


# 加载数据
data = xr.open_dataset("~/Documents/equator_waves/data/olr.day.mean.nc").olr.sel(
    lat=5, lon=170, time="2000"
)

# 过滤 3 ~ 10 日周期的信号
butter_filtered = mew.butter_filter(
    data, 'bandpass',
    freq=(1/3 * mew.unit('cpd'), 1/10 * mew.unit('cpd')), # 指定频率范围
    axis='time'
)

fft_filtered = mew.fft_time_filter(
    data, 'bandpass',
    period=(3 * mew.unit('day'), 10 * mew.unit('day')), # 指定周期范围
    axis='time'
)

# 画图
fig = plt.figure(figsize=(6, 6), constrained_layout=True)
gs = fig.add_gridspec(3, 1, hspace=0.1)
for i, (d, title) in enumerate(zip(
    [data, butter_filtered, fft_filtered],
    ['Original Data', 'Butterworth Filtered', 'FFT Filtered']
)):
    ax = fig.add_subplot(gs[i, 0])
    im = ax.plot(d.time, d, color='b')
    ax.set_title(title, loc='left')

fig.savefig("docs/images/time_filter.png", dpi=100)