Разложение Фурье, нахождение гармоник, прогноз

[B@R5uk]

Schrodinger's cat в сообщении #1565572 писал(а):

В чем проявляется неустойчивость?

Ну вот, например, для трёх зашумлённых синусоид нужно как минимум от 40 до 60 компонент, чтобы эти самые синусоиды выловить в зависимости от величины шума. Причём, между величиной шума и количеством компонент нет ровно никакой корреляции. С реальным сигналом вообще затык. Может быть я, конечно, что-то не так делаю...

Schrodinger's cat в сообщении #1565572 писал(а):

И как вы ее решаете если не секрет?

Да никак. Я когда-то давно, в студенческие годы по наивности своей мечтал сделать программку, которая будет как раз резать сигнал на элементарные гармоники. Потом с сигналом в таком представлении (по моему видению) можно что угодно делать: сжимать, распознавать ноты, менять длительность, не меняя тона, или наоборот, тон, не меняя длительности. Что хочешь, короче. В результате долгих, мучительных и безуспешных экспериментов я таки допёр, что фаза у гармоник плывёт! Вот это было для меня откровение! (Ни преподы по физике, ни преподы по матану не могли меня на путь истинный наставить! Кучу времени на это угробил). После чего быстро пришло понимание того, что каждая из выше перечисленных задач потребует не только удачной модели, по типу приведённой выше (не синусами едиными аналитики живы), но и заточенных под тип сигнала и тип задачи алгоритмов. Более того, даже представление данных в виде такой модели — та ещё плохо поставленная задача. И интерес у меня к этой проблеме как-то поутих. Ну, просто потому что мои умственные способности ограничены, а то, что я хотел решить и запрограммировать, потребует годы работы для целого коллектива исследователей и программистов.

[Schrodinger's cat]

svv в сообщении #1565554 писал(а):

Просто напишите, удалось ли проверить, или были проблемы.

"Нули" первого уравнения. Довольно таки отличаются от нуля. Вокруг 2.5 крутятся, это я бы сказал проблема а не ноль.

(Оффтоп)

"Нули" второго уравнения. Вполне себе нули. Самые "плохие" $10^{-13}$ вполне нормально.

(Оффтоп)

$a_3$ и $\lambda_3$ у нас не определены, определим их например так:
$a_3=5ie^{-0.7}$
$\lambda_3=+i\frac{2\pi}{15}$
Суммируем к уже имеющимся 2м экспонентам, график принимает вид:

(Оффтоп)

И соответственно "нули" третьего уравнения, тоже вполне нормальные:

(Оффтоп)

С тремя экспонентами "нули" первого уравнения испортились до 3.

(Оффтоп)

B@R5uk в сообщении #1565574 писал(а):

Да никак.

Ясно. Ну такое за милисекунды вряд ли можно сделать, ну за секунды думаю да. У вас целевая (фитнес) функция не такая тяжелая 2 переменных в каждом отсчете, отсчетов у вас сколько? Что то вроде 1024 и того 2048 переменных. Достаточно немного. В нейронных сетях на порядки больше и отжиг справляется нормально и довольно быстро. Именно отжиг, потому что всевозможная муть типа генетического алгоритма, алгоритм роя, муравьиный, пчелиный и т.д. за счет лишней мудрености справится за большее время.
Можно даже стартуя от классического Фурье просто рандомом с мелким разбросом "дрейфовать" к лучшему новому решению.
...
Или я неправильно посчитал. 2048 это у Вас только для одной частоты... тогда переменных порядка миллиона, уже средненько, если Вы эти амплитуды задаете в виде таблицы, а если задать их при помощи некой функции, то в функции будет гораздо меньше переменных и тогда мы опять возвращаемся примерно к 2000 переменных.

Размышления в слух, я бы пробовал так делать.

[svv]

Прекрасно, что Вы это подтвердили в численном эксперименте.
Я выше писал, что
$x_{n+1}-e^\lambda x_n=0$ выполняется, если сигнал состоит из одной экспоненты,
$x_{n+2}-(e^{\lambda_1}+e^{\lambda_2}) x_{n+1}+e^{\lambda_1}e^{\lambda_2} x_n=0$ — если из двух,
$x_{n+3}-(e^{\lambda_1}+e^{\lambda_2}+e^{\lambda_3}) x_{n+2}+(e^{\lambda_1}e^{\lambda_2}+e^{\lambda_2}e^{\lambda_3}+e^{\lambda_3}e^{\lambda_1}) x_{n+1}- e^{\lambda_1}e^{\lambda_2}e^{\lambda_3} x_n=0$ — если из трёх.
Если Вы вместо нужного числа экспонент взяли меньшее — это не страшно: всегда можно считать, что недостающие экспоненты есть, просто у них нулевая "амплитуда". Это законный случай, и соотношения будут выполняться.
А вот если, наоборот, соотношение рассчитано на такое-то число экспонент, а Вы взяли большее, оно выполняться не будет.

Но вообще я имел в виду чисто теоретическую выкладку. Я покажу, что подразумевалось, для первого случая, а Вы попробуйте аналогично проделать для второго и третьего.

svv в сообщении #1565513 писал(а):

1) Пусть $x_n=ae^{\lambda n}$.
Проверьте, что тогда для любого $n$ выполняется
$x_{n+1}-e^\lambda x_n=0$

