Filtre de lissage Savitzky-Golay pour des données pas également espacées

16

J'ai un signal mesuré à 100 Hz et je dois appliquer le filtre de lissage Savitzky-Golay sur ce signal. Cependant, en y regardant de plus près, mon signal n'est pas mesuré à une vitesse parfaitement constante, le delta entre les mesures se situe entre 9,7 et 10,3 ms.

Existe-t-il un moyen d'utiliser le filtre Savitzky-Golay sur des données pas également espacées? Y a-t-il d'autres méthodes que je pourrais appliquer?

VLC
la source
Un article de 1991 de Gorry porte à peu près sur ce sujet exact datapdf.com/… . Mais tldr, la réponse de datageist est la bonne idée principale (moindres carrés locaux). Ce que Gorry observe, c'est que les coefficients dépendent uniquement des variables indépendantes et sont linéaires dans les variables dépendantes (comme Savitzky-Golay). Ensuite, il donne un moyen de les calculer, mais si vous n'écrivez pas une bibliothèque optimisée, n'importe quel ancien monteur de moindres carrés pourrait être utilisé.
Dave Pritchard

Réponses:

5

Une méthode consiste à rééchantillonner vos données de manière à ce qu'elles soient également espacées, puis vous pouvez effectuer le traitement que vous souhaitez. Le rééchantillonnage à bande limitée utilisant le filtrage linéaire ne sera pas une bonne option car les données ne sont pas uniformément espacées, vous pouvez donc utiliser une sorte d'interpolation polynomiale locale (par exemple des splines cubiques) pour estimer les valeurs du signal sous-jacent à "exactes" Intervalles de 10 millisecondes.

Jason R
la source
J'avais cette solution en tête en dernier recours. Je me demande si à la fin cette approche donne une meilleure solution que de simplement supposer que mon signal est mesuré à un taux constant.
VLC
Je pense que même s'il est échantillonné de manière non uniforme, vous pouvez toujours utiliser l'interpolation sinc () (ou un filtre passe-bas hautement échantillonné différent). Cela peut donner de meilleurs résultats qu'une spline ou une puce
Hilmar
1
@ Hilmar: Vous avez raison. Il existe plusieurs façons de rééchantillonner les données; une interpolation sinc approximative serait la méthode "idéale" pour un rééchantillonnage à bande limitée.
Jason R
15

En raison de la façon dont le filtre Savitzky-Golay est dérivé (c'est-à-dire sous forme d'ajustements polynomiaux des moindres carrés locaux), il y a une généralisation naturelle à l'échantillonnage non uniforme - c'est juste beaucoup plus cher en termes de calcul.

Filtres Savitzky-Golay en général

Pour le filtre standard, l'idée est d'adapter un polynôme à un ensemble local d'échantillons [en utilisant les moindres carrés], puis de remplacer l'échantillon central par la valeur du polynôme à l'indice central (c'est-à-dire à 0). Cela signifie que les coefficients de filtre SG standard peuvent être générés en inversant une matrice Vandermonde d'indices d'échantillon. Par exemple, pour générer un ajustement parabolique local sur cinq échantillons (avec les indices locaux -2, -1,0,1,2), le système d'équations de conception A c = y serait le suivant:y0y4UNEc=y

[-20-21-22-10-11-12000102101112202122][c0c1c2]=[y0y1y2y3y4].

c0c2c0+c1X+c2X2X=0c0c=(UNETUNE)-1UNETy  ) produira les coefficients du filtre SG dans la ligne supérieure. Dans ce cas, ils seraient

[c0c1c2]=[-3121712-3-sept-404sept5-3-5-35][y0y1y2y3y4].

c0+c1X+c2X2c1+2c2Xc1 ) sera un filtre dérivé lissé. Le même argument s'applique aux lignes successives - elles donnent des dérivées d'ordre supérieur lissées. Notez que j'ai redimensionné la matrice de 35 afin que la première ligne corresponde aux coefficients de lissage donnés sur Wikipedia (ci-dessus). Les filtres dérivés diffèrent chacun par d'autres facteurs d'échelle.

Échantillonnage non uniforme

Xntn0

t-2=X-2-X0t-1=X-1-X0t0=X0-X0t1=X1-X0t2=X2-X0

alors chaque matrice de conception aura la forme suivante:

UNE=[t-20t-21t-22t-10t-11t-12t00t01t02t10t11t12t20t21t22]=[1t-2t-221t-1t-121001t1t121t2t22].

UNE c0

datageist
la source
on dirait qu'il passe de O (log (n)) à O (n ^ 2).
EngrStudent
Voici une implémentation de Scala décrite par datageist vers le haut.
Noyau moyen
1
@Mediumcore Vous n'avez pas ajouté de lien à votre message d'origine. De plus, je l'ai supprimé car il ne répondait pas à la question. Veuillez essayer de modifier la publication de datageist pour ajouter un lien; il sera modéré après examen.
Peter K.
4


F sur toute la largeur duN la fenêtre de points est inférieure à N/2 fois le bruit de mesure sur un seul point, alors la méthode bon marché peut être utilisée. "
- Recettes numériques pp. 771-772

(dérivation quelqu'un?)

("Faire semblant d'être également espacés" signifie:
prendre le±N/2 points autour de chacun t où vous voulez SavGol (t),
ne pas tout cassertjeje. C'est peut-être évident, mais ça m'a pris un moment.)

denis
la source
1

J'ai découvert qu'il y a deux façons d'utiliser l'algorithme savitzky-golay dans Matlab. Une fois comme filtre, et une fois comme fonction de lissage, mais fondamentalement, ils devraient faire de même.

  1. yy = sgolayfilt (y, k, f): Ici, les valeurs y = y (x) sont supposées être également espacées en x.
  2. yy = lisse (x, y, span, 'sgolay', degré): Ici, vous pouvez avoir x comme entrée supplémentaire et en se référant à l'aide de Matlab, x n'a pas besoin d'être également espacé!
Jochen
la source
0

Si cela peut vous aider, j'ai fait une implémentation C de la méthode décrite par datageist. Libre d'utiliser à vos risques et périls.

/**
 * @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;
}

lissage

techwinder
la source