Search code examples
pythonnumpysignal-processinggaussiansmoothing

how to smooth a curve in python


I have an entropy curve (1d numpy array) but this curve has a lot of noise. I would like to delete the noise with a smoothing.

This is the plot of my curve: curve and noise

I have tried to solve this issue making a convolution product with a Kaiser-Bessel filter:

gaussian_curve = window_kaiser(windowLength, beta=20)  # kaiser filter
gaussian_curve = gaussian_curve / sum(gaussian_curve)

for i in range(0, windows_number):
     start = (i * step) + 1
     end = (i * step) + windowLength
     convolution[i] = (np.convolve(entropy[start:end + 1], gaussian_curve, mode='valid'))
     entropy[i] = convolution[i][0]

but this code returns this error:

File "/usr/lib/python2.7/dist-packages/numpy/core/numeric.py", line 822, in convolve
    raise ValueError('v cannot be empty')
ValueError: v cannot be empty

the numpy.convolve operator with 'valid' mode, returns the central element in the overlap but, in this case, returns an empty element.

is there a simple way to apply a smoothing?

thanks!


Solution

  • ok, I solved. I have used another approach: Savitzky-Golay filter

    The code:

    def savitzky_golay(y, window_size, order, deriv=0, rate=1):
    
        import numpy as np
        from math import factorial
    
        try:
            window_size = np.abs(np.int(window_size))
            order = np.abs(np.int(order))
        except ValueError, msg:
            raise ValueError("window_size and order have to be of type int")
        if window_size % 2 != 1 or window_size < 1:
            raise TypeError("window_size size must be a positive odd number")
        if window_size < order + 2:
            raise TypeError("window_size is too small for the polynomials order")
        order_range = range(order+1)
        half_window = (window_size -1) // 2
        # precompute coefficients
        b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
        m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
        # pad the signal at the extremes with
        # values taken from the signal itself
        firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
        lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
        y = np.concatenate((firstvals, y, lastvals))
        return np.convolve( m[::-1], y, mode='valid')
    

    now, I can type:

    entropy = np.array(entropy)
    entropy = savitzky_golay(entropy, 51, 3) # window size 51, polynomial order 3
    

    the result is this:

    enter image description here