15
$\begingroup$

I have a signal that is measured at 100Hz and I need to apply the Savitzky-Golay smoothing filter on this signal. However, on closer inspection my signal is not measured at perfectly constant rate, the delta between measurements ranges between 9.7 and 10.3 ms.

Is there a way to use the Savitzky-Golay filter on not equally spaced data? Are there other methods that I could apply?

$\endgroup$
  • $\begingroup$ A 1991 paper by Gorry is on pretty much this exact subject datapdf.com/…. But tldr, datageist's answer is the right main idea (local least-squares). What Gorry observes is that the coefficients depend on the independent variables only and are linear in the dependent variables (like Savitzky-Golay). Then he gives a way to compute them, but if you're not writing an optimized library, any old least-squares fitter could be used. $\endgroup$ – Dave Pritchard Aug 24 at 7:09
5
$\begingroup$

One method would be to resample your data so that it is equally spaced, then you can do whatever processing you like. Bandlimited resampling using linear filtering isn't going to be a good option since the data isn't uniformly spaced, so you could use some sort of local polynomial interpolation (e.g. cubic splines) to estimate what the underlying signal's values are at "exact" 10-millisecond intervals.

$\endgroup$
  • $\begingroup$ I had this solution in mind as a last resort. I'm wondering if in the end this approach gives a better solution than just assuming that my signal is measured at a constant rate. $\endgroup$ – VLC Mar 9 '12 at 15:13
  • $\begingroup$ I think even if it's non-uniformly sampled you can still use sinc() interpolation (or a a different highly sampled low pass filter). This may give better results than a spline or a pchip $\endgroup$ – Hilmar Mar 10 '12 at 16:42
  • 1
    $\begingroup$ @Hilmar: You're correct. There are a number of ways you could resample the data; approximate sinc interpolation would be the "ideal" method for bandlimited resampling. $\endgroup$ – Jason R Mar 10 '12 at 16:56
15
$\begingroup$

Because of the way the Savitzky-Golay filter is derived (i.e as local least-squares polynomial fits), there's a natural generalization to nonuniform sampling--it's just much more computationally expensive.

Savitzky-Golay Filters in General

For the standard filter, the idea is to fit a polynomial to a local set of samples [using least squares], then replace the center sample with the value of the polynomial at the center index (i.e. at 0). That means the standard SG filter coefficients can be generated by inverting a Vandermonde matrix of sample indicies. For example, to generate a local parabolic fit across five samples $y_0\dots y_4$ (with local indicies -2,-1,0,1,2), the system of design equations $Ac = y$ would be as follows:

$$ \begin{bmatrix}-2^0 & -2^1 & -2^2 \\ -1^0 & -1^1 & -1^2 \\ 0^0 & 0^1 & 0^2 \\ 1^0 & 1^1 & 1^2 \\ 2^0 & 2^1 & 2^2 \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ y_4 \end{bmatrix}. $$

In the above, the $c_0 \dots c_2$ are the unknown coefficients of the least squares polynomial $c_0 + c_1x + c_2x^2$. Since the value of the polynomial at $x = 0$ is just $c_0$, computing the pseudoinverse of the design matrix (i.e. $c = (A^TA)^{-1}A^T y\space $) will yield the SG filter coefficients in the top row. In this case, they would be

$$ \begin{bmatrix}c_0 \\ c_1 \\ c_2 \end{bmatrix} = \begin{bmatrix} -3 & 12 & 17 & 12 & -3 \\ -7 & -4 & 0 & 4 & 7 \\ 5 & -3 & -5 & -3 & 5 \\ \end{bmatrix} \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ y_4 \end{bmatrix}. $$

Note that since the derivative of $c_0 + c_1x + c_2x^2$ is $c_1 + 2c_2x$, the second row of the matrix (which evaluates $c_1$) will be a smoothed derivative filter. The same argument applies for successive rows--they give smoothed higher-order derivatives. Note that I scaled the matrix by 35 so the first row would match the smoothing coefficients given on Wikipedia (above). The derivative filters each differ by other scaling factors.

Nonuniform Sampling

When the samples are evenly spaced, the filter coefficients are translation-invariant, so the result is just an FIR filter. For nonuniform samples, the coefficients will differ based on the local sample spacing, so the design matrix will need to be constructed and inverted at each sample. If the nonuniform sample times are $x_n$, and we construct local coordinates $t_n$ with each center sample time fixed at $0$, i.e.

$$ \begin{align} t_{-2} & = x_{-2} - x_0 \\ t_{-1} & = x_{-1} - x_0 \\ t_0 & = x_0 - x_0 \\ t_1 & = x_1 - x_0 \\ t_2 & = x_2 - x_0 \end{align} $$

then each design matrix will be of the following form:

$$ A = \begin{bmatrix} t_{-2}^0 & t_{-2}^1 & t_{-2}^2 \\ t_{-1}^0 & t_{-1}^1 & t_{-1}^2 \\ t_0^0 & t_0^1 & t_0^2 \\ t_1^0 & t_1^1 & t_1^2 \\ t_2^0 & t_2^1 & t_2^2 \end{bmatrix} = \begin{bmatrix} 1 & t_{-2} & t_{-2}^2 \\ 1 & t_{-1} & t_{-1}^2 \\ 1 & 0 & 0 \\ 1 & t_1 & t_1^2 \\ 1 & t_2 & t_2^2 \end{bmatrix}. $$

