17  Savitzky–Golay

The Savitzky-Golay filter smoothes a noisy signal by performing a polynomial fit over a sliding window. It is a close relative of LOESS, which also fits local polynomials, but gives more weight to the points near the center of the window.

Polynomial fit of order 3, window size = 51 pts

Polynomial fit of order 2, window size = 51 pts

The simulations look different because the order of the polynomial makes a very different impression on us, but in reality the outcome of the two filterings is exactly the same, except near the edges of the data:

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.dates import DateFormatter
import matplotlib.dates as mdates
import datetime as dt
import matplotlib.ticker as ticker
from scipy.signal import savgol_filter
import os
import scipy
import seaborn as sns
sns.set(style="ticks", font_scale=1.5)  # white graphs, with large and legible letters
# %matplotlib widget

# dirty trick to have dates in the middle of the 24-hour period
# make minor ticks in the middle, put the labels there!
# from https://matplotlib.org/stable/gallery/ticks/centered_ticklabels.html

def centered_dates(ax):
    date_form = DateFormatter("%d %b")  # %d 3-letter-Month
    # major ticks at midnight, every day
    ax.xaxis.set_major_locator(mdates.DayLocator(interval=1))
    ax.xaxis.set_major_formatter(date_form)
    # minor ticks at noon, every day
    ax.xaxis.set_minor_locator(mdates.HourLocator(byhour=[12]))
    # erase major tick labels
    ax.xaxis.set_major_formatter(ticker.NullFormatter())
    # set minor tick labels as define above
    ax.xaxis.set_minor_formatter(date_form)
    # completely erase minor ticks, center tick labels
    for tick in ax.xaxis.get_minor_ticks():
        tick.tick1line.set_markersize(0)
        tick.tick2line.set_markersize(0)
        tick.label1.set_horizontalalignment('center')

df = pd.read_csv('shani_2022_january.csv', parse_dates=['date'], index_col='date')
start = "2022-01-02"
end = "2022-01-05"
df = df.loc[start:end]

df['sg_3_51'] = savgol_filter(df['TD'], window_length=51, polyorder=3)
df['sg_2_51'] = savgol_filter(df['TD'], window_length=51, polyorder=2)

fig, ax = plt.subplots(figsize=(8,5))

plot_data, = ax.plot(df['TD'], color='black')
plot_sg2, = ax.plot(df['sg_2_51'], color='xkcd:hot pink')
plot_sg3, = ax.plot(df['sg_3_51'], color='xkcd:mustard')

ax.legend(handles=[plot_data, plot_sg2, plot_sg3],
          labels=['data', 'sg order 2', 'sg order 3'],
          frameon=False)

plot_settings = {
    'ylim': [5, 17.5],
    'xlim': [df.index[0], df.index[-1]],
    'ylabel': "Temperature (°C)",
    'title': "Yatir Forest, 2022",
    'yticks': [5, 10, 15]
}

ax.set(**plot_settings)
centered_dates(ax)

In the middle of the data the pink curve is hidden under the mustard one: the two filters are identical there. This is a general property of the filter: orders 2 and 3 give the same smoothing, and so do orders 4 and 5, and so on. They only differ in the first and last 25 points, where the window does not fit inside the data.

To really see the difference between window width and polynomial order, we need to play with their ratio,

\text{ratio} = \frac{w}{p} = \frac{\text{window width}}{\text{polynomial order}}

start = "2022-01-02 00:00:00"
end = "2022-01-02 23:50:00"
df = df.loc[start:end]

# window_length, polyorder
df['sg_1'] = savgol_filter(df['TD'], 5, 3)
df['sg_2'] = savgol_filter(df['TD'], 11, 2)
df['sg_3'] = savgol_filter(df['TD'], 25, 3)

fig, ax = plt.subplots(figsize=(8,5))

plot_data, = ax.plot(df['TD'], color='black')
plot_sg1, = ax.plot(df['sg_1'], color='xkcd:hot pink')
plot_sg2, = ax.plot(df['sg_2'], color='xkcd:mustard')
plot_sg3, = ax.plot(df['sg_3'], color='xkcd:royal blue')

ax.legend(handles=[plot_data, plot_sg1, plot_sg2, plot_sg3],
          labels=['data', r'$w/p=1.7$', r'$w/p=5.5$', r'$w/p=8.3$'],
          frameon=False)

plot_settings = {
    'ylim': [5, 17.5],
    'xlim': [df.index[0], df.index[-1]],
    'ylabel': "Temperature (°C)",
    'title': "Yatir Forest, 2022",
    'yticks': [5, 10, 15]
}

ax.set(**plot_settings)

locator = mdates.AutoDateLocator(minticks=7, maxticks=11)
formatter = mdates.ConciseDateFormatter(locator)

ax.xaxis.set_major_locator(locator)
ax.xaxis.set_major_formatter(formatter)

The higher the ratio, the more aggressive the smoothing.

17.1 it’s a kernel too

The Savitzky-Golay filter is also a sliding weighted average, just like the ones we saw in the previous pages: for each window, the result is the sum of the data points multiplied by weights, that is, by a kernel. Let’s compare the kernel of the Savitzky-Golay filter (51 points, order 3) with a gaussian kernel of the same width:

from scipy.signal import savgol_coeffs
from scipy.signal.windows import gaussian

window_width = 51
points = np.arange(window_width) - window_width // 2   # -25, ..., 25

kernel_sg = savgol_coeffs(window_width, 3)        # the weights of the Savitzky-Golay filter
kernel_gauss = gaussian(window_width, 6)
kernel_gauss = kernel_gauss / kernel_gauss.sum()  # weights that add up to 1, like the Savitzky-Golay ones

fig, ax = plt.subplots(figsize=(8,5))
ax.plot(points, kernel_gauss, color='xkcd:royal blue', label='gaussian')
ax.plot(points, kernel_sg, color='xkcd:mustard', label='Savitzky-Golay, order 3')
ax.axhline(0, color='gray', lw=0.5)
ax.set(xlabel="points from the center of the window",
       ylabel="kernel weights",
       title="two kernels")
ax.legend(frameon=False);

Look at the tails of the Savitzky-Golay kernel: some of the weights are negative! This is how the filter manages to follow peaks and valleys much better than an ordinary running average, which flattens them.

There is a lot more about the Savitzky-Golay filter, but for our purposes this is enough. If you want some more discussion about how to choose the parameters of the filter, read this.