Russian Qt Forum
Апрель 29, 2024, 07:40 *
Добро пожаловать, Гость. Пожалуйста, войдите или зарегистрируйтесь.
Вам не пришло письмо с кодом активации?

Войти
 
  Начало   Форум  WIKI (Вики)FAQ Помощь Поиск Войти Регистрация  

Страниц: [1] 2   Вниз
  Печать  
Автор Тема: Восстановление синусоиды  (Прочитано 7003 раз)
Racheengel
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 2679


Я работал с дискетам 5.25 :(


Просмотр профиля
« : Март 18, 2021, 11:04 »

Всем Привет,

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

What is the 11 in the C++11? It’s the number of feet they glued to C++ trying to obtain a better octopus.

COVID не волк, в лес не уйдёт
Igors
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 11445


Просмотр профиля
« Ответ #1 : Март 18, 2021, 12:07 »

Делал такое для кардиологии в начале 90-х, и помню очень немного. Есть простая и удобная штука: "автокорреляционная ф-ция". Два вектора (одинакового размера, положительные) просто множатся один на другой, макс сумма произведений получается при макс совпадении. Ну пытаемся взять меньше/больше точек на период и сравниваем произведения деленные на число точек. Если частота плавает - там хужее
« Последнее редактирование: Март 18, 2021, 12:46 от Igors » Записан
m_ax
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 2094



Просмотр профиля
« Ответ #2 : Март 18, 2021, 12:27 »

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

Над водой луна двурога. Сяду выпью за Ван Гога. Хорошо, что кот не пьет, Он и так меня поймет..

Arch Linux Plasma 5
Racheengel
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 2679


Я работал с дискетам 5.25 :(


Просмотр профиля
« Ответ #3 : Март 18, 2021, 12:29 »

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

What is the 11 in the C++11? It’s the number of feet they glued to C++ trying to obtain a better octopus.

COVID не волк, в лес не уйдёт
Racheengel
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 2679


Я работал с дискетам 5.25 :(


Просмотр профиля
« Ответ #4 : Март 18, 2021, 12:47 »

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

What is the 11 in the C++11? It’s the number of feet they glued to C++ trying to obtain a better octopus.

COVID не волк, в лес не уйдёт
m_ax
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 2094



Просмотр профиля
« Ответ #5 : Март 18, 2021, 13:03 »

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

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

Над водой луна двурога. Сяду выпью за Ван Гога. Хорошо, что кот не пьет, Он и так меня поймет..

Arch Linux Plasma 5
Igors
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 11445


Просмотр профиля
« Ответ #6 : Март 18, 2021, 13:07 »

Да, спасибо, это уже нечто близкое к моей задаче.
Но проблема в том, что у меня "неоткуда" брать второй вектор - входной сигнал только один ("битый"), а оригинальные его период и амплитуда мне не известны.
Известно только, что это некий синус.
Вот есть 500 самплов. Берете первые 100, считаете суммы (0..99) * (100..199),
(0..99) * (200..399) и.т.д  Осреднили все суммы, получили какое-то значение

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

А можете посоветовать "хорошую" С или С++ либу для фиттинга синусоиды?
А то "в интернетах" всё завалено питоном Грустный
Да-да, "программист - это тот кто пользуется готовым" (тут один чувак говорил, кажется из Ярославля)

Да там проще самому написать) Завтра постараюсь опубликовать) Мне тож это нужно будет позже)
А для Вас есть соседняя тема  Улыбающийся
« Последнее редактирование: Март 18, 2021, 13:10 от Igors » Записан
ssoft
Программист
*****
Offline Offline

Сообщений: 579


Просмотр профиля
« Ответ #7 : Март 18, 2021, 21:52 »

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


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

Сообщений: 2679


Я работал с дискетам 5.25 :(


Просмотр профиля
« Ответ #8 : Март 19, 2021, 12:39 »

Есть же обычное преобразование Фурье, которое выделит и период и фазовый сдвиг синусоиды. Тем более, что начало гарантированно совпадает с концом по амплитуде (периоду). А если количество сэмплов будет кратно степени 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 поможет выделить синусоиду из сигнала...
Записан

What is the 11 in the C++11? It’s the number of feet they glued to C++ trying to obtain a better octopus.

COVID не волк, в лес не уйдёт
m_ax
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 2094



Просмотр профиля
« Ответ #9 : Март 19, 2021, 13:03 »

Цитировать
Хорошо, но как его к данному конкретному случаю применить?
Да, всё хорошо будет) Сейчас допишу) Делов то там..
Записан

Над водой луна двурога. Сяду выпью за Ван Гога. Хорошо, что кот не пьет, Он и так меня поймет..

Arch Linux Plasma 5
Igors
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 11445


Просмотр профиля
« Ответ #10 : Март 19, 2021, 13:10 »

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

И такие задачи часто называются "curve fitting" или "fit curve", добавьте "sine" или "sinusoidal", там должны были рыться многие
Записан
Racheengel
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 2679


Я работал с дискетам 5.25 :(


Просмотр профиля
« Ответ #11 : Март 19, 2021, 13:20 »

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

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

What is the 11 in the C++11? It’s the number of feet they glued to C++ trying to obtain a better octopus.

COVID не волк, в лес не уйдёт
m_ax
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 2094



Просмотр профиля
« Ответ #12 : Март 20, 2021, 14:17 »

Работает)
Держите:
Код
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;
}
 
« Последнее редактирование: Март 22, 2021, 13:20 от m_ax » Записан

Над водой луна двурога. Сяду выпью за Ван Гога. Хорошо, что кот не пьет, Он и так меня поймет..

Arch Linux Plasma 5
m_ax
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 2094



Просмотр профиля
« Ответ #13 : Март 20, 2021, 17:50 »

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

Над водой луна двурога. Сяду выпью за Ван Гога. Хорошо, что кот не пьет, Он и так меня поймет..

Arch Linux Plasma 5
Racheengel
Джедай : наставник для всех
*******
Offline Offline

Сообщений: 2679


Я работал с дискетам 5.25 :(


Просмотр профиля
« Ответ #14 : Март 22, 2021, 15:31 »

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

What is the 11 in the C++11? It’s the number of feet they glued to C++ trying to obtain a better octopus.

COVID не волк, в лес не уйдёт
Страниц: [1] 2   Вверх
  Печать  
 
Перейти в:  


Страница сгенерирована за 0.084 секунд. Запросов: 22.