在对数字序列进行滤波的过程中,我们发现对某些序列滤波后,在首端和尾端出现异常值(主要为异常大值)。具体的原因我还不是很清楚,期待得到数学解释。
消除上述问题的一个途径是在滤波之前给序列加窗函数。窗函数当然有很多了,Hanning窗,Hamming窗,Dolph-Chebyshev窗,Gaussing窗等,各个窗函数的优点与缺点需要根据实际处理需求才便于评价。下面就Hanning窗为例,对序列加窗函数及滤波进行说明。
- 根据序列长度(采样点数),计算得到窗函数序列。假设我们有一个序列
data
,长度为N
,那么我们需要生成一个长度为N*leng_rate
的窗函数。为什么不生成长度为N的窗函数呢?因为使用一个完整的窗函数来乘以原始序列会改变原始序列的全部振幅,实际上我们的目的是在序列(或信号)的两端各加一个衰减的窗。
from scipy import signal
import numpy as np
# === PARA ====
length_rate = 0.1;
dt = 0.01;# sample rate
# ==============
mywin = np.ones((N,));
lenwinhalf = int(N*length_rate);
lenwin = int(lenwinhalf*2);
han_window = signal.hanning(lenwin);
mywin[0:lenwinhalf] = han_window[0:lenwinhalf];
mywin[N-lenwinhalf:N] = han_window[lenwinhalf:lenwin];
- 对原始序列加窗函数。
data = np.multiply(han_window,data);
- 制作Butter带通滤波器,滤波频带为拐角频率下界(
freq_lowcut
),拐角频率上界(freq_highcut
)
nyquest = 0.5*(1/dt);
digilow = freq_lowcut/nyquest;
digihigh = freq_highcut/nyquest;
[b, a] = signal.butter(4, (digilow,digihigh), btype='bandpass');
- 对原始序列进行带通滤波
data_filt = signal.filtfilt(b,a,data,axis=0);