Filtro de suavização Savitzky-Golay para dados não igualmente espaçados

16

Eu tenho um sinal medido em 100Hz e preciso aplicar o filtro de suavização Savitzky-Golay neste sinal. No entanto, em uma inspeção mais minuciosa, meu sinal não é medido a uma taxa perfeitamente constante, o delta entre as medições varia entre 9,7 e 10,3 ms.

Existe uma maneira de usar o filtro Savitzky-Golay em dados não igualmente espaçados? Existem outros métodos que eu poderia aplicar?

VLC
fonte
Um artigo de 1991 de Gorry aborda praticamente esse assunto exato datapdf.com/… . Porém, a resposta do datagrama é a idéia principal correta (mínimos quadrados locais). O que Gorry observa é que os coeficientes dependem apenas das variáveis ​​independentes e são lineares nas variáveis ​​dependentes (como Savitzky-Golay). Depois, ele fornece uma maneira de calculá-las, mas se você não estiver escrevendo uma biblioteca otimizada, qualquer instalador antigo de mínimos quadrados poderá ser usado.
Dave Pritchard

Respostas:

5

Um método seria reamostrar seus dados para que sejam igualmente espaçados, para que você possa fazer o processamento que quiser. A reamostragem ilimitada de banda usando filtragem linear não será uma boa opção, pois os dados não são uniformemente espaçados, portanto, você pode usar algum tipo de interpolação polinomial local (por exemplo, splines cúbicos) para estimar quais são os valores "exatos" do sinal subjacente Intervalos de 10 milissegundos.

Jason R
fonte
Eu tinha essa solução em mente como último recurso. Gostaria de saber se, no final, essa abordagem fornece uma solução melhor do que apenas supor que meu sinal seja medido a uma taxa constante.
VLC
Acho que, mesmo que seja amostrada de maneira não uniforme, você ainda pode usar a interpolação sinc () (ou um filtro passa-baixa diferente e altamente amostrado). Isso pode dar melhores resultados do que um spline ou um pchip
Hilmar
11
@ Hilmar: Você está correto. Existem várias maneiras pelas quais você pode reamostrar os dados; uma interpolação sinc aproximada seria o método "ideal" para reamostragem ilimitada de banda.
Jason R
15

Devido à maneira como o filtro Savitzky-Golay é derivado (ou seja, como o polinômio local dos mínimos quadrados se encaixa), existe uma generalização natural para a amostragem não uniforme - é apenas muito mais caro em termos computacionais.

Filtros Savitzky-Golay em geral

Para o filtro padrão, a idéia é ajustar um polinômio a um conjunto local de amostras [usando mínimos quadrados] e substituir a amostra central pelo valor do polinômio no índice central (ou seja, 0). Isso significa que os coeficientes de filtro SG padrão podem ser gerados invertendo uma matriz de Vandermonde de indicações de amostra. Por exemplo, para gerar um ajuste parabólico local em cinco amostras (com indicações locais -2, -1,0,1,2), o sistema de equações de design A c = yy0 0...y4UMAc=y seria o seguinte:

[-20 0-21 1-22-1 10 0-1 11 1-1 120 00 00 01 10 021 10 01 11 11 1220 021 122][c0 0c1 1c2]=[y0 0y1 1y2y3y4].

Acima, são os coeficientes desconhecidos do polinômio de mínimos quadrados c 0 + c 1 x + c 2 x 2 . Como o valor do polinômio em x = 0 é apenas c 0 , o cálculo do pseudo - inverso da matriz de projeto (ou seja, c = ( A T A ) - 1 A T y ) produzirá os coeficientes do filtro SG na linha superior. Nesse caso, eles seriamc0 0...c2c0 0+c1 1x+c2x2x=0 0c0 0c=(UMATUMA)-1 1UMATy 

[c0 0c1 1c2]=[-3121712-3-7-40 0475-3-5-35][y0 0y1 1y2y3y4].

Observe que, como a derivada de é c 1 + 2 c 2 x , a segunda linha da matriz (que avalia cc0 0+c1 1x+c2x2c1 1+2c2x ) será um filtro de derivada suavizada. O mesmo argumento se aplica a linhas sucessivas - elas fornecem derivadas de ordem superior suavizadas. Observe que dimensionei a matriz em 35 para que a primeira linha correspondesse aos coeficientes de suavização fornecidos na Wikipedia (acima). Os filtros derivativos diferem por outros fatores de escala.c1 1

Amostragem não uniforme

Quando as amostras são espaçadas uniformemente, os coeficientes do filtro são invariantes à conversão, portanto, o resultado é apenas um filtro FIR. Para amostras não uniformes, os coeficientes diferem com base no espaçamento local da amostra, portanto, a matriz de projeto precisará ser construída e invertida em cada amostra. Se os tempos de amostra não uniformes são , e construímos coordenadas locais t n com cada tempo de amostra central fixo em 0 , ou seja,xntn0 0

t-2=x-2-x0 0t-1 1=x-1 1-x0 0t0 0=x0 0-x0 0t1 1=x1 1-x0 0t2=x2-x0 0

cada matriz de design terá a seguinte forma:

UMA=[t-20 0t-21 1t-22t-1 10 0t-1 11 1t-1 12t0 00 0t0 01 1t0 02t1 10 0t1 11 1t1 12t20 0t21 1t22]=[1 1t-2t-221 1t-1 1t-1 121 10 00 01 1t1 1t1 121 1t2t22].

A primeira linha do pseudoinverso de UMA pontilhado com os valores da amostra local produzirác0 0, o valor suavizado nessa amostra.

datagrama
fonte
Parece que ele se move de O (log (n)) para O (n ^ 2).
EngrStudent - Restabelece Monica
Aqui está uma implementação do Scala descrita por datageist para cima.
Núcleo médio
11
@Mediumcore Você não adicionou um link à sua postagem original. Além disso, excluí-o porque não forneceu uma resposta para a pergunta. Por favor, tente editar a publicação do datageist para adicionar um link; será moderado após a revisão.
Peter K.
4

"Como uma alternativa barata, pode-se simplesmente fingir que os pontos de dados estão igualmente espaçados ...
se a mudança nof em toda a largura do N janela de ponto é menor que N/2 vezes o ruído de medição em um único ponto, o método barato pode ser usado ".
- Receitas Numéricas pp. 771-772

(derivação alguém?)

("Finja igualmente espaçado" significa:
pegue a±N/2 pontos em torno de cada t onde você quer SavGol (t),
não encaixe tudotEuEu. Isso pode ser óbvio, mas me pegou por um tempo.)

denis
fonte
1

Eu descobri que existem duas maneiras de usar o algoritmo savitzky-golay no Matlab. Uma vez como filtro e outra como função de suavização, mas basicamente eles devem fazer o mesmo.

  1. yy = sgolayfilt (y, k, f): Aqui, os valores y = y (x) são assumidos como igualmente espaçados em x.
  2. yy = suave (x, y, span, 'sgolay', degree): Aqui você pode ter x como uma entrada extra e, consultando a ajuda do Matlab, x não precisa ser igualmente espaçado!
Jochen
fonte
0

Se for de alguma ajuda, fiz uma implementação em C do método descrito pelo datageist. Livre para usar por sua conta e risco.

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

alisamento

techwinder
fonte