The first row of the pseudoinverse of $A$ dotted with the local sample values will yield $c_0$, the smoothed value at that sample.

$\endgroup$
  • $\begingroup$ sounds like it moves from O(log(n)) to O(n^2). $\endgroup$ – EngrStudent - Reinstate Monica Feb 12 at 21:29
  • $\begingroup$ Here's an implementation of Scala described by datageist upwards. $\endgroup$ – Medium core May 9 at 2:12
  • 1
    $\begingroup$ @Mediumcore You didn't add a link to your original post. Also, I deleted it because it didn't provide an answer to the question. Please try to edit datageist's post to add a link; it'll be moderated in after review. $\endgroup$ – Peter K. May 9 at 13:41
4
$\begingroup$

"As a cheap alternative, one can simply pretend that the data points are equally spaced ...
if the change in $f$ across the full width of the $N$ point window is less than $\sqrt{N/2}$ times the measurement noise on a single point, then the cheap method can be used."
$\qquad -$ Numerical Recipes pp. 771-772

(derivation anyone ?)

("Pretend equally spaced" means:
take the nearest $\pm N/2$ points around each $t$ where you want SavGol($t$),
not snap all $t_i \to i$ . That may be obvious, but got me for a while.)

$\endgroup$
1
$\begingroup$

I found out, that there are two ways to use the savitzky-golay algorithm in Matlab. Once as a filter, and once as a smoothing function, but basically they should do the same.

  1. yy = sgolayfilt(y,k,f): Here, the values y=y(x) are assumed to be equally spaced in x.
  2. yy = smooth(x,y,span,'sgolay',degree): Here you can have x as an extra input and referring to the Matlab help x does not have to be equally spaced!
$\endgroup$
0
$\begingroup$

If it's of any help, I've made a C implementation of the method described by datageist. Free to use at your own risk.

/**
 * @brief smooth_nonuniform
 * Implements the method described in  https://dsp.stackexchange.com/questions/1676/savitzky-golay-smoothing-filter-for-not-equally-spaced-data
 * free to use at the user's risk
 * @param n the half size of the smoothing sample, e.g. n=2 for smoothing over 5 points
 * @param the degree of the local polynomial fit, e.g. deg=2 for a parabolic fit
 */
bool smooth_nonuniform(uint deg, uint n, std::vector<double>const &x, std::vector<double> const &y, std::vector<double>&ysm)
{
    if(x.size()!=y.size()) return false; // don't even try
    if(x.size()<=2*n)      return false; // not enough data to start the smoothing process
//    if(2*n+1<=deg+1)       return false; // need at least deg+1 points to make the polynomial

    int m = 2*n+1; // the size of the filter window
    int o = deg+1; // the smoothing order

    std::vector<double> A(m*o);         memset(A.data(),   0, m*o*sizeof(double));
    std::vector<double> tA(m*o);        memset(tA.data(),  0, m*o*sizeof(double));
    std::vector<double> tAA(o*o);       memset(tAA.data(), 0, o*o*sizeof(double));

    std::vector<double> t(m);           memset(t.data(),   0, m*  sizeof(double));
    std::vector<double> c(o);           memset(c.data(),   0, o*  sizeof(double));

    // do not smooth start and end data
    int sz = y.size();
    ysm.resize(sz);           memset(ysm.data(), 0,sz*sizeof(double));
    for(uint i=0; i<n; i++)
    {
        ysm[i]=y[i];
        ysm[sz-i-1] = y[sz-i-1];
    }

    // start smoothing
    for(uint i=n; i<x.size()-n; i++)
    {
        // make A and tA
        for(int j=0; j<m; j++)
        {
            t[j] = x[i+j-n] - x[i];
        }
        for(int j=0; j<m; j++)
        {
            double r = 1.0;
            for(int k=0; k<o; k++)
            {
                A[j*o+k] = r;
                tA[k*m+j] = r;
                r *= t[j];
            }
        }

        // make tA.A
        matMult(tA.data(), A.data(), tAA.data(), o, m, o);

        // make (tA.A)-¹ in place
        if (o==3)
        {
            if(!invert33(tAA.data())) return false;
        }
        else if(o==4)
        {
            if(!invert44(tAA.data())) return false;
        }
        else
        {
            if(!inverseMatrixLapack(o, tAA.data())) return false;
        }

        // make (tA.A)-¹.tA
        matMult(tAA.data(), tA.data(), A.data(), o, o, m); // re-uses memory allocated for matrix A

        // compute the polynomial's value at the center of the sample
        ysm[i] = 0.0;
        for(int j=0; j<m; j++)
        {
            ysm[i] += A[j]*y[i+j-n];
        }
    }

    std::cout << "      x       y       y_smoothed" << std::endl;
    for(uint i=0; i<x.size(); i++) std::cout << "   " << x[i] << "   " << y[i]  << "   "<< ysm[i] << std::endl;

    return true;
}

smoothing

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.