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

[svv]

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

Если компонент больше появляются нулевые амплитуды, тут понятно, можно отсеять.

Отсеивайте не только нулевые, но и очень малые на всём отрезке.

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

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

Это можно сделать ещё до составления и решения второй СЛАУ.
Так как у вас сигнал вещественный, характеристический полином будет иметь вещественные коэффициенты. То есть первая СЛАУ решается в вещественных числах.

На втором этапе Вы решаете полиномиальное уравнение. Каждый его корень $r_k$ либо вещественный, либо составляет пару с комплексно сопряжённым корнем $\overline{r}_k$. Допустим, функция-решатель этого уравнения перечисляет найденные корни вразброс. Я бы наводил порядок так. Пройдитесь по списку, ища корни с $\operatorname{Im}r_k>0$. Найдя такой $r_k$, внесите в новый список его и сразу за ним $\overline{r}_k$. Этот последний тоже где-то есть в списке, но специально искать его не надо, поскольку его значение и так известно. После этого внесите в новый список корни с $\operatorname{Im}r_k=0$, т.е. вещественные. Теперь в новом списке будут сначала идти комплексно сопряжённые пары (и известно, сколько их), а потом одиночные вещественные корни. На последующих этапах сохраняйте этот порядок.

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

Получились 2 амплитуды $+10$, $+10$ и начальные фазы $+2.642$, $+0.5$.

Понятно, точное значение первой фазы равно $\pi-0.5$.
Если сравнить вот с этими фазами:

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

Это и есть амплитуды (комплексные, т.е. учитывающие и фазу). Их можно записать и так:
$\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}$
В показателе экспонент в скобках — фазы.

— то Ваши больше на $\pi/2$. Так может получиться, если считать, что нулевую фазу имеет не косинус, а синус. Кто нарушает стандарт, Вы или Mathcad?

[Schrodinger's cat]

Пардон, само решение то не написал второго СЛАУ:
При заказе 2х компонент получается $h=4.794-i 8.776$, $4.794+i 8.776$

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

Кто нарушает стандарт, Вы или Mathcad?

А что я? я тут ни при чем
Похоже что я, у меня сигнал задан просто 0.5 без $\frac{\pi}{2}$

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

уже бы давно радовался бы

Метод замены первого этапа Прони на SVD и правый сингулярный вектор работает. Радуюсь. Спасибо.
Метод решения первого СЛАУ через МНК с QR (мою реализацию Грамм-Шмадта) вида $x=R^{-1}Q^TB$ Работает. Радуюсь. Спасибо.
Думаю, что 3 шага как Вы написали тоже должен работать, сейчас попробую, просто руки не дошли.
Через Холецкого еще не попробовал.
Через полином в целом разобрался, благодаря svv. Спасибо. Радуюсь.

Из непонятного у меня остается:
1) Вычисление спектра. По формулам Марпла получаются какие то мутанты. Да и там 2 графика, у меня только один.
2) Различия метода для затухающих и не затухающих экспонент. Тут мне прежде самому надо почитать, поэкспериментировать.
3) Там есть еще необязательный четвертый этап, оценивание направления прихода (пеленга) сигнала и определение резонансных пиков в диаграммах эффективной площади рассеяния радиолокационных целей. Стр. 391, самый конец. Но пока непонятно интересно это мне или нет. Непонятно в каком виде это получается. Будет ли это работать например со звуковым сигналом. Понятно, что постановка абсурдна, но если это какой то вид представления сигнала, то возможно это интересно для машинного обучения, возможно эта форма будет более ясна для нейронных сетей чем исходный сигнал.
Это вот как идентифицируют типы целей радары и головки самонаведения?
А это можно применить скажем для идентификации патернов в сигнале? Или это не о том?

[B@R5uk]

Вот интересно, как интерпретировать скачки в положениях максимума автокорреляционной функции?



Понятно, что это высокие гармоники шалят, когда у низких в амплитуде провал. Но по идее сигнал то один, и он меняется по частоте (и периоду) плавно. В чём же дело? Ошибка в периоде на величину скачка вызовет уже в 20-х гармониках сдвиг индекса на единицу. Забавно, так же, что расстояние между скачками равно периоду сигнала.

study_cycle.m

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

