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

[Schrodinger's cat]

Подытожу:
По методу предложенному ilghiz продвижение 0. Я клюнул на то, что часть с разложением SVD мне знакома, но к сожалению на этом все и закончилось. И хоть в книге Марпла это упоминается, сдвинуться с места не удалось.

svv
С комплексными сигналами доселе не работал.
Разобрался как их представлять на графике, координат то теперь две. Надеюсь правильно.

(Оффтоп)

https://postimages.org/ru/

Самое большое мое продвижение заключается в осознании, что в методе 2 СЛАУ и первая моя цель заключается в отыскании решения 11.18 по книге, стр. 372:
$\begin{bmatrix} x[p] & x[p-1] & ... & x[1] \\ x[p+1] & x[p] & ... & x[2] \\ \vdots & \vdots & \ddots & \vdots \\ x[2p-1] & x[2p-2] & ... & x[p] \\ \end{bmatrix} \begin{bmatrix} a[1] \\ a[2] \\ \vdots \\ a[p] \end{bmatrix} = \begin{bmatrix} x[p+1] \\ x[p+2] \\ \vdots \\ x[2p] \end{bmatrix}$

Потом второй целью является СЛАУ 11.12 (стр. 371):
$\begin{bmatrix} z_{1}^{0} & z_{2}^{0} & ... & z_{p}^{0} \\ z_{1}^{1} & z_{2}^{1} & ... & z_{p}^{1} \\ \vdots & \vdots & \ddots & \vdots \\ z_{1}^{p-1} & z_{2}^{p-1} & ... & z_{p}^{p-1} \\ \end{bmatrix} \begin{bmatrix} h_1 \\ h_2 \\ \vdots \\ h_p \end{bmatrix} = \begin{bmatrix} x[1] \\ x[2] \\ \vdots \\ x[p] \end{bmatrix}$

Вы пишите $n+1$ т.е. точка из будущего, которой у нас нет.
В книге же пишут, что у нас всего $2p$ точек, а суммы в индексах относительно $p$ в СЛАУ 11.18 это видно. Таким образом мы не выходим за пределы известного ряда.

Первое СЛАУ составил, решил (по крайней мере думаю что решил, получился вектор неких чисел...).
Куда теперь эти коэффициенты совать?... подумаю некоторое время. MathCAD мне все равно только вечером доступен, так что до следующего вечера.

[Pphantom]

Тема перемещена из форума «Помогите решить / разобраться (М)» в форум «Карантин»
по следующим причинам:

- да, набирайте (краткие инструкции: «Краткий FAQ по тегу [math]» и видеоролик Как записывать формулы), но кажется, что тут проще сначала подумать, нужна ли соответствующая информация в теме.

Исправьте все Ваши ошибки и сообщите об этом в теме Сообщение в карантине исправлено.
Настоятельно рекомендуется ознакомиться с темами Что такое карантин и что нужно делать, чтобы там оказаться и Правила научного форума.

[Pphantom]

Тема перемещена из форума «Карантин» в форум «Помогите решить / разобраться (М)»

[svv]

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

Самое большое мое продвижение

А ещё большим продвижением будет осознание, что мы с ilghiz говорим совершенно об одном и том же. Чуть в разных обозначениях.

[B@R5uk]

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

На этом и основан метод Прони.

Клёвая штука я смотрю! Даже позволяет затухающие по экспоненте сигналы фитить. Большое спасибо за упоминание и объяснение. А то я всё голову ломал: чем бы сумму синусов не кратных частот на МК фитить? Так и чувствовал, что что-то здесь такое возможно, потому что для одинокого синуса отношение амплитуды производной к амплитуде сигнала без проблем позволяет частоту посчитать.

[svv]

Рад, что понравилось.
МК — это микроконтроллер? Меня вот что тогда беспокоит. Чтобы получить коэффициенты линейного разностного уравнения, надо решить СЛАУ. Чтобы по известным коэффициентам восстановить параметры экспонент ($\lambda_k$), надо решить одно алгебраическое уравнение высокой степени (равной числу экспонент). Это микроконтроллеру под силу? (Скорее всего, я просто отстал от жизни )

[ilghiz]

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

МК — это микроконтроллер?

хоть я не B@R5uk, но попытаюсь ответить. Да, МК — это микроконтроллер, у которого обычно не больше 512КБайт оперативной памяти (часто бывает меньше).

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

Меня вот что тогда беспокоит. Чтобы получить коэффициенты линейного разностного уравнения, надо решить СЛАУ. Чтобы по известным коэффициентам восстановить параметры экспонент ($\lambda_k$), надо решить одно алгебраическое уравнение высокой степени (равной числу экспонент). Это микроконтроллеру под силу?

Матрица в СЛАУ после наименьших квадратов довольно маленькая, и ее размер зависит только от числа экспонент. 50 экспонент - это только 10КБ (если использовать эрмитовость и одинарную точность), а сам сигнал можно на лету втыкать в эту матрицу не храня всю исходную постановку.

Задача нахождения корней - решается через задачу на собственные значения матрицы той же размерности.

То есть в 512К реально вмещать задачи на 350 экспонент, но вот о том как быстро и устойчиво они будут сходятся я тихо промолчу.

То есть у МК за 2 бакса, например, esp32, можно ожидать, что где-то раз в секунду можно фиттить те самые 350 экспонент, или решать более приземленные задачи с 20 экспонентами примерно каждую миллисекунду.

Я сам решал под тысячу экспонент на довольно жирной FPGA (у которой было около 10МБ памяти), но на ней был не метод Прони, а похожий, но через тензорное разложение, и, фактически на FPGA строился только тензор, а потом он в довольно компактном виде пересылался на небольшой микрокомпьютер с линуксом, в котором уже все и досчитывалось.

[svv]

ilghiz, спасибо за ответ и информацию.

[ilghiz]

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

Самое большое мое продвижение заключается в осознании, что в методе 2 СЛАУ и первая моя цель заключается в отыскании решения 11.18 по книге, стр. 372:

во, а теперь попользуйте полученные $a[1], ..., a[p]$ в качестве коэффициентов полинома $p$-ой степени (свободный член будет тогда единицей) и найдите его корни и будет вам счастье.

С тензорным разложением, просто строите 3-х мерный тензор, сдвигая по второй и третьей координате ваши данные. Тогда если данных у вас было $N$ и вы планируете найти $R$ компонент, то надо строить тензо размерности $(N-2R) \times R \times R$, и напускаете на него любое доступное тензорное разложение.

Вариантов для тензорного разложения много, можно почитать кратко и подобрать один из них у меня в doi: 10.1002/nla.297 или в вытащить код в нашем патенте US десять семь семь три ноль девяноста три B2, а про то, как получать очень многомерные тензоры, можно почитать опять же у меня в doi: 10.1038/nmeth900

[Schrodinger's cat]

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

2) $x_{n+2}-(e^{\lambda_1}+e^{\lambda_2}) x_{n+1}+e^{\lambda_1}e^{\lambda_2} x_n=0$

Как получилось данное уравнение?
Из 11.17?

$$\sum_{m=0}^{p}a[m]x[n-m]=\sum_{i=0}^{p}h_iz_i^{n-p}\sum_{m=0}^{p}a[m]z_i^{p-m-1}=0$$

Левая часть понятна, из нее образованы левые матрицы СЛАУ.
И то, что
$a[0]=1$,
$a[1]=-(e^{\lambda_1}+e^{\lambda_2})$,
$a[2]=e^{\lambda_1}e^{\lambda_2}$
Понятно.
Непонятно $-(e^{\lambda_1}+e^{\lambda_2})$ и $e^{\lambda_1}e^{\lambda_2}$ как получились эти формулы?
Средняя часть с двумя сумматорами 11.17 нам не нужна? Кроме того непонятно, что она выражает.

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

во, а теперь попользуйте

Вот сюда 11.14

$$\Phi (z)=\sum_{m=0}^{p}a[m]z^{p-m}$$

А что есть левая часть?

[svv]

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

Как получилось данное уравнение?

Я собирался рассказать (не опираясь на формулы из книги), но хотелось бы, чтобы сначала Вы проверили, что выполняются те соотношения между отсчётами. Просто напишите, удалось ли проверить, или были проблемы.

[ilghiz]

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

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

во, а теперь попользуйте

Вот сюда 11.14

$$\Phi (z)=\sum_{m=0}^{p}a[m]z^{p-m}$$

А что есть левая часть?

А зачем вам левая часть? Schrodinger's cat вот я вам написал

Цитата:

во, а теперь попользуйте полученные $a[1], ..., a[p]$ в качестве коэффициентов полинома $p$-ой степени (свободный член будет тогда единицей) и найдите его корни и будет вам счастье.

а вы опять что-то да напридумали, ай-ай-ай!

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

[Schrodinger's cat]

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

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

Сделаю это сегодня вечером, вчера не успел.

Вижу, что для следующего шага понадобится алгоритм Дженкинса-Трауба для отыскания корней полинома степени n.
В качестве альтернативы вроде есть методы: Лагерра, 2d приближений и через собственные числа видимо что пытался донести уважаемый ilghiz.
Но я надеюсь, что для начала эту задачу осилит метод polyroots, mathcad'а.

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

И еще, я конечно все понимаю...

Я конечно премного благодарен Вам за попытки помочь, но перечитав 20 раз вдумчиво я так и не понимаю что вообще делать.
Пока что больше понятно из того что написано в книге Марпла, спасибо svv, и собственно из объяснений svv.
Из Вашего объяснения мне даже тут непонятно в какой полином это все ставить?! В канонический? А может 11.17, непонятно в какую его часть, а может 11.14, там еще с десяток походу главы...
Каким образом найти корни?! Похоже вообще отдельная тема.
Кроме того корни ищутся относительно чего то, искать когда левая часть равна нулю или еще чему то? Как я буду искать корни вообще без левой части?
Извините я не академик, а вы говорите загадками. Поэтому и написал, что очевидное Вам не очевидно тому кто первый раз это видит.

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

и когда вы таким повелительным тоном требуете от меня

Где Вы такое вычитали? Даже в мыслях такого не было.
Если я кого то чем то обидел, прошу прощения, совершенно не хотел этого делать.

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

А ещё большим продвижением будет осознание, что мы с ilghiz говорим совершенно об одном и том же.

Да вижу я вижу, но слишком тут до много нужно догадываться.

Догадался, спасибо svv, полином в моем случае будет выглядеть так:
$a_0z^{100}+a_1z^{99}+...+a_{99}z^{1}+a_{100}z^{0}=0$
Все же это формула 11.14.

[B@R5uk]

Поигрался я тут с этим методом Прони:

Используется синтаксис Matlab M

clc
clearvars
format compact

if 0
    num = 700;
    filename = 'fragment_1.wav';
    [wdata, fs, nbits] = wavread (filename);
    %wdata = wdata (2000 : end - 10000);
    wdata = wdata (5001 : 10000);
else
    num = 60;
    ww = 2 * pi * (1 : 2000)';
    wdata = 0.5 * sin (ww / 39.3 + 0.2) + 0.3 * sin (ww / 63.1 + 0.1) + 0.1 * sin (ww / 98.1 + 0.5);
    wdata = wdata + 0.001 * randn (size (ww));
end
dsize = size (wdata, 1);

A = repmat (wdata, num, 1);
A = reshape (A (1 : end - num), dsize - 1, num);
A = A (num : end, :);
X = (A' * A) ^ -1 * A' * wdata (num + 1 : end);
Y = roots ([-1; X]);
B = repmat (Y', dsize, 1) .^ repmat ((1 : dsize)', 1, num);
Z = (B' * B) ^ -1 * B' * wdata;
wdata2 = real (B * Z);

%disp (Y)

subplot (211)
plot (wdata)
hold on
%plot (wdata2, 'g')
plot (wdata - wdata2, 'r')
grid on
hold off
subplot (212)
stem (2 * pi ./ imag (Y), abs (Z))
%stem (imag (Y) / (2 * pi), abs (Z))
%ylim ([0, 0.30000000001])
grid on
 

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

-- 28.09.2022, 16:42 --

Более того, модель

$$f(n)=\sum\limits_{k=1}^N\left(A_k\cos(\omega_k n)+B_k\sin(\omega_k n)\right)$$

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

$$f(n)=\sum\limits_{k=1}^N\left(A_k(n)\cos(\omega_k n)+B_k(n)\sin(\omega_k n)\right)$$

где $A_k(n)$ и $B_k(n)$ — это (относительно) медленно меняющиеся амплитуды. Это, правда, делает задачу недоопределённой, в связи с чем принятие решения о том, как именно разрезать спектр на гармоники становится эмпирической задачей для машинного обучения.

[Schrodinger's cat]

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

Жутко неустойчивая штука

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

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

становится эмпирической задачей

И как вы ее решаете если не секрет?
Метод отжига?

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