matlab数据平滑 smooth函数的python实现

数据平滑是数据处理中的一种简单且常用的操作，matlab中比较好用的一个函数是smooth函数，

其调用方式为

smoothed_y = smooth(y, window_length, method)

其中y为待处理的信号，window_length为平滑作用的窗口长度，method为平滑的方法，可选方法有如下六种

 其中默认的平滑方式为滑动平均，假设window_length为w，则平滑后的值为每w个的算术平均值，很简单也很直观。

不过遗憾的是python中似乎没有这个函数，在scipy中有函数savgol_filter可以实现平滑滤波的目的，但是其原理跟matlab中的‘moving’方法不同。

这里先说一下python中savgol_filter函数的用法

from scipy.signal import savgol_filter

smoothed_y = savgol_filter(y, window_len, order_num, mode='nearest')

其中y为待处理信号，window_len为窗口大小，order_num为滤波器阶数，mode为滤波方式，科工选择的滤波方式有如下四种：

上图中以输入[1,2,3,4,5,6,7,8]为例，说明了四种不同的滤波方式所对应的添加辅助数据的方式。

在matlab的smooth函数中，我们可以令参数method='sgolay',实现这种滤波方式。

也就是说matlab中的smooth函数可以向下兼容python中的savgol_filter函数，但是savgol_filter函数却无法实现简单的滑动平均。

于是就自己动手写了一个，虽然简陋，但是经过实验验证，与matlab中的输出结果是一致的。

代码如下：

def smooth(signal, window_len=5):
#     # matlab smooth function implement for moving average method in python
#     in matlab smooth function, smoothed signal yy is calculated by the following rules (window_len=5):
#     yy(1) = y(1)
#     yy(2) = (y(1) + y(2) + y(3))/3
#     yy(3) = (y(1) + y(2) + y(3) + y(4) + y(5))/5
#     yy(4) = (y(2) + y(3) + y(4) + y(5) + y(6))/5
#     ....
#     yy(n) = y(n)

    # signal: input signal, one dimensional, 一维数组
    # window_len: dimension of the smoothing window, 窗口长度, 
    # must be an odd integer, default value: 5
    # If you specify window_len as an even number, window_len is automatically reduced by 1.
    # output:
        # the smoothed signal

    if not window_len % 2:
        if window_len == 0:
            raise Exception("Window length needs to be a positive odd number.")
        else:
            # 窗口宽度需为奇数，若为偶数则默认减1
            window_len = window_len - 1
    if window_len > len(signal):
        raise Exception("Input signal length needs to be bigger than window size.")

    half_win = window_len //2
    sig_len = len(signal)
    smoothed_signal = np.zeros(sig_len)
    start_critical_index = half_win
    end_critical_index = sig_len - half_win - 1
    count = 0
    for index in range(0, sig_len):
        if index < start_critical_index:
            smoothed_signal[index] = np.mean(signal[0:2*index+1])
        elif index > end_critical_index:
            count += 1
            smoothed_signal[index] = np.mean(signal[-window_len+2*count:])
        elif index >=start_critical_index and index <= end_critical_index:
            startInd = max([0, index-half_win])
            endInd = min(index+half_win+1, sig_len)
            smoothed_signal[index] = np.mean(signal[startInd:endInd])

    return smoothed_signal 

 也欢迎python大神批评指正，指出优化方向。