Savitzky-Golay-Glättungsfilter für nicht gleichmäßig verteilte Daten

Ich habe ein Signal, das bei 100 Hz gemessen wird, und ich muss das Savitzky-Golay-Glättungsfilter auf dieses Signal anwenden. Bei näherer Betrachtung wird mein Signal jedoch nicht mit einer vollkommen konstanten Rate gemessen, das Delta zwischen den Messungen liegt zwischen 9,7 und 10,3 ms.

Gibt es eine Möglichkeit, den Savitzky-Golay-Filter für nicht gleichmäßig verteilte Daten zu verwenden? Gibt es andere Methoden, die ich anwenden könnte?

Eine Methode wäre, Ihre Daten so neu abzutasten, dass sie gleichmäßig verteilt sind, und dann können Sie jede beliebige Verarbeitung durchführen. Das bandbegrenzte Resampling mit linearer Filterung ist keine gute Option, da die Daten nicht gleichmäßig verteilt sind. Sie können daher eine lokale Polynominterpolation (z. B. kubische Splines) verwenden, um die Werte des zugrunde liegenden Signals auf "genau" zu schätzen. 10-Millisekunden-Intervalle.

Aufgrund der Art und Weise, wie das Savitzky-Golay-Filter abgeleitet wird (dh als lokale Polynomanpassungen der kleinsten Quadrate), gibt es eine natürliche Verallgemeinerung der ungleichmäßigen Abtastung - es ist nur viel rechenintensiver.

Savitzky-Golay-Filter im Allgemeinen

Für den Standardfilter besteht die Idee darin, ein Polynom an eine lokale Menge von Abtastwerten anzupassen [unter Verwendung der kleinsten Quadrate] und dann das mittlere Abtastwert durch den Wert des Polynoms am mittleren Index (dh bei 0) zu ersetzen. Das heißt, die Standard-SG-Filterkoeffizienten können durch Invertieren einer Vandermonde-Matrix von Abtastwertangaben erzeugt werden. Um beispielsweise eine lokale parabolische Anpassung über fünf Stichproben (mit lokalen Angaben -2, -1, 0 , 1, 2) zu erzeugen, wäre das System der Entwurfsgleichungen A c = y wie folgt:$y_0\dots y_4$$Ac = y$

$$\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}.$$

Oben sind die die unbekannten Koeffizienten des Polynoms der kleinsten Quadrate c 0 + c 1 x + c 2 x 2 . Da der Wert des Polynoms bei x = 0 gerade c 0 ist , ergibt die Berechnung der Pseudoinverse der Entwurfsmatrix (dh c = ( A T A ) - 1 A T y ) die SG-Filterkoeffizienten in der oberen Reihe. In diesem Fall wären sie$c_0 \dots c_2$$c_0 + c_1x + c_2x^2$$x = 0$$c_0$$c = (A^TA)^{-1}A^T y\space$

$$\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}.$$

$c_0 + c_1x + c_2x^2$$c_1 + 2c_2x$$c_1$

Ungleichmäßige Probenahme

$x_n$$t_n$$0$

$$\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}$$

dann hat jede Entwurfsmatrix die folgende 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}.$$

$A$ $c_0$

$f$$N$$\sqrt{N/2}$
$\qquad -$

(Ableitung jemand?)

$\pm N/2$$t$$t$
$t_i \to i$

Ich habe herausgefunden, dass es zwei Möglichkeiten gibt, den Savitzky-Golay-Algorithmus in Matlab zu verwenden. Einmal als Filter und einmal als Glättungsfunktion, aber im Grunde sollten sie dasselbe tun.

1. yy = sgolayfilt (y, k, f): Hier wird angenommen, dass die Werte y = y (x) in x gleich beabstandet sind.
2. yy = smooth (x, y, span, 'sgolay', degree): Hier können Sie x als zusätzliche Eingabe haben und in Bezug auf die Matlab-Hilfe muss x nicht gleichmässig verteilt sein!

Wenn es hilfreich ist, habe ich eine C-Implementierung der von datageist beschriebenen Methode durchgeführt. Die Benutzung erfolgt auf eigenes Risiko.

/**
 * @brief smooth_nonuniform
 * Implements the method described in  /signals/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;
}