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

[ilghiz]

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

Чудеса математики

это не чудеса математики, а кривизна вашей реализации. Той нестабильности, о которой вы писали выше быть не должно от слова совсем.

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

[svv]

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

Вектор $c = 0.671, -0.314, -1.285$

Это правильно! После этого надо решить полиномиальное уравнение
$r^p+c_{p-1}r^{p-1}+...+c_1r+c_0=0$,
или, что то же самое,
$\sum\limits_{k=0}^p c_k r^k=0$,
где по определению $c_p=1$.
Вы его решили правильно и $\lambda_k$ тоже нашли правильно (кроме, кажется, знака в одном месте).

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

Первое значение непонятно что

Если бы Вы задали $p=2$, вектор $c$ получился бы: $c_0=1,\; c_1\approx -1.956295$ (точное теоретическое значение $c_1=-2\cos\frac{2\pi}{30}$).
Составляем характеристическое уравнение: $r^2-1.956295r+1=0$.
Его корни $0.97815 \pm 0.20791 i$ (точные значения $e^{\pm\frac{2\pi}{30}i}$).
Отсюда "лямбды" равны $\pm 0.20944 i$$ (точные значения $\pm\frac{2\pi}{30}i$).

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

Вот что произошло в Вашем конкретном случае. На уровне ошибок округления метод выявил экспоненту, которая аннулируется оператором $\mathsf L+0.67145$. Число $0.67145$ приближённое и объяснить его очень трудно. Но зато, зная его, можно объяснить всё остальное. Из сказанного выше ясно, что две "настоящие" экспоненты аннулируются оператором
$(\mathsf L-0.97815 -0.20791 i)(\mathsf L-0.97815 +0.20791 i)=\mathsf L^2-1.956295\mathsf L+1$
Значит, все три экспоненты будут аннулироваться оператором
$(\mathsf L+0.67145)(\mathsf L^2-1.956295\mathsf L+1)\approx\mathsf L^3 - 1.285 \mathsf L^2 - 0.314 \mathsf L + 0.671$
А это и есть компоненты Вашего вектора $c$.

Для этой шумовой экспоненты получаем $\lambda=\ln(-0.67145)\approx -0.398+3.142i$ (вот тут, мне кажется, Вы потеряли знак у вещественной части). Ясно, что мнимая часть обусловлена знаком "минус" числа $-0.67145$, и её точное значение равно $\pi$.

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

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

И кстати говоря, вопрос: Все сверточные разложения типа Хартли, Косинусное,... , будут обладать теми же недостатками что и Фурье?
Z, Лапласа сюда можно отнести или это из другой оперы?

Хартли - способ посчитать Фурье, если на комплексные числа аллергия ;)
Но по сути - тот же Фурье (ещё может быть небольшой выигрыш в скорости счёта, если сравнивать с обычными реализациями Фурье, но и реализации Фурье можно оптимизировать, так что выигрыш не гарантирован).
Косинусное - Фурье, если функция представляется чётная (изначально чётная, или мы искусственно дополняем симметричным отрезком для чётности).
То есть все недостатки и достоинства Фурье в наличии, частные удобства для отдельных приложений могут быть важны, но принципиально ничего не меняют.
Лаплас - в общем, та же опера. Двустороннее преобразование Лапласа это Фурье с комплексным аргументом.

