Получить характеристики сигнала с помощью МНК

[Kevsh]

У меня есть набор точек, причём я знаю общий вид функции: $s(t)=B+A\sin{(2 \pi ft + \varphi)}$. Проблема в том, что я не знаю ни $A$, ни $B$, ни $f$, ни $\varphi$. Я хочу их определить с помощью МНК, но не понимаю, что мне делать с $\varphi$ и $f$ (эти коэффициенты постоянны и не зависят от $t$). Как-то вынести эти два параметра и привести это всё к стандартному для МНК виду ($f(x)=a_1\varphi_1(x)+a_2\varphi_2(x)...$) нельзя, как мне кажется. Поиск в интернете мне выдал только что-то типа ряда Фурье, но мне нужна только одна синусоида и постоянная составляющая (я точно знаю, что точки по этому закону отмечены). Как тут можно поступить?

[alisa-lebovski]

В явном виде все равно не получится. Каким-нибудь методом последовательных приближений минимизируете сумму квадратов отклонений.

[ipgmvq]

Kevsh в сообщении #1585869 писал(а):

Как тут можно поступить?

Можно найти $B$ как выборочную среднюю $s(t)$, $A$ через амплитуду (через максимальное значение выборки $s(t)$ минус минимальное), $f$ как частоту с максимумом на спектре, а потом, если нужно, искать $\varphi$. Или Вы что-то более точное, несмещенное, с гипотезой об распределении ошибок хотели?

[Kevsh]

ipgmvq
Дело в том, что эта функция– результат измерения в сумме с синусоидальной погрешностью. То есть по сути меня вообще интересует как можно более точно вычислить коэффициент $B$. Выборочная средняя даёт достаточно точный результат, но я подумал, раз я знаю вид зависимости, то можно как-то найти саму функцию и по ней уже определить $B$.

[Евгений Машеров]

Общий подход - нелинейная регрессия. Метод Левенберга-Марквардта, скажем.
Удобно перепараметризовать, избавившись от одного нелинейного параметра, как $s(t)=B+A_1\sin(2\pi f t)+A_2\cos(2\pi ft)$
Если ошибка, отягощающая измерения, невелика, а отрезок достаточно длинен, можно перейти ко вторым разностям и построить регрессию сигнала на вторые разности. Свободный член будет B, а как связан коэффициент при регрессоре и частота - предоставляется получить самому (на самом деле мне просто лень выписывать, и я прячусь за правилом "не давать полного решения", но там действительно просто). Если нужна ещё и фаза - регрессию сигнала на синусы и косинусы уже оцененной на прошлом шаге частоты.

[ipgmvq]

Kevsh в сообщении #1585884 писал(а):

результат измерения в сумме с синусоидальной погрешностью

А измерение происходит через обсолютно одинаковые отрезки времени? И выборка многократно превосходит по продолжительности длину цикла?

-- 18.03.2023, 21:17 --

Потому что если на оба вопроса ответы да (и там нет других случайных погрешностей, или есть?), то хотя распредедение ошибки будет иметь плотность вероятности
$f(x) = \frac{A}{\sqrt{1-x^2}}$
но это будет подарок, а не ошибка, потому что она регулярная (случайность обусловлена только фазой). Думаю, можно смело брать среднюю.

[ipgmvq]

Точнее плотность:
$f(s) = \frac{1}{A \pi \sqrt{1-(\frac{s}{A})^2}}$

[Александрович]

Kevsh в сообщении #1585869 писал(а):

У меня есть набор точек, причём я знаю общий вид функции: $s(t)=B+A\sin{(2 \pi ft + \varphi)}$.

Покажите этот набор.

[Евгений Машеров]

Какого рода ошибкой отягощён сигнал?

[Kevsh]

Александрович
[64.5049139938209 67.8223436757140 72.4760976391098 76.9773303280885 79.4198744610976 78.5818892925456 74.8971689391827 70.0909793980171 65.9570193799335 63.6174765994403 63.5354596305803 65.7249719562536 69.7603564745536 74.5720681788433 78.4000811411012 79.4734695741806 77.2380697049508 72.8197571958683 68.1194380402119 64.6751617108084 63.2943193678395 64.2253086891689 67.3053514922505 71.8566206632193 76.4820148895542 79.2805816151615 78.8717749876296 75.4678852479630 70.6949575707295 66.3976558680346 63.7936998809994 63.4166729510651 65.3311905279433 69.1761069601301 73.9752021858035 78.0373688733821 79.5265289248018 77.6788634001248 73.4348858978174 68.6704620839186 65.0078483339424 63.3424176569383 63.9808833921220 66.8130172868168 71.2391825334004 75.9566992674927 79.0873700072592 79.1114498701974 76.0164222631263 71.3075732059835 66.8664415438455 64.0062823981233 63.3352177453742 64.9692406271896 68.6082629370231 73.3668520855100 77.6317868513622 79.5234034929759 78.0798477537931 74.0421632883020 69.2402641151579 65.3733945905512 63.4280325825031 63.7723146322078 66.3472757823477 70.6273619315867 75.4054907923039 78.8420012181945 79.2987489179218 76.5386055456489 71.9254179737826 67.3616147031868 64.2546522137271 63.2912970195925]

