Russian Qt Forum

Всем Привет,

появилась такая задача: на вход поступает сигнал (цифровой) синусоидальной формы, но достаточно зашумленный (случайные пики, пропуски данных и т.д.)
Известно, что исходный сигнал это синусоида с неким постоянным периодом.
Объем данных - массив из примерно 500-1000 сэмплов. В пределах этого массива сигнал "закольцован", т.е. его начало гарантированно совпадает с концом по  значениям амплитуды.
Вот требуется восстановить исходную синусоиду (период в общем случае не известен).
Есть ли какие-либо известные алгоритмы для этого?
"гугл не помог" пока :(

Делал такое для кардиологии в начале 90-х, и помню очень немного. Есть простая и удобная штука: "автокорреляционная ф-ция". Два вектора (одинакового размера, положительные) просто множатся один на другой, макс сумма произведений получается при макс совпадении. Ну пытаемся взять меньше/больше точек на период и сравниваем произведения деленные на число точек. Если частота плавает - там хужее

Это фитинг называется: минимизируем дисперсию.
Проходная задача для экспериментаторов) 

Да, спасибо, это уже нечто близкое к моей задаче.
Но проблема в том, что у меня "неоткуда" брать второй вектор - входной сигнал только один ("битый"), а оригинальные его период и амплитуда мне не известны.
Известно только, что это некий синус.

А можете посоветовать "хорошую" С или С++ либу для фиттинга синусоиды?
А то "в интернетах" всё завалено питоном :(

Да там проще самому написать) Завтра постараюсь опубликовать) Мне тож это нужно будет позже)

Вот есть 500 самплов. Берете первые 100, считаете суммы (0..99) * (100..199),
(0..99) * (200..399) и.т.д  Осреднили все суммы, получили какое-то значение

Потом все то же для 101 точки, для 102 и.т.п. Выбираете тот вариант где сумма больше. Если есть фаза (первая синусоида неполная), то пробуете начать с точки 0, потом с 1 и.т.д.

Да-да, "программист - это тот кто пользуется готовым" (тут один чувак говорил, кажется из Ярославля)

А для Вас есть соседняя тема  :)

Есть же обычное преобразование Фурье, которое выделит и период и фазовый сдвиг синусоиды. Тем более, что начало гарантированно совпадает с концом по амплитуде (периоду). А если количество сэмплов будет кратно степени 2, то можно и fft   алгоритм  применить.

Хорошо, но как его к данному конкретному случаю применить? Т.е. есть массив данных (сигнал) и его длина, типа std::vector<char>. И всё.
Я смотрел, к примеру, это: https://docs.opencv.org/master/de/dbc/tutorial_py_fourier_transform.html
Но как-то мало релевантно...

То есть, в идеале что надо бы:

bool fit_sine(const std::vector<char>& input_signal, std::vector<char> &output_sine_wave);

где output_sine_wave - "чистая" синусоида такой же частоты и амплитуды, как и input_signal.

Я не совсем понимаю, как FFT поможет выделить синусоиду из сигнала...

Да, всё хорошо будет) Сейчас допишу) Делов то там..

Да, похоже здесь лучше искать "готовое проверенное". Может не стоит завязываться на Фурье, там смысл - представить сигнал суммой синусоид, у каждой своя частота, фаза и амплитуда, Вам нужна только первая. И да, на входе нужен только набор точек.

И такие задачи часто называются "curve fitting" или "fit curve", добавьте "sine" или "sinusoidal", там должны были рыться многие

http://mariotapilouw.blogspot.com/2012/03/sine-fitting.html

похоже на правду... как затестю, сообщу...

Работает)
Держите:

C++ (Qt)
#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <random>
 
#include <specmath/minsearch.h>
#include <specmath/grid.h>
 
constexpr double pi = M_PI;
constexpr size_t Sample_Size = 1000;
 
int main()
{
    specmath::grid1d<double> grid(0.0, 2*pi, Sample_Size);
 
    std::random_device rd;
    std::mt19937 gen(rd());
    const double dispersion = 0.2;
    std::normal_distribution<double> dist(0.0, dispersion);
 
    const double A = 1.0; // Amplitude
    const double omega = 1.0; // frequency
 
    grid.apply([&](double x)
    {
        // generate random data
       return A * sin(omega * x) + dist(gen);
    });
 
 
    specmath::minsearch<double> minsearch(1.0, 10);
 
    std::vector<double> x = {0.0, 0.0}; // init point
 
    double sqr_sigma = minsearch.find_minimum([&](const std::vector<double> & v)
    {
        auto it = grid.first_point(); // x = 0.0, .. 2*pi
        double res = 0.0;
 
        double m_A  = v[0];
        double m_omega = v[1];
 
        for (const auto & val : grid)
        {
            double xi = *(it++);
            double yi = val;
 
            double d = m_A * sin(m_omega * xi) - yi;
 
            res += d*d;
        }
 
        return res;
 
    }, x, specmath::breaker<double>(1e-6));
 
    std::cout << "Amplitude = " << x[0] <<  " omega = " << x[1] << std::endl;
 
    std::cout << "sqr_sigma = " << sqr_sigma << std::endl;
 
    {
        std::ofstream out("data.txt");
        grid.serialize(out);
    }
 
 
    return 0;
}
 

Да и ещё.. minsearch - многопоточный, так что по хорошему нужно лочить вызываемую функцию...
Конкретно в этом примере проблем быть не должно, но имейте в виду..

спасибо)
а specmath это что такое?

Это такая штука, которая минимум функции многих переменных находит)
В данном случае размерность = 2: Amplitude и omega)

Я имею в виду, там просто в исходниках нигде ни лицензии не указано, ни авторства... поэтому и интересуюсь :)

А это мои, так сказать, наработки) 

Ну, выглядит довольно серьёзно) а планы на будущее в развитии есть? на гитхаб-гитлаб там положить, или как либу выпускать?

Да что-то пока руки не доходят.. Это больше по работе.. Такое расширение boost'а.. Пишу понемногу для своих нужд  :)

Вам шашечки или ехать?  :)

Мне то ехать, а кастомерам - еще и с шашечками...

Да и просто интересно, откуда код был (думал - GSL-like), а оказалось - авторский)