Filter penghalusan Savitzky-Golay untuk data yang tidak berjarak sama

Saya memiliki sinyal yang diukur pada 100Hz dan saya perlu menerapkan filter smoothing Savitzky-Golay pada sinyal ini. Namun, pada pemeriksaan lebih dekat sinyal saya tidak diukur pada laju konstan sempurna, delta antara pengukuran berkisar antara 9,7 dan 10,3 ms.

Apakah ada cara untuk menggunakan filter Savitzky-Golay pada data yang tidak berjarak sama? Apakah ada metode lain yang bisa saya terapkan?

Salah satu metode adalah dengan menguji ulang data Anda sehingga spasi sama, maka Anda dapat melakukan pemrosesan apa pun yang Anda suka. Resampling tanpa batas menggunakan filter linear tidak akan menjadi pilihan yang baik karena data tidak diberi spasi secara seragam, jadi Anda bisa menggunakan semacam interpolasi polinomial lokal (misalnya splines kubik) untuk memperkirakan berapa nilai sinyal yang mendasarinya "tepat" Interval 10 milidetik.

Karena cara Savitzky-Golay filter diturunkan (yaitu sebagai cocok polinomial kuadrat-lokal), ada generalisasi alami untuk pengambilan sampel tidak seragam - itu jauh lebih mahal secara komputasi.

Savitzky-Golay Filter secara Umum

Untuk filter standar, idenya adalah untuk menyesuaikan polinomial ke set sampel lokal [menggunakan kuadrat terkecil], lalu ganti sampel pusat dengan nilai polinomial pada indeks tengah (yaitu pada 0). Itu berarti koefisien filter SG standar dapat dihasilkan dengan membalikkan matriks Vandermonde indeks sampel. Misalnya, untuk menghasilkan parabola fit lokal di lima sampel (dengan indicies lokal -2, -1,0,1,2), sistem persamaan desain A c = y akan menjadi sebagai berikut:$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}.$$

Dalam contoh di atas, adalah koefisien yang tidak diketahui dari polinomial kuadrat terkecil c 0 + c 1 x + c 2 x 2 . Karena nilai polinomial pada x = 0 hanya c 0 , menghitung pseudoinverse dari matriks desain (yaitu c = ( A T A ) - 1 A T y ) akan menghasilkan koefisien filter SG di baris atas. Dalam hal ini, mereka akan menjadi$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}.$$

Perhatikan bahwa karena turunan dari adalah c 1 + 2 c 2 x , baris kedua dari matriks (yang mengevaluasi c 1 ) akan menjadi filter derivatif yang dihaluskan. Argumen yang sama berlaku untuk baris berturut-turut - mereka memberikan turunan tingkat tinggi yang dihaluskan. Perhatikan bahwa saya menskalakan matriks dengan 35 sehingga baris pertama akan cocok dengan koefisien smoothing yang diberikan di Wikipedia (di atas). Filter derivatif masing-masing berbeda oleh faktor penskalaan lainnya.$c_0 + c_1x + c_2x^2$$c_1 + 2c_2x$$c_1$

Pengambilan sampel yang tidak seragam

Ketika sampel ditempatkan secara merata, koefisien filter adalah terjemahan-invarian, sehingga hasilnya hanyalah filter FIR. Untuk sampel tidak seragam, koefisien akan berbeda berdasarkan jarak sampel lokal, sehingga matriks desain perlu dibangun dan dibalik di setiap sampel. Jika waktu sampel tidak seragam adalah , dan kami membangun koordinat lokal t n dengan setiap waktu sampel pusat ditetapkan pada 0 , yaitu$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}$$

maka setiap matriks desain akan berbentuk sebagai berikut:

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

"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.)

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!

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