$x_{n+1}-e^\lambda x_n=ae^{\lambda (n+1)}-e^\lambda\,ae^{\lambda n}=ae^{\lambda (n+1)}-ae^{\lambda (n+1)}=0$
Сделайте для себя дома на листочке, здесь писать не обязательно, мне важно — получилось или нет. Следующим шагом будет пояснение, откуда именно такие соотношения.

[Schrodinger's cat]

svv в сообщении #1565591 писал(а):

Сделайте для себя дома на листочке

Да, я понял почему такие выражения. Все слагаемые взаимно аннигилируют.
Второй сделал на листочке руками, а третий честно сделал в maple

[svv]

Отлично. Представим это по-другому.
У нас $x_n$ означает $n$-й элемент последовательности отсчётов. Символ $x$ без индекса обозначает всю бесконечную (в обе стороны) последовательность $...,x_{-2},x_{-1},x_0, x_1,x_2,...$.
Определим несколько действий над последовательностями. Пусть $x,y,z$ — последовательности, $c$ — число (вещественное или комплексное).

Поэлементное сложение и вычитание.
Запись $z=x+y$ означает $z_n=x_n+y_n$ для любого $n$.
Запись $z=x-y$ означает $z_n=x_n-y_n$ для любого $n$.

Поэлементное умножение на число.
Запись $y=cx$ означает $y_n=cx_n$ для любого $n$.

Сдвиг влево.
Запись $y=\mathsf L x$ означает $y_n=x_{n+1}$ для любого $n$.

Пример. Запись $\mathsf L x-c x$ означает: из $x$, сдвинутой влево, вычтем $x$, умноженную на $c$. Это же можно записать так: $(\mathsf L-c)x$. Подробно:
$\begin{tabular}{|c|c|c|c|c|c|c|c|}\hline &индекс&$-2$&$-1$&$0$&$1$&$2$&... \\ послед-ть &&&&&&& \\ \hline $x$&...&$x_{-2}$&$x_{-1}$&$x_0$&$x_1$&$x_2$&...\\ \hline $\mathsf Lx$&...&$x_{-1}$&$x_{0}$&$x_1$&$x_2$&$x_3$&...\\ \hline $ cx$&...&$cx_{-2}$&$cx_{-1}$&$cx_0$&$cx_1$&$cx_2$&...\\ \hline $(\mathsf L-c)x$&...&$x_{-1}-cx_{-2}$&$x_0-cx_{-1}$&$x_1-cx_0$&$x_2-cx_1$&$x_3-cx_2$&...\\ \hline\end{tabular}$

Теперь вот это

svv в сообщении #1565513 писал(а):

1) Пусть $x_n=ae^{\lambda n}$.
Проверьте, что тогда для любого $n$ выполняется
$x_{n+1}-e^\lambda x_n=0$

можно сформулировать так:
Оператор $\mathsf L-e^\lambda$ убивает экспоненциальную последовательность $x$ вида $x_n=ae^{\lambda n}$ (именно с тем же $\lambda$; с другим показателем не убивает). То есть превращает $x$ в последовательность из одних нулей.

Всё ли понятно?

