滤波(Butterworth滤波器)

Function

分析气候变化趋势中常用的方法还有滤波技术，主要分为低通滤波，高通滤波以及带通滤波。Scipy库中提供了构造Butterworth滤波器的函数。

1
scipy.signal.butter(N, Wn, btype='low', analog=False, output='ba')

参数：

　　　　N: 滤波器阶数

　　　　Wn: 归一化的截止频率。计算公式为Wn=2*截止频率/采样频率。（注意：采样频率要大于两倍的信号本身最大的频率，才能还原信号。截止频率一定小于信号本身最大的频率，所以Wn一定在0和1之间）。当构造带通滤波器或者带阻滤波器时，Wn为长度为2的列表。

　　　　btype:{‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’}参数可选

　　　　　　‘lowpass’　默认值，低通滤波

　　　　　　‘highpass’　高通滤波

　　　　　　‘bandpass’　带通滤波

　　　　　　‘bandstop’　阻通滤波

analog：若为True，返回模拟滤波器，否则为数字滤波器

默认设置下，构造滤波器函数会返回滤波器的分子系数向量b和滤波器的分母系数向量a。

滤波器构建成功后，即可对数据进行滤波。滤波函数：

1
scipy.signal.filtfilt(b, a, x, axis=-1, padtype='odd', padlen=None, method='pad', irlen=None)

参数：

　　　　b: 滤波器的分子系数向量

　　　　a: 滤波器的分母系数向量

　　　　x: 要过滤的数据数组

　　　　x: 指定要过滤的数据数组的x维

　　　　padtype: {‘odd’, ‘even’, ‘constant’, None} 参数可选，默认为’odd’

　　　　padlen：整形，在应用滤波器之前在轴两端延伸的元素数。此值必须小于要滤波元素个数-1

　　　　method：确定处理信号边缘的方法，{‘pad’ ，’gust’}参数可选

　　　　　　‘pad’　填充信号；填充类型padtype和padlen决定，irlen被忽略

　　　　　　‘gust’　使用古斯塔夫森方法，而忽略padtype和padlen

　　　　irlen：当method为’gust’时，irlen指定滤波器的响应长度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
from scipy import signal
import numpy as np
y = np.array([6.08, 4.56, 5.63, 5.31, 5.15, 5.44, 4.65, 4.24, 7.3,  5.86, 4.51, 6.28, 5.55, 5.35,
              5.12, 4.76, 4.35, 3.76, 4.74, 5.55, 4.54, 5.74, 5.54, 3.67, 4.77, 4.9,  3.06, 3.9,
              4.18, 5.44, 5.21, 3.86, 3.96, 4.47, 4.37, 4.86, 4.43, 3.63, 3.98, 3.94, 5.09, 4.48,
              4.05, 4.81, 4.07, 4.48, 4.46, 3.95, 5.24, 3.54, 3.11, 5.07, 6.09, 4.59, 4.55, 4.7,
              3.43, 4.37, 4.79, 3.64, 4.3,  3.5 ])
####高通滤波####
b, a = signal.butter(3, 2/3, 'highpass')
high_series = signal.filtfilt(b, a, y)
plt.plot(np.arange(0,62,1),y,'k',label='y')
plt.plot(np.arange(0,62,1),high_series+y.mean(),'r',label='highpass')
plt.legend()
plt.show()

image-20200519214619977

1
2
3
4
5
6
7
####低通滤波####获得10年以上的信号
b, a = signal.butter(3, 2/10, 'lowpass')
low_series = signal.filtfilt(b, a, y)
plt.plot(np.arange(0,62,1),y,'k',label='y')
plt.plot(np.arange(0,62,1),low_series,'r',label='lowpass')
plt.legend()
plt.show()

image-20200519214647083

1
2
3
4
5
6
7
####带通滤波####获得3-10年间的信号
b, a = signal.butter(3, [2/10,2/3], 'bandpass')
low_series = signal.filtfilt(b, a, y)
plt.plot(np.arange(0,62,1),y,'k',label='y')
plt.plot(np.arange(0,62,1),low_series+y.mean(),'r',label='bandpass')
plt.legend()
plt.show()

image-20200519214659491