[Schrodinger's cat]

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

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

Возможно, какая-то стандартная функция в Mathcad решает (в некотором смысле) и переопределённые системы. Если нет, я скажу, что делать. Пока меня интересует только решение этой системы

Озвучьте пожалуйста, если нет, что тогда? Когда я захочу это переложить на C# у меня не будет волшебной функции mathcad'a. Неплохо в этом случае знать варианты.

Еще смущает то, что lsolve, решатель СЛАУ mathcad'a, который получил этот правильный ответ, как написано в справке основан на LU, а как оно может разложить прямоугольную матрицу?! Скорее всего (я это вечером проверю) он просто берет верхушку прямоугольной матрицы, раскладывает, а весь низ зануляется. В этом смысле таки получается точный ответ но для первых строк (уравнений).

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

Но Вы "заказали" не две, а три экспоненты.

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

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

Но это не страшно, так как у этой экспоненты на последнем этапе метода (до которого Вы ещё не дошли) получится очень малый "амплитудный" коэффициент, на уровне точности вычислений. По этому признаку, как я писал выше, она просто отбрасывается.

На самом деле я уже пробовал решать и вторую СЛАУ, не знаю насколько верно я ее составляю, но амплитуды получились 0, 10, -10 (на память знаки могу путать), начальные фазы не помню, но число 0.209 там есть. Они похожи на задержки. Так и должно быть?
Так же пробовал строить спектр, почти ничего не получилось.
И допер, что:
$x_n=\sum\limits_{k=0}^{p-1} h_k z_k^n$
$z$ - это корни первого этапа.
$h$ - результат второго этапа.
Строит аппроксимацию и экстраполяцию. Вроде работает.
Смущает то, что z в степени n, из-за длинного ряда и длинной экстраполяции это число возможно может выходить за пределы регистров.
Но! Это все забегания вперед, хочется основательно сперва разобраться с текущими проблемами.

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

а кривизна вашей реализации.

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

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

Не в обиду.

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

[ilghiz]

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

Озвучьте пожалуйста, если нет, что тогда? Когда я захочу это переложить на C# у меня не будет волшебной функции mathcad'a. Неплохо в этом случае знать варианты.

Так я же вам уже писал как все вручную написать! Грамм-Шмидт, одно умножение матрицы на вектор и одно решение линейной системы в верхней треугольной матрицей. Каждое из этого пишется в 5-10 строк на любом языке программирования.

То есть если у вас есть задача

$$\min_c ||Ac - x||_2^2$$

вам надо выполнить только следующие шаги
1. методо Грамма-Шмидта разложить матрицу $A$ на ортогональную $Q$ и верхнюю треугольную $R$ как $A = QR$,
2. вычислить $b = Q^* x$,
3. найти решение линейной системы с верхней треугольной матрицей вида $Rc = b$ методом обратной подстановки.

Если у вас при очередном шаге Грамма-Шмидта получается почти нулевая норма, то вы взяли слишком много сдвигов (то есть столбцов) в исходной матрице и надо просто урезать саму матрицу. Алгоритм будет полностью детерминистический и всегда будет давать предсказуемые результаты.

[B@R5uk]

Вот интересно, можно ли построить аналог Метода Прони, который будет искать чисто комплексные экспоненты, без затухания/нарастания?

[Schrodinger's cat]

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

без затухания/нарастания?

Да. Второй вариант как раз для НЕ затухающих экспонент.
В моем вымышленном примере они, я так понимаю, не затухающие. В прогнозе по крайней мере не затухают.

ПС: Похоже, что метод Прони по варианту ilghiz называется matrix pencil method (MPM), случайно нашел.
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2473-y
Пишут, что он менее чувствителен к шуму чем полиномиальный.

[svv]

B@R5uk, да. В книге Марпла-младшего он называется модифицированным методом Прони и описывается в пункте 11.6 на стр.376:

Цитата:

Обычный метод наименьших квадратов Прони может быть модифицирован для аппроксимации последовательности комплексных данных с помощью модели, состоящей из незатухающих ($\alpha=0$) комплексных синусоид. Ниже мы рассмотрим лишь случай четного числа комплексных экспонент.

-- Вт окт 04, 2022 13:35:56 --

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

Озвучьте пожалуйста, если нет, что тогда?

Мне нечего добавить к тому, что уже рассказал уважаемый ilghiz.

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

Мне казалось, что все осмысленное должно "всплывать" на первые места, а все бессмысленное на последние.

Да, тут так не происходит, по двум причинам. Первая: при известных частотах (точнее, комплексных коэффициентах нарастания/затухания) рекуррентное соотношение никак не зависит от амплитуд, а следовательно, и корни характеристического полинома тоже. Вторая: порядок перечисления этих корней функцией, которая их нашла, ничем не регламентирован, и зависит от реализации.

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

На самом деле я уже пробовал решать и вторую СЛАУ, не знаю насколько верно я ее составляю, но амплитуды получились 0, 10, -10 (на память знаки могу путать), начальные фазы не помню, но число 0.209 там есть. Они похожи на задержки. Так и должно быть?

Похоже. Амплитуда той непонятной экспоненты $0$, и о ней можно забыть. А $\pm 10$ там, конечно, есть. Но надо ещё уточнить, что именно получилось, и проверить.
Я тут, да и раньше, использую слово "амплитуда" в смысле "коэффициент перед экспонентой".

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

Правильно так:
$\begin{array}{ll}a_1=-10ie^{+0.5{\color{magenta}i}}\\[0.8ex]a_2=+10ie^{-0.5{\color{magenta}i}}\end{array}$

Это и есть амплитуды (комплексные, т.е. учитывающие и фазу). Их можно записать и так:
$\begin{array}{ll}a_1=10e^{(-\frac{\pi}2+0.5)i}\\[0.8ex]a_2=10e^{(+\frac{\pi}2-0.5)i}}\end{array}$
В показателе экспонент в скобках — фазы.

[B@R5uk]

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

В книге Марпла-младшего...

А можно полностью название книги с инициалами автора?

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

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

Если на промежутке данных, на 30...50% более длинном, чем период первой гармоники, сделать разложение по этой модели, то получаются довольно неплохие результаты. Если это проделать для последовательности перекрывающихся промежутков, а потом (с помощью огибающей) их гладко сшить, то получится довольно хорошее разложение сигнала с меняющейся частотой и амплитудой на гармоники. Основная проблема: угадать частоту на каждом промежутке. В эксперименте для одного промежутка я делаю перебор в лоб:



На первом графике — отклонение (в определённом смысле) аппроксимации по модели от исходных данных, на втором — исходные данные (зелёная линия), аппроксимация (синяя линия), оптимальная на данный момент аппроксимация (чёрная линия) и первая гармоника (красная линия). В разложении используется 40 гармоник и константа. Если бы только можно было найти правильную частоту $\omega=0,01866$ без перебора!



Mi.wav (48 кГц, 56 КБ)

study_freq_3.m

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

clc
clearvars
format compact

tnum = 40;
slen = 600;
spos = 15800;
w1 = round (1.3 * 2000 * pi / slen) / 1000;
w2 = 0.022;
wnum = 1001;
frmstep = 20;

[wdata, fs, nbits] = wavread ('Mi.wav');
sdata = wdata (spos + (0 : slen - 1));

ww = linspace (w1, w2, wnum);
vals = zeros (wnum, 1);
nn = (1 : slen)';
pp = repmat (nn, 1, tnum) .* repmat (1 : tnum, slen, 1);
env = 1 - cos (2 * pi * (nn - 0.5) / slen);
bval = Inf;
fh = gcf ();
frmdata = getframe (fh);
frmdata = frmdata .cdata;
imdata = zeros ([size(frmdata, 1), size(frmdata, 2), 1, ceil(wnum / frmstep)]);
map = [0.8 0.8 0.8; 0 0 0; 1 1 1; 1 0 0; 0 1 0; 0 0 1];
%map = [0 0 0; 0.8 0.8 0.8; 1 1 1; 0 1 0; 0 0 1; 1 0 0];
for l = 1 : wnum
    B = [ones(slen, 1), cos(pp * ww (l)), sin(pp * ww (l))];
    X = B \ sdata;
    adata = B * X;
    val = sqrt (sum ((sdata - adata) .^ 2));
    vals (l) = val;
    if bval > val
        bval = val;
        bdata = adata;
        bl = l;
    end
    harm1 = B (:, 2) * X (2) + B (:, 2 + tnum) * X (2 + tnum);
    
    subplot (211)
    plot (ww (1 : l), vals (1 : l))
    hold on
    if 1 ~= l
        plot (ww (bl), vals (bl), 'sk')
        plot (ww (l), vals (l), 's')
    end
    hold off
    grid on
    xlim ([w1, w2])
    ylim ([0, 4])
    
    subplot (212)
    plot (sdata, 'g', 'LineWidth', 2)
    hold on
    plot (bdata, 'k')
    plot (adata)
    plot (harm1, 'r')
    hold off
    grid on
    ylim ([-0.4, 0.6])
    
    if 1 == mod (l, frmstep)
        frmdata = getframe (fh);
        frmdata = rgb2ind (frmdata .cdata, map);
        imdata (:, :, 1, 1 + (l - 1) / frmstep) = frmdata;
    end
    
    pause (0.01)
end

imwrite (1 + imdata, map, 'study_freq_3.gif', 'GIF', 'LoopCount', Inf, 'TransparentColor', 1)

subplot (211)
plot (ww (1 : l), vals (1 : l))
hold on
plot (ww (bl), vals (bl), 'sk')
hold off
grid on
xlim ([w1, w2])

subplot (212)
plot (sdata, 'g', 'LineWidth', 2)
hold on
plot (bdata, 'k')
hold off
grid on
 

[svv]

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

С. Л. Марпл-мл. Цифровой спектральный анализ и его приложения. Глава 11.

[B@R5uk]

svv, спасибо. Скачал, будем почитать.

Я правильно понимаю, что отсутствие затухания/нарастания в экспонентах приводит к симметричности коэффициентов порождающего полинома?

[svv]

Да, но с оговоркой. Итак, у Вас вещественный сигнал, а все экспоненты имеют вид $e^{i\omega n}$. Потребуем для однозначности $\omega\in(-\pi,\pi]$. Разобьем экспоненты в зависимости от $\omega$ на три группы:
1) $\omega=0$
2) $\omega=\pi$
3) все остальные.
Экспоненты третьей группы могут входить в сигнал только парами $ae^{i\omega n}+\bar ae^{-i\omega n$. Если все экспоненты в сигнале таковы, то да, коэффициенты характеристического полинома симметричны: $c_k=c_{p-k}$.
Экспонента c $\omega=\pi$ даёт вклад в полином в виде множителя $r+1$, который не нарушает симметричность.
Но вот экспонента с $\omega=0$ даёт вклад в виде множителя $r-1$, который меняет симметричность на антисимметричность, $c_k=-c_{p-k}$.

[Schrodinger's cat]

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

кривизна вашей реализации

Перепроверил кривизну моей реализации метода Грамм-Шмидта и выяснил что на тестовых данных, матрицах любой размерности, пробовал так же как в методе до 100 строк, 2...10 столбцов, ничего не меняется, разложение стабильно.
А вот применительно к методу Прони, там где мы заказываем излишние компоненты и он их вынужденно изобретает, вот в этой части наблюдается некий разброс. При этом компоненты с осмысленной частью стабильны.

Решение первого СЛАУ через SVD, при заказе компонент больше необходимого, получается неверным.
Видимо он дополнительно разлагает компоненты на большее число компонент.

Решение первого СЛАУ mathcad'овским QR невозможно по причине того, что он не может вычислить QR.

Но зато специализированная функция lsolve работает, стабильно и правильно, ей пока и буду пользоваться. А нюансы переноса на C# учту потом.

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

вот тут, мне кажется, Вы потеряли знак у вещественной части

Потерял, пока писал. Прошу прощения.

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

Но Вы "заказали" не две, а три экспоненты

Интересно так же поведение когда мы заказываем меньше компонент чем есть. Решение так же получается не совсем ожидаемым.
Откуда возникает вопрос, а как это программно анализировать?
Если компонент меньше чем есть амплитуда не 0.
Если компонент больше появляются нулевые амплитуды, тут понятно, можно отсеять.
Как правильно отлавливать компоненты образующие пары? Тупо по совпадению модулей их аплитуд, частот...

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

Если бы Вы задали $p=2$, вектор $c$ получился бы: $c_0=1,\; c_1\approx -1.956295$

Получается.

Первое СЛАУ создаю как написали Вы, это несколько отличается от Марпла.
Второе СЛАУ я создаю по Марплу, оно решается, получаются амплитуды и начальные фазы.

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

Но надо ещё уточнить, что именно получилось, и проверить.

Получились 2 амплитуды $+10$, $+10$ и начальные фазы $+2.642$, $+0.5$. При заказанных 2х компонетах.
Если заказать больше, то результат такой же. Могут только местами меняться.

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

Еще смущает то, что lsolve, решатель СЛАУ mathcad'a, который получил этот правильный ответ, как написано в справке основан на LU, а как оно может разложить прямоугольную матрицу?!

Разобрался. Он допускает прямоугольную L.

Готов пойти дальше.
Спектр по прежнему непонятно как вычислять.

[ilghiz]

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

Решение первого СЛАУ через SVD, при заказе компонент больше необходимого, получается неверным.

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

Я бы на вашем месте чем писать такие длинные отчеты тут, таки запрограммировал бы те три шага, что я написал (да хоть на C#) и уже бы давно радовался бы, там все в 20-30 строк кода помещается.

[B@R5uk]

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

(Оффтоп)

...Если бы только можно было найти правильную частоту $\omega=0,01866$ без перебора!

Помедитировал над анимацией (для разных участков сигнала). Обратил внимание на то, что подгоняемая модель (синяя и чёрная линии) на левом и правом концах обрабатываемого участка одинаковы (потому что период первой гармоники меньше длины этого участка) и являются усреднённым значением сигнала справа и слева (тоже очевидный результат). А минимум "дисперсии" наступает, когда эти усредняемые концы совместятся. Ну, дык! Это ж почти автокорреляционная функция на первом графике! (Не зря я запас по длине отрезка в 30...50% взял — интуиция). Её и перебрать в лоб куда проще, чем один раз коэффициенты модели рассчитать для фиксированного омега. И не надо никакие полиномиальные уравнения решать. Кажется путь намечен...