-- 26.03.2023, 11:52 --

Евгений Машеров
Ошибка возникает из-за эффекта растекания спектра при измерении фазы сигнала с помощью ДПФ

[Евгений Машеров]

Я о том, надо ли значения $s(t)$ считать точными $s(t)=B+A\sin{(2 \pi ft + \varphi)}$, или же следует рассматривать $\hat{s(t)}=s(t)+ \varepsilon$, и что известно тогда об эпсилонах?

[ilghiz]

Евгений Машеров в сообщении #1585907 писал(а):

Общий подход - нелинейная регрессия. Метод Левенберга-Марквардта, скажем.
Удобно перепараметризовать, избавившись от одного нелинейного параметра, как $s(t)=B+A_1\sin(2\pi f t)+A_2\cos(2\pi ft)$

Плюсую, но с маленьким дополнением, задачу ТС можно и нужно свести к задаче поиска по одному параметру, которую решить или Ньютоном, или методом деления отрезка пополам, а именно исходную задачу надо привести к виду

$$\min_{\phi} \min_{A,B,C} \int || s(t) - A - B \sin(2\pi \phi t) - C \cos(2\pi \phi t)||_2^2 dt$$

причем внутреннюю минимизацию по $A$, $B$, $C$ надо аналитически решить и полученный крокодил подставить в минимизацию по $\phi$. Крокодил, если лениво, можно на мепле/математике/максмие посчитать, но он реально считается и это, ИМХО, единственно адекватно быстрый способ решения этой задачи, которую мне уже много раз приходилось решать.

[Евгений Машеров]

Можно, конечно, и так. Линейная часть - попросту регрессия, сумма квадратов, как ЦФ для минимизации по параметру "частота". Но мне сдаётся, что Левенберг-Марквардт лучше. Там тоже линейная регрессия считается, но меньше вызовов. Собственно, она для уточнения частоты. Может пары вызовов хватить.

[ilghiz]

Евгений Машеров в сообщении #1586892 писал(а):

Но мне сдаётся, что Левенберг-Марквардт лучше.

Не соглашусь. ТС решает аппроксимацию, значит у него есть шум. Какой именно там у него шум - хз, и в этой задаче будет много локальных минимумов. Как их разруливать Левенберг-Марквардт будет - реально можно только повеситься, программируя всякие симулейтед-анниелинги с регионами поверх оного, чтобы он из минимума вылез. Если исходных оцифрованных данных $N$, то мной предложенный метод на одну итерацию будет требовать $kN$ арифметических операций, где $k$ примерно равно $30$, что в общем-то не есть супер много. И Ньютон там одномерный пишется, и можно честно обойти все регионы, оценивая по сходимости Ньютона регион, куда он сходится, а куда - нет, так как фактически область задания можно и сверху и снизу ограничить.

[ilghiz]

Еще дополню, если ТС по своим исходным данным сделает Фурье и посчитает нормы Фурье коэффициентов, то каждый максимум этого спектрального представления можно воткнуть в начальное приближение одномерной решалки, как я предложил, и оно сойдется Ньютоном за пару итераций (ЛМ с таким начальным приближением и допосчитанными линейными коэффициентами тоже хорошо сойдется, но мой метод лучше :))) , так как можно случайно поиметь супер плохо скалированные параметры минимизации и поиметь большие численные ошибки у ЛМ, а не пользуя регулярно методы типа ЛМ это можно и не заметить и долго потом расстраиваться, что не работает.

Тогда у ТС будет набор решений и невязок к ним и ему надо будет решить стоит ли только на невязку каждого решения смотреть, или он, исходя из свойств своей задачи, еще поверх что-то да прикрутит и найдет как отсортировать все его получаемые решения.

Понятно, что самые яркие пики в спектральной области будут представлять решения с минимальной невязкой, но вот уже два пика с одинаковой яркостью могут дать на наименьших квадратах слегка различные невязки в силу того, что число периодов может не супер хорошо лечь на Фурье сетку.