Сглаживающий фильтр Савицкого-Голея для неравномерно распределенных данных

16

У меня есть сигнал, который измеряется при 100 Гц, и мне нужно применить сглаживающий фильтр Савицкого-Голея к этому сигналу. Однако при ближайшем рассмотрении мой сигнал не измеряется с совершенно постоянной частотой, разница между измерениями составляет от 9,7 до 10,3 мс.

Есть ли способ использовать фильтр Савицкого-Голея для неравномерно распределенных данных? Есть ли другие методы, которые я мог бы применить?

VLC
источник
А 1991 статья Gorry на довольно много точной теме datapdf.com/... . Но tldr, ответ datageist - правильная основная идея (локальные наименьшие квадраты). Что замечает Горри, так это то, что коэффициенты зависят только от независимых переменных и являются линейными по зависимым переменным (как Савицкий-Голей). Затем он дает способ вычислить их, но если вы не пишете оптимизированную библиотеку, можно использовать любой старый инструмент для наименьших квадратов.
Дэйв Притчард

Ответы:

5

Одним из методов будет повторная выборка ваших данных, чтобы они были равномерно распределены, тогда вы можете выполнять любую обработку по своему усмотрению. Ограниченная полоса передискретизации с использованием линейной фильтрации не будет хорошим вариантом, так как данные не распределены равномерно, поэтому вы можете использовать какую-то локальную полиномиальную интерполяцию (например, кубические сплайны), чтобы оценить, какие значения базового сигнала имеют «точное» 10-миллисекундные интервалы.

Джейсон Р
источник
Я имел в виду это решение в качестве последнего средства. Мне интересно, если в итоге этот подход даст лучшее решение, чем просто предположить, что мой сигнал измеряется с постоянной скоростью.
VLC
Я думаю, что даже если он неравномерно дискретизирован, вы все равно можете использовать интерполяцию sinc () (или другой высокочастотный фильтр низких частот). Это может дать лучшие результаты, чем сплайн или чип
Хильмар
1
@Hilmar: Ты прав. Существует несколько способов повторной выборки данных; Приблизительная интерполяция sinc была бы «идеальным» методом для ограниченного диапазона дискретизации.
Джейсон Р
15

Из-за того, как получен фильтр Савицкого-Голея (т. Е. В качестве локальных полиномов наименьших квадратов), существует естественное обобщение неравномерной выборки - она ​​просто намного дороже в вычислительном отношении.

Савицкий-Голей Фильтры в целом

Для стандартного фильтра идея состоит в том, чтобы подогнать полином к локальному набору выборок [используя наименьшие квадраты], а затем заменить центральную выборку значением полинома в центральном индексе (т. Е. В 0). Это означает, что стандартные коэффициенты фильтра SG могут быть сгенерированы путем инвертирования матрицы Вандермонда выборочных признаков. Например, чтобы сгенерировать локальное параболическое соответствие для пяти выборок (с локальными показателями -2, -1,0,1,2), система уравнений расчета A c = y будет выглядеть следующим образом:Y0...Y4Aсзнак равноY

[-20-21-22-10-11-12000102101112202122][с0с1с2]знак равно[Y0Y1Y2Y3Y4],

Выше, - неизвестные коэффициенты полинома наименьших квадратов c 0 + c 1 x + c 2 x 2 . Поскольку значение полинома при x = 0 равно c 0 , вычисление псевдообратной матрицы проектирования (т. Е. C = ( A T A ) - 1 A T y ) даст коэффициенты фильтра SG в верхней строке. В этом случае они будутс0...с2с0+с1Икс+с2Икс2Иксзнак равно0с0сзнак равно(ATA)-1ATY 

[с0с1с2]знак равно[-3121712-3-7-40475-3-5-35][Y0Y1Y2Y3Y4],

с0+с1Икс+с2Икс2с1+2с2Иксс1

Неравномерный отбор проб

ИксNTN0

T-2знак равноИкс-2-Икс0T-1знак равноИкс-1-Икс0T0знак равноИкс0-Икс0T1знак равноИкс1-Икс0T2знак равноИкс2-Икс0

тогда каждая матрица дизайна будет иметь следующую форму:

Aзнак равно[T-20T-21T-22T-10T-11T-12T00T01T02T10T11T12T20T21T22]знак равно[1T-2T-221T-1T-121001T1T121T2T22],

A с0

специалист по обработке данных
источник
звучит так, как будто он движется от O (log (n)) к O (n ^ 2).
EngrStudent - Восстановить Монику
Вот реализация Scala, описанная датейстом вверх.
Среднее ядро
1
@Mediumcore Вы не добавили ссылку на свой оригинальный пост. Кроме того, я удалил его, потому что он не дал ответа на вопрос. Пожалуйста, попробуйте отредактировать запись в datageist, чтобы добавить ссылку; это будет модерироваться после проверки.
Питер К.
4


еNN/2
-

(вывод кого-нибудь?)


±N/2TT
Tяя

Денис
источник
1

Я обнаружил, что есть два способа использовать алгоритм Савицкого-Голея в Matlab. Один раз как фильтр, а другой как функция сглаживания, но в основном они должны делать то же самое.

  1. yy = sgolayfilt (y, k, f): Здесь предполагается, что значения y = y (x) находятся на одинаковом расстоянии по x.
  2. yy = smooth (x, y, span, 'sgolay', Степень): здесь вы можете использовать x в качестве дополнительного ввода, и, ссылаясь на справку Matlab, x не обязательно должен находиться на одинаковом расстоянии!
Jochen
источник
0

Если это поможет, я сделал реализацию C метода, описанного datageist. Бесплатное использование на свой страх и риск.

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

сглаживающий

techwinder
источник