[Schrodinger's cat]

svv спасибо, пока понятно.
Да, если мы теперь оператор поставим под мультипликатор то при раскрытии скобок будем получать как раз такие выражения. Верно?

$$\prod_{k=1}^{p}(L-e^{\lambda _{k}})$$

Попробую например для второго получить выражение:
$(L-e^{\lambda _{1}})(L-e^{\lambda _{2}})=LL-L(e^{\lambda _{1}}+e^{\lambda _{2}})+e^{\lambda _{1}}e^{\lambda _{2}}$
Да, все понял, спасибо.
Дальше Марпл предлагает это представить в виде степенной последовательности, собственно то, что Вы называете полиномом, корни которого надо найти методом Дженкинса-Трауба.
Но как из разностного уравнения получили степенной... пока недопираю.

[ilghiz]

B@R5uk в сообщении #1565569 писал(а):

Жутко неустойчивая штука. Старый добрый Фурье куда надёжнее, особенно для реальных данных.

так а что же вы хотите, она же не $l_2$ норму минимизирует, а что-то похожее, но, все-таки другое. Да и отсчетов надо брать существенно больше, чем компонент в сигнале. У этого метода есть положительные стороны - если частота может медленно изменяться во времени, то Фурье дает очень размазанный пик, а этот метод эту частоту хорошо улавливает. Но надо на многие нюансы обращать внимание, или таки приводить исходную задачу к честным наименьшим квадратам и минимизировать как есть, но тут только через трехмерные формулировки и выше, правда численно это довольно дорогое удовольствие.

-- 29.09.2022, 10:23 --

B@R5uk в сообщении #1565574 писал(а):

Schrodinger's cat в сообщении #1565572 писал(а):

В чем проявляется неустойчивость?

Ну вот, например, для трёх зашумлённых синусоид нужно как минимум от 40 до 60 компонент, чтобы эти самые синусоиды выловить в зависимости от величины шума. Причём, между величиной шума и количеством компонент нет ровно никакой корреляции. С реальным сигналом вообще затык. Может быть я, конечно, что-то не так делаю...

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

[Schrodinger's cat]

ilghiz в сообщении #1565638 писал(а):

Да и отсчетов надо брать существенно больше, чем компонент в сигнале.

Марпл пишет, что для полного восстановления сигнала нужно $p$ экспонент и соответственно $2p$ комплексных отсчетов.
Поясните пожалуйста $2p$ это достаточно для "существенно больше" или нужно еще больше, что кажется избыточным?
Для Фурье так же нужно в 2 раза больше отсчетов чем детектируется гармоник.

Вопрос немного не по теме.
Т.е. по сути алгоритм Дженкинса-Трауба находит собственные значения? Это вытекает из того что задачу можно решить SVD разложением или я так понимаю любым другим методом на собственные значения. Но тогда почему я ни разу не видел этот алгоритм в списке алгоритмов позволяющих находить собственные значения?

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

Schrodinger's cat в сообщении #1564911 писал(а):

ФЧХ симметрична относительно середины.

Антисимметрична

-- 29 сен 2022, 15:07 --

Schrodinger's cat в сообщении #1565644 писал(а):

Т.е. по сути алгоритм Дженкинса-Трауба находит собственные значения? Это вытекает из того что задачу можно решить SVD разложением или я так понимаю любым другим методом на собственные значения. Но тогда почему я ни разу не видел этот алгоритм в списке алгоритмов позволяющих находить собственные значения?

Алгоритм Дженкинса-Трауба служит для нахождения корней полиномов. Если мы начинаем с построения характеристического полинома, то затем его можно решать и этим алгоритмом. Просто никто по крайней мере для матриц, составленных из чисел, собственные значения через характеристические полиномы не выражает. Есть более эффективные вычислительные алгоритмы.

[Schrodinger's cat]

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

Антисимметрична

Спасибо, буду знать.

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

Просто никто по крайней мере для матриц, составленных из чисел, собственные значения через характеристические полиномы не выражает.

Выражают.
Кроме непосредственного развертывания определителя в характеристический полином, есть метод Данилевского, метод Крылова, метод Леверрье-Фадеева и прочие.

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

Есть более эффективные вычислительные алгоритмы.

Не спорю.

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

Schrodinger's cat в сообщении #1564911 писал(а):

Хочу при помощи разложение Фурье найти все гармоники сигнала а так же потом зная их построить прогноз.

Если мы найдём все гармоники для заданного отрезка сигнала, а потом попробуем на их основе получить прогноз - не то, чтобы ничего не выйдет, но выйдет не то. Мы представляем сигнал суммой функций, заданных на бесконечном отрезке (синусов и косинусов), но используем данные на конечном, нам известном. Можно это понимать так, что строится бесконечно длинный отрезок, составленный из повторений заданного (поезд из одинаковых вагончиков). И если построено точное выражение сигнала на заданном отрезке через синусоиды, то за пределами отрезка те же синусоиды, точно так же точно представляющие сигнал других "вагончиков". Иначе говоря, прогноз совпадёт со значениями заданного нам сигнала. В точности.
Смысл в таком разложении появится, если мы откажемся от точного представления. Может быть, какие-то слагаемые с малой амплитудой будут отброшены, поскольку амплитуда их того же порядка, что шум. Или будут введены веса для компонентов, зависящие от амплитуды. Или некоторые компоненты будут отброшены на основе содержательного анализа, скажем, данная частота связана с работой некоего механизма, который далее работать не будет. Но механически считать Фурье и пытаться построить прогноз - не работает.
Да, и ещё. Алгоритмы ДПФ (БПФ) работают с конечным набором частот, некоей "гребёнкой". А в реальных процессах попадание частоты колебательного процесса в точности на "зубец гребёнки" крайне маловероятно.

-- 29 сен 2022, 15:38 --

Schrodinger's cat в сообщении #1565668 писал(а):

Кроме непосредственного развертывания определителя в характеристический полином, есть метод Данилевского, метод Крылова, метод Леверрье-Фадеева и прочие.

Не то, чтобы их нет. Но они имеют скорее историческое значение. Данилевский предложил метод своего имени в 30е годы ХХ века, Крылов в 1931, в каком Леверье - не вем, но умер он в 1877, а замечательная книга Фадеевых вышла в 1963 (вторым изданием).

-- 29 сен 2022, 15:47 --

Schrodinger's cat в сообщении #1565630 писал(а):

Но как из разностного уравнения получили степенной... пока недопираю.

А тут очень похожая с линейными дифуравнениями теория. Такой же вспомогательный полином, через корни которого выражается решение. Только показатели степеней соответствуют не порядкам производной, а сдвигам. А в решении вместо членов вида $A_ie^{\lambda_i t}$ появляются $A_i\lambda_i^t$, где лямбды - корни вспомогательного полинома (некратные, если есть кратные - тоже похоже получается)

[Schrodinger's cat]

Евгений Машеров где Вы раньше были? Такой содержательный ответ бы в начале темы... Но нет я бы все равно сам лично в этом убедился...
Вот ровно то, что вы написали, я и наблюдал в своих экспериментах, чуть дальше если вы почитаете и собственно написал что с Фурье я закончил.
В методе Прони тоже пишут, что он ищет аппроксимацию и это несколько напрягает, может быть такая же бессмысленность как и у Фурье, но интересное свойство точного выделения частот заставляет меня все же подробно изучить метод.
А вот способ через разложение SVD (спасибо ilghiz) уже напротив работает с сингулярными числами, векторами и выделяет закономерности из ряда (знаю по методу АГК, но к Прони еще не применял) и это так же интересно изучить. Пусть даже только линейные закономерности, но это уже лучше чем Фурье.
Это насколько я смог понять.

Кстати кроме нейронных сетей есть методы для нелинейного анализа рядов?

[ilghiz]

Schrodinger's cat в сообщении #1565644 писал(а):

ilghiz в сообщении #1565638 писал(а):

Да и отсчетов надо брать существенно больше, чем компонент в сигнале.

Марпл пишет, что для полного восстановления сигнала нужно $p$ экспонент и соответственно $2p$ комплексных отсчетов.

это утверждение верно только в отсутствие шума в сигнале и неустойчивости алгоритма из-за ошибок округления в машинной арифметике.

То, где я это использую:

пример номер 1(МРТ): данных отсчетов около 3000-10000, компонент (число экспонент) около 20, сигнал/шум существенно меньше 1 для второго десятка компонент, Фурье достоверно видит первые 10, а этот метод немного с бубнами (поиском оптимального числа компонент) находит все, что надо.

пример номер 2(ЯМР): данных много, около миллиона отсчетов, частота у всех плывет на сотни PPM за время сбора данных, Фурье видит несколько самых ярких пиков, но они довольно размыты, остальные не видны вообще, а этот метод видит где-то до 50 пиков и устойчив к слабому изменению частоты.

Классический пример: 100 компонент, 1000 отсчетов, шум в виде еще около 100 или даже больше компонент с меньшей интенсивностью, чем основные. В этом случае, этот метод со свистом продувает Фурье.

Schrodinger's cat в сообщении #1565644 писал(а):

Т.е. по сути алгоритм Дженкинса-Трауба находит собственные значения? Это вытекает из того что задачу можно решить SVD разложением или я так понимаю любым другим методом на собственные значения. Но тогда почему я ни разу не видел этот алгоритм в списке алгоритмов позволяющих находить собственные значения?

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

[Schrodinger's cat]

ilghiz познавательно, понятно, спасибо!

[svv]

Schrodinger's cat в сообщении #1565630 писал(а):

Да, если мы теперь оператор поставим под мультипликатор то при раскрытии скобок будем получать как раз такие выражения. Верно?

$$\prod_{k=1}^{p}(L-e^{\lambda _{k}})$$

Попробую например для второго получить выражение:
$(L-e^{\lambda _{1}})(L-e^{\lambda _{2}})=LL-L(e^{\lambda _{1}}+e^{\lambda _{2}})+e^{\lambda _{1}}e^{\lambda _{2}}$

Всё верно Кстати, "произведение" двух или нескольких сдвигов можно обозначить степенью, например, $\mathsf{LLL}=\mathsf L^3$. Ясно, что это сдвиг на 3 шага влево.
Обозначим ещё для краткости $r_k=e^{\lambda_k}$. Тогда оператор

Schrodinger's cat в сообщении #1565630 писал(а):

$(L-e^{\lambda _{1}})(L-e^{\lambda _{2}})=LL-L(e^{\lambda _{1}}+e^{\lambda _{2}})+e^{\lambda _{1}}e^{\lambda _{2}}$

можно записать так:
$(\mathsf L-r_1)(\mathsf L-r_2)=\mathsf L^2-(r_1+r_2)\mathsf L+r_1r_2$
Из левой части видно, что этот оператор убивает (или "аннулирует", если так больше нравится) последовательность $x_n=a_1e^{\lambda_1 n}+a_2e^{\lambda_2 n}$.
При этом оператор $\mathsf L-r_1$ убивает $a_1e^{\lambda_1 n}$ и умножает на некоторое число $a_2e^{\lambda_2 n}$.
Ну, а оператор $\mathsf L-r_2$ убивает $a_2e^{\lambda_2 n}$ и умножает на некоторое число $a_1e^{\lambda_1 n}$. А будучи применены последовательно (в любом порядке), эти операторы аннулируют сумму.

Очень легко решить уравнение $\bigl(r-4\bigr)\bigl(r-5\bigr)\bigl(r-(2+3i)\bigr)\bigl(r-(2-3i)\bigr)=0$ в комплексных числах. Очевидно, его корни
$r_1=4\quad r_2=5\quad r_3=2+3i\quad r_4=2-3i$
Если же раскрыть скобки и привести подобные, получим второе уравнение
$r^4-13r^3+69r^2-197r+260=0$,
которое решить намного труднее, чем первое (а также намного труднее, чем раскрыть скобки в первом уравнении). Это тем более справедливо для реалистичных случаев большой степени и нецелых коэффициентов.

Совершенно то же у нас. Сначала мы анализируем длинную последовательность отсчётов и находим, что она удовлетворяет разностному уравнению
$x_{n+4}-13x_{n+3}+69x_{n+2}-197x_{n+1}+260x_n=0$
Или, в наших обозначениях, последовательность $x$ аннулируется оператором
$\mathsf L^4-13\mathsf L^3+69\mathsf L^2-197\mathsf L+260$
Но как из этого извлечь гармоники? Оператор надо разложить на множители вида $\mathsf L-r_k$. Для этого берём (характеристическое) алгебраическое уравнение $r^4-13r^3+69r^2-197r+260=0$, решаем (или, что то же, раскладываем на множители левую часть). Получаем корни
$r_1=4\quad r_2=5\quad r_3=2+3i\quad r_4=2-3i$.
Значит, и оператор допускает разложение
$(\mathsf L-r_1)(\mathsf L-r_2)(\mathsf L-r_3)(\mathsf L-r_4)$
с теми же $r_k$. После этого $\lambda_k=\ln r_k$.