Smoothing Example with Savitzky-Golay Filter in Python

     Savitzky-Golay filter is used in signal processing to eliminate noise in a signal and improve the smoothness of a signal trend. The filter calculates a polynomial fit of each window based on polynomial degree and window size. 

    SciPy API provides the savgol_filter() function to implement Savitzky-Golay filter in Python. In this tutorial, we'll briefly learn how to smooth the signal data by using savgol_filter() function in Python. 

    The tutorial covers:

  1. Preparing signal data
  2. Smoothing with Savitzky-Golay filter
  3. Source code listing

    We'll start by loading the required libraries.

 
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt 
 
 
 
Preparing signal data
 
    First, we'll prepare signal data for this tutorial and visualize it in a graph.
 

x = np.arange(-2, 2, 0.02)
y = np.array(x**2 + 2*np.sin(x*np.pi)) 
y = y + np.array(np.random.random(len(x))*2.3)
 
plt.figure(figsize=(12, 9))
plt.plot(x, y, label="y_signal")
plt.legend()
plt.grid(True)
plt.show() 
 
 


    Our task is to smooth above signal by using  savgol_filter() function.
 
 
Smoothing with Savitzky-Golay filter
 
    To smooth the signal data, Savitzky-Golay filter calculates a polynomial fit of each window based on polynomial degree and window size. Scipy savgol_filter() function requires one dimensional array data, window length, polynomial order, and other optional parameters. We can set the parameters as shown below.  
 

y_smooth = signal.savgol_filter(y, window_length=11, polyorder=3, mode="nearest")
 


'mode' is an optional parameter and it determines the type of extension to use for the padded signal to which the filter is applied. Mode can be defined as 'mirror', 'nearest', 'constant', and 'wrap' type. You can check the definition of each mode and other optional parameters of the function in SciPy reference page. 

    Now, we can visualize smoothed data in a graph. 

 
plt.figure(figsize=(12, 9))
plt.plot(x, y, label="y_signal")
plt.plot(x, y_smooth, linewidth=3, label="y_smoothed")
plt.legend()
plt.grid(True)
plt.show() 
 

 


    You can change the window size and polynomial order parameters to smooth signal trend well. Here, we'll test multiple window size and check smoothing effect in a given data.  


ig, ax = plt.subplots(4, figsize=(8, 14))
i = 0
 
# define window sizes 5, 11, 21, 31
for w_size in [5, 11, 21, 31]:    
    y_fit = signal.savgol_filter(y, w_size, 3, mode="nearest")
    ax[i].plot(x, y, label="y_signal", color="green")
    ax[i].plot(x, y_fit, label="y_smoothed", color="red")
    ax[i].set_title("Window size: " + str(w_size))
    ax[i].legend()
    ax[i].grid(True)
    i+=1

plt.tight_layout()        
plt.show() 
 
 
 
    As you have seen in above graph, we can change smoothness of signal trend by changing the window size. This helps us to decide optimal window size value for our target data. 
 
    In this tutorial, we've briefly learned how to smooth signal data by using savgol_filter() function in Python. The full source code is listed below. 
 
 
Source code listing

 
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

x = np.arange(-2, 2, 0.02)
y = np.array(x**2+2*np.sin(x*np.pi)) 
y = y + np.array(np.random.random(len(x))*2.3)
 
plt.figure(figsize=(12, 9))
plt.plot(x, y, label="y_signal")
plt.legend()
plt.grid(True)
plt.show()


y_smooth = signal.savgol_filter(y, window_length=11, polyorder=3, mode="nearest")
    
plt.figure(figsize=(12, 9))
plt.plot(x, y, label="y_signal")
plt.plot(x, y_smooth, linewidth=3, label="y_smoothed")
plt.legend()
plt.grid(True)
plt.show()


ig, ax = plt.subplots(4, figsize=(8, 14))
i = 0
 
# define window sizes 5, 11, 21, 31
for w_size in [5, 11, 21, 31]:    
    y_fit = signal.savgol_filter(y, w_size, 3, mode="nearest")
    ax[i].plot(x, y, label="y_signal", color="green")
    ax[i].plot(x, y_fit, label="y_smoothed", color="red")
    ax[i].set_title("Window size: " + str(w_size))
    ax[i].legend()
    ax[i].grid(True)
    i+=1

plt.tight_layout()        
plt.show() 
   
 
 
References:

No comments:

Post a Comment