clc
clearvars
format compact

cnum = 600;
step = 30;
csize = 750;

[wdata, fs, nbits] = wavread ('Mi.wav');
dsize = numel (wdata);

dnum = ceil (dsize / step);

env = 1 - cos (2 * pi * (0.5 : csize)' / csize);
correl = zeros (cnum, dnum);
for k = 1 : cnum
    tmp = [wdata .* [wdata(k : end); zeros(k - 1, 1)]; zeros(step * dnum + csize - dsize, 1)];
    for l = 1 : dnum
        correl (k, l) = sum (tmp ((1 : csize) + (l - 1) * step) .* env);
        %correl (k, l) = sum (tmp ((1 : csize) + (l - 1) * step));
    end
end
maxpos = zeros (dnum, 1);
for k = 1 : dnum
    correl (:, k) = correl (:, k) / correl (1, k);
    maxval = max (correl (220 : 500, k));
    maxpos (k) = find (correl (:, k) == maxval, 1, 'first');
end

correl = min (1, max (-1, 0.85 * correl));

subplot (111)
imshow (round (255 * (1 + correl)) / 2, jet, 'InitialMagnification', 50)
axis on
axis xy
%image (round (255 * (1 + flipud (correl))) / 2)
hold on
plot (maxpos, '.k')
hold off

%surf (correl)
%surf (log (abs (correl)))
%shading interp
%axis vis3d
 

[Schrodinger's cat]

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

Это и есть амплитуды

$a_1={\color{red}-}10e^{(-\frac{\pi}2+0.5)i}$
$a_2=+10e^{(+\frac{\pi}2-0.5)i}$
Тут кажется минус должен быть.

$z=0.978+i 0.208$ и $0.978-i 0.208$
Задержки $+i 0.209$ и $-i 0.209$
Частоты $+0.209$ и $-0.209$

$h=-8.776-i 4.794$ и $-8.776+i 4.794$
Амплитуды $10$ и $10$
Частоты $+0.5$ и $-0.5$

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

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

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

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

Работает. Спасибо. Радуюсь.
Только Q прямоугольная. Но все равно транспонированная ей дает ответ.

Холецкий в общем тоже работает. Сделал не до конца, но очевидно, что работает. Спасибо. Радуюсь.

[ilghiz]

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

\
Холецкий в общем тоже работает. Сделал не до конца, но очевидно, что работает. Спасибо. Радуюсь.

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

[svv]

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

$a_1={\color{red}-}10e^{(-\frac{\pi}2+0.5)i}$
$a_2=+10e^{(+\frac{\pi}2-0.5)i}$
Тут кажется минус должен быть.

Спасибо за бдительность, но я проверил в Вольфраме, он подтвердил, что всё правильно:
10*e^((-pi/2+0.5)*i)*e^(2*pi*i/30*n) + 10*e^((pi/2-0.5)*i)*e^(-2*pi*i/30*n) = 20*sin(2*pi/30*n+0.5)
Просто я записал $a_1$ и $a_2$ в показательной форме, избавившись от множителей $\pm i$ с помощью формулы $\pm i=e^{\pm\frac{\pi}2i}$. Так сказать, перевёл эти множители в фазу.

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

При заказе 2х компонент получается $h=4.794-i 8.776$, $4.794+i 8.776$

И это правильно.
Точные значения $10e^{(-\frac\pi 2+0.5)i}$ и $10e^{(\frac\pi 2-0.5)i}$. То есть как раз $a_1$ и $a_2$. Согласно "моим" определениям,

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

В показателе экспонент в скобках — фазы.

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

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

Вот интересно, как интерпретировать скачки в положениях максимума автокорреляционной функции?

Понятно, что это высокие гармоники шалят, когда у низких в амплитуде провал. Но по идее сигнал то один, и он меняется по частоте (и периоду) плавно. В чём же дело? Ошибка в периоде на величину скачка вызовет уже в 20-х гармониках сдвиг индекса на единицу. Забавно, так же, что расстояние между скачками равно периоду сигнала.

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

[Schrodinger's cat]

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

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

В общем то нет, решив всеми методами просто потренировался.
Но вопрос, а чем же решают на практике, актуален. Какой вариант посоветуете отправлять в C#?

[Schrodinger's cat]

Перечитал 11.5 "Метод наименьших квадратов Прони".
1) Правильно ли я понял, что незатухающая гармоника может быть представлена исходным методом Прони суммарно за 2 комопненты (собственно мой пример, что мы разбирали) а модифицированным получается может за одну?
2) Я так понял, что для модифицированного полином удлиняется в 2 раза. А для СЛАУ надо градиенты составлять? Как правильно составить СЛАУ?
3) Как сделать искусственную комплексную компоненту шума? Для моего примера, чтобы он с шумом был.

По поводу спектра. Вот так же нужно:

$$X_1(z)=\sum_{k=0}^{p-1}(\frac{h_k}{1-z_kz^{-1}})$$

$$S_1(f)=\left \| TX_1 e^{i 2\pi fT} \right \|^2$$

T - Период вообще нужен в этой формуле?
Пределы спектра $-1/2T$ и $1/2T$?
Если выбрать такие то ничего не видно, а если больше то картина повторяется.

[ilghiz]

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

Но вопрос, а чем же решают на практике, актуален.

От задачи зависит. Вот, например, если вы можете исходный вектор в памяти держать, и памяти на QR тоже есть, то QR разложение исходной матрицы с ранк-ревеалингом - это самый устойчивый и численно быстрый алгоритм.

Но, например, у меня часто возникает ситуация, когда данные сыпятся на ПЛИС, и у ПЛИС нет возможности их всех сохранять, а вот на матрицу для Холецкого памяти есть, тогда я решаю через Холецкого, вернее конечно "правоверного" Холецкого с выбором ведущего блочного элемента и уже не на ПЛИСке, но его в 10 строк не реализуешь (у меня около 3 тысяч строк + комментарии на С++).

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

Какой вариант посоветуете отправлять в C#?

в три шага как я уже советовал - он примерно в 30 строк кода на C# уместится, а проблем с памятью наверное у вас не будет, так как C# обычно на МК или ПЛИСках не ходит.

[Schrodinger's cat]

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

От задачи зависит.

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

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

проблем с памятью наверное у вас не будет, так как C# обычно на МК или ПЛИСках не ходит.

Именно.

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

(у меня около 3 тысяч строк + комментарии на С++)

Сколько простите? Это полностью класс со всякими плюшками? С трудом представляю себе функцию в 3к строк, хотя я видывал не мало (под ПК, с МК опыта довольно мало, но тоже имеется). Да и алгоритм не такой сложный.

ПС: Вспоминается ассемблер, времена когда выжимали лишние байты и скорость.

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

Кто нарушает стандарт, Вы или Mathcad?

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

Похоже что я, у меня сигнал задан просто 0.5 без $\frac{\pi}{2}$

Обнаружил у себя ошибку. Похоже я чист, это все mathcad

[ilghiz]

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

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

(у меня около 3 тысяч строк + комментарии на С++)

Сколько простите? Это полностью класс со всякими плюшками? С трудом представляю себе функцию в 3к строк, хотя я видывал не мало (под ПК, с МК опыта довольно мало, но тоже имеется).

конечно класс, но из плюшек там только Холецкий - одна функция для разложения, и одна - для решения, ничего лишнего. В самом классе конечно есть сколько-то функций - где-то около сотни функций и лямбд.

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

Да и алгоритм не такой сложный.

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

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

С трудом представляю себе функцию в 3к строк, хотя я видывал не мало (под ПК, с МК опыта довольно мало, но тоже имеется).

Исходники Lapackа посмотрите. Ах, да, у меня тоже есть где-то на 5 тысяч строк ОДНА функция, правда с несколькими десятками лямбд. Это очень распространенный стиль программирования в вычислительной математике, ибо на входе - данные, на выходе - тоже данные, и смысла заводить под это все класс обычно нет, а, в некоторых случаях удобно все через лямбды вытащить, но, наверное, в этой теме это оффтопик.

[ilghiz]

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

ПС: Вспоминается ассемблер, времена когда выжимали лишние байты и скорость.

так и сейчас есть куча задач, где надо, чтобы все до последнего флопа было, просто на C# это никто не делает, да и не на ассемблере сейчас эта битва идет, а на уровне правильного распределения данных по памяти, чтобы в кеши легло, по ядрам распараллелилось и все конвейеры заполнило. Вот возьмите хорошо оптимизированную лапаковскую умножаловку двух матриц, например из пакета Intel Math Kernel Library, и в три цикла напишите свою и сравните на первом попавшемся более-менее современном компе эти две умножаловки, например на квадратных матрицах размера $5000 \times 5000$, гарантирую вам хотя бы в 10 раз разницу по скорости. А если вы на CUDA полезете, то там разница между наивным и обычным методом, и хорошо оптимизированным, часто бывает около миллиона раз.

[Schrodinger's cat]

(Оффтоп)

ilghiz
Не верю, знаю, проверял недавно один свой алгоритм с применением PCA.
SVD брались из разных либ. Результаты:
- AlgLib 12.916 сек
- NumPy (C#) 5.061 сек
- MathNet (C#) 2.618 сек
- MathNet (MKL) 0.654 сек
Своих реализаций SVD у меня тогда еще не было, не думаю что смог бы победить MKL.
Если сравнивать самую худшую реализацию и MKL выигрыш почти в 20 раз.

Сам давно пробовал оптимизировать всякие функции, в основном математического характера, и быстро понял, что современные компиляторы настолько умны и они выдают настолько хороший код, что и пробовать обогнать их бессмысленно.
Смотрел одну конференцию, там человек тоже заморочился оптимизацией и к своему удивлению обнаружил, что его код в 10 раз хуже. Он не сдался и потратил неделю на это дело. Итоговый выигрыш у него составил 8% (быстрее чем сделал компилятор).
Поэтому все эти ручные ассемблерные оптимизации для десктопных программ давно закончились, просто нужно установить оптимизацию в настройках компилятора и он все сделает.
А вот написать изначально годный алгоритм никто не отменял.

На CUDA сам не писал, но видел их рекламный код, они кажется сумму матриц делали, выигрыш в 500 раз.
Но использую TensorFlow (фреймворк для нейронных сетей) при переключении с CPU на GPU выигрыш не измерял, но он очевиден.

Когда я только начинал знакомиться с Python и TensorFlow мне понадобилось перенести логику программы с C# на Python еще даже до завершения этой процедуры предварительные тесты показывали то, что работает на C# минуты, десятки минут, часы, на классической реализации Python работает часы, дни, недели.

Поэтому с Вашими выводами по скорости полностью согласен.

[ilghiz]

Про MKL - да, на десктопах он должен всех обыгрывать, хотя его иногда слегка обыгрывают ATLAS и AOCL.

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

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

Вот же Фома не верующий!!! Это в корне не верное утверждение. Компилятор не умеет и не имеет право менять порядок операций, или менять место расположения данных в алгоритме.

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

У меня эти два варианта на моем компе с максимальным уровнем оптимизации отличаются по скорости выполнения в 35 раз. Если умножать блочно, то быстрая реализация "проиграет" правильной блочной еще в 20 раз, но вот блочная реализация уже в 20 строк не вместится, то есть слегка кривая реализация алгоритма будет в 700 раз медленнее более-менее правильной, и это все без ассемблера (чтобы посмотреть что из себя представляет блочная реализация, можно поискать исходники dgemm в ATLAS-BLAS)!

Уверен, что и на C# будет также. Может размерность матрицы только надо немного по больше сделать, чтобы сильнее почувствовать разницу.

Используется синтаксис C++

#include <stdio.h>

#define N 2048

double A[N][N], BT[N][N], C[N][N];

int main()
{ for(int i=0; i<N; i++)
    for(int j=0; j<N; j++)
      A[i][j]=BT[i][j]=C[i][j]=i+j; // нужно только чтобы проинициализировать
  for(int i=0; i<N; i++)
    for(int j=0; j<N; j++)   // вариант 1 - оставить как есть
      for(int k=0; k<N; k++) // вариант 2 - поменять местами циклы по j и k - ведь тут никакого ассемблера и все должно быть точно также?
        C[j][i]+=A[j][k]*BT[i][k];
  double s=0;
  for(int i=0; i<N; i++)
    for(int j=0; j<N; j++)
      s+=C[i][j]; // нужно только чтобы компилятор не выбросил то, что выше
  printf("%lg\n", s);
  return 0;
}