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

[Schrodinger's cat]

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

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

А может вам проще попользовать метод линейного предсказания?

Это какой? Типа Бёрг?
А смысл, когда у нас уже есть гармоники из разложения, которые уже можно экстраполировать?

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

Если у вас мало шума и только две гармоники, то вам достаточно основной сигнал $f_i, i=1,...,N$ пересохранить в матрицу $P = [p_1,...,p_5]$, $p_k = (f_k,...,f_{N-5+k})^T$ и у этой матрицы гарантированно будет ранг 4, а, соответствующий нулевому сингулярному значению правый сингулярный вектор будет сразу давать вам частоты с фазами, и все заморочки с тем, что ваша сетка будет не кратна вашим частотам уйдут в небытье, как и логарифмическая сложность от Фурье, то есть все будет линейно, с довольно маленьким коэффициентом.

Описание похоже на SSA (Singular Spectrum Analysis) типа Гусеницы.
Но это же совсем другая история
Да и там не все так однозначно с выделением отдельных гармоник. Главные компоненты могут быть сами по себе а могут (и часто это делают) образовывать коррелированные пары с другими компонентами. В написании программы по анализу каких либо данных пойди пойми это все еще одна гармоника или это уже другая... В большинстве случаев корреляция конечно показывает их связь, но иногда встречаются неоднозначные случаи.
Или я чего то не понял? Как может сингулярный вектор давать частоты с фазами?

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

[svv]

Метод Прони (ударение на «и»)
С. Л. Марпл-мл. Цифровой спектральный анализ и его приложения. Глава 11.
Метод Прони — это разложение сигнала на экспоненциальные составляющие. Коэффициенты затухания экспонент могут быть и комплексными, так что задачу выделения гармонических составляющих метод тоже решает. Они не обязаны быть кратны чему бы то ни было.
На первом этапе по сигналу находится характеристический полином (осуществляющий то самое линейное предсказание).
Второй этап — определение коэффициентов затухания, сводится к решению одного алгебраического уравнения высокой степени.
На третьем этапе определяются коэффициенты при экспонентах.
Я этим методом когда-то давно занимался, метод нормально работал.

[ilghiz]

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

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

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

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

Описание похоже на SSA (Singular Spectrum Analysis) типа Гусеницы.
Но это же совсем другая история
Да и там не все так однозначно с выделением отдельных гармоник. Главные компоненты могут быть сами по себе а могут (и часто это делают) образовывать коррелированные пары с другими компонентами.

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

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

Описание похоже на SSA (Singular Spectrum Analysis) типа Гусеницы.
Но это же совсем другая история
...
Как может сингулярный вектор давать частоты с фазами?

Раньше, в 90-е годы это называлось линейным предсказанием, но, по современной википедии - описание как для SSA.
Ссылки под рукой нет, да и 80-90-х годов она, ибо уже в 1999 я это во всю применял.

Там смысл такой:

Пусть у вас сигнал - $f_i$ одна комплексная экспонента, тогда
$a f_i = f_{i+1}$ можно выразить через $f_i$, и $a^{-1}$ - как раз комплексное число, в котором норма - амплитуда, а все экспоненциальная часть - частота.

Для двух комплексных экспонент надо взять уже $a f_i$ и $a^2 f_i$, и тогда при выражении двух сдвинутых векторов через исходный вы получите в качестве коэффициентов - коэффициенты полинома второй степени, корни которого соответствуют этим экспонентам.

То есть фактически, набор $f_i, f_{i+1}, f_{i+2},...$ умноженный на коэффициенты $a^k, a^{k-1}, ..., a, 1$ стремиться к нулю, когда $a^{-1}$ соответствует комплексным экспонентам, и как раз надо взять для $k$ экспонент полином $k$-ой степени.

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

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

-- 25.09.2022, 22:29 --

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

Метод Прони (ударение на «и»)

а, точно, значит оно самое!

[ozheredov]

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

Метод Прони

Слегка причастился. Есть вопросы:

1. Верно ли, что характеристический полином имеет степень порядка числа отсчетов сигнала?
2. И чем искать корни полинома 100500-й степени? Риторический вопрос -- будут ли корни устойчивыми? И нериторический -- а нам эта устойчивость нужна? (Интуитивно кажется, что нужна.)

[svv]

1. Нет, полином имеет степень, равную числу экспоненциальных составляющих вида $Ae^{\lambda t}$, которые нужно выявить в сигнале. Здесь и $A$, и $\lambda$ могут быть комплексными, поэтому метод фактически обеспечивает и поиск синусоидальных составляющих. Число отсчётов обычно избыточно, но это не увеличивает степень полинома (число отсчётов лишь ограничивает степень полинома сверху).

Вопрос "а сколько искать синусоид?" не такой простой. Идеальная ситуация — когда их количество в сигнале заранее известно и нет шумов. С шумами можно бороться, "заказав" поиск большего числа синусоид, чем реально есть в сигнале (скажем, искать 20 вместо 10 заведомо имеющихся). Синусоидальные составляющие шумового происхождения можно потом распознать по некоторым признакам (малая амплитуда и высокая частота) и отбросить. Если этого не делать, шумы будут снижать точность.

2. В книге Марпла приводится программа на Фортране для факторизации комплексного полинома методом Дженкинса и Трауба. Я её когда-то переписывал на C++, но сейчас уже вряд ли найду. Но должны быть библиотеки, кто-то добрый мог реализовать факторизацию (возможно, другим методом) и выложить на GitHub. Может, подойдёт и матпакет.

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

[ilghiz]

Добавлю свои пять копеек.

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

Сколько я пробовал решать, больше 30-50 экспонент там все обычно упирается в тот факт, что надо
1. мало шума,
2. точное понимание числа этих экспонент,
и, из-за этого метод начинает сильно затыкаться и дает плохое решение.

Еще у метода есть положительное свойство, если экспонента во времени слегка "плывет", это все еще позволяет решать эту задачу, в то время как Фурье сильно размывает пик и теряет его среди шума.

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

[Schrodinger's cat]

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

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

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

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

Первое за что зацепился глаз, так это под конец главы предлагается перейти к SVD разложению, что делает метод крайне похожим на Анализ Главных Компонент (АГК, он же PCA), SSA и подобных. Только тут в голом виде SVD не нормированной и не центрированной матрицы. Матрица Теплица, в АГК любят Ганкеля, что по сути не особо важно. С этими методами я уже разбирался.
Но неясно, как имея SVD разложение определить амплитуды и фазы компонент?

В книге речь все про комплексный сигнал. А у нас обычно имеется вещественный. Как быть?

Там упоминаются всякие процедуры типа разложения Холецкого, решение СЛАУ, соответственно разложение SVD это все я так понимаю нужно иметь в варианте для комплексных чисел?

Цитата:

Есть еще один метод

Спасибо. Попозже воспользуюсь Вашим предложением. Сперва с этим разберусь.

[ilghiz]

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

Но неясно, как имея SVD разложение определить амплитуды и фазы компонент?

я тут распинался, формулы выводил, рассказывал, интересно, для кого я это делал?

[B@R5uk]

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

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

А вы пробовали с помощью Фурье (а точнее, БПФ) очищать реальный звук от шумов? Вот, например, вот так:

fragment_edit_3.m

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

clc
clearvars
format compact

filename = 'fragment';
blocksize = 512;
fade = 20;
val_a = 2.8;
val_b = 1.5;
blur_map = ones (5, 3) / 15;
width = 7;
height = 91;

[wdata, fs, nbits] = wavread ([filename, '.wav']);
dsize = size (wdata, 1);
bnum = floor (dsize / blocksize) - 4;
env = 0.5 * (1 - cos (2 * pi * ((1 : 4 * blocksize)' - 0.5) / (4 * blocksize)));

spec = zeros (4 * blocksize, bnum);
for k = 1 : bnum
    spec (:, k) = fft (env .* wdata ((k - 1) * blocksize + (1 : 4 * blocksize)));
end

spec_amp = abs (spec);
spec_phs = angle (spec);

smean = sqrt (sum (spec_amp (:, 1 : 200) .^2, 2) / 200);
spec_amp_blur = conv2 (spec_amp, blur_map, 'same');
smean_blur = conv2 (smean, sum (blur_map, 2), 'same');

addr_signal_a = spec_amp_blur > val_a * repmat (smean_blur, 1, bnum);
addr_signal_a_ex = conv2 (1 * addr_signal_a, ones (height, width), 'same');
addr_signal_b_val = zeros (size (spec_amp));
addr_signal_b_val (addr_signal_a_ex == 0) = val_a;
addr_signal_b_val (addr_signal_a_ex ~= 0) = val_b;
addr_signal_b = spec_amp_blur > repmat (smean_blur, 1, bnum) .* addr_signal_b_val;
addr_signal_b = 0 ~= conv2 (1 * addr_signal_b, blur_map, 'same');

spec_amp_ed = spec_amp;
spec_amp_ed (~addr_signal_b) = spec_amp_ed (~addr_signal_b) / fade;
spec_ed = spec_amp_ed .* exp (i * spec_phs);

wdata_ed = zeros (dsize, 1);
for k = 1 : bnum
    addr = (k - 1) * blocksize + (1 : 4 * blocksize);
    wdata_ed (addr) = wdata_ed (addr) + env .* real (ifft (spec_ed (:, k)));
end
wdata_ed = wdata_ed / 1.5;
wavwrite (wdata_ed, fs, nbits, [filename, '_ed_3_vals_2.8_1.5.wav'])
 

fragment_view_2.m

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

clc
clearvars
format compact

filename = 'fragment.wav';
blocksize = 512;
val_a = 2.8;
val_b = 1.5;
blur_map = ones (5, 3) / 15;
width = 7;
height = 91;

[wdata, fs, nbits] = wavread (filename);
dsize = size (wdata, 1);
bnum = floor (dsize / blocksize) - 4;
env = 0.5 * (1 - cos (2 * pi * ((1 : 4 * blocksize)' - 0.5) / (4 * blocksize)));

spec = zeros (4 * blocksize, bnum);
for k = 1 : bnum
    spec (:, k) = fft (env .* wdata ((k - 1) * blocksize + (1 : 4 * blocksize)));
end

spec_amp = abs (spec (1 : 2 * blocksize, :));
%smean = sqrt (sum (spec_amp .^2, 2) / bnum);
smean = sqrt (sum (spec_amp (:, 1 : 200) .^2, 2) / 200);

spec_amp_blur = conv2 (spec_amp, blur_map, 'same');
smean_blur = conv2 (smean, sum (blur_map, 2), 'same');

addr_signal_a = spec_amp_blur > val_a * repmat (smean_blur, 1, bnum);
addr_signal_a_ex = conv2 (1 * addr_signal_a, ones (height, width), 'same');
addr_signal_b_val = zeros (size (spec_amp));
addr_signal_b_val (addr_signal_a_ex == 0) = val_a;
addr_signal_b_val (addr_signal_a_ex ~= 0) = val_b;
addr_signal_b = spec_amp_blur > repmat (smean_blur, 1, bnum) .* addr_signal_b_val;
addr_signal_b = 0 ~= conv2 (1 * addr_signal_b, blur_map, 'same');
spec_noise = spec_amp_blur < 2.5 * repmat (smean_blur, 1, bnum);
smean_noise = sqrt (sum (spec_amp .^2 .* spec_noise, 2) ./ sum (spec_noise, 2));
%spec_noise = spec_amp_blur < 2 * repmat (smean_noise, 1, bnum);
spec_sgnl = 1 - spec_noise;

imdata = log10 (spec_amp) * 0.5 + 0.3;
%imdata = log10 (spec_amp_blur) * 0.5 + 0.3;
imdata = max (0, imdata);
imdata = min (1, imdata);
%imdata = 1 - imdata;

subplot (1, 5, 1 : 4)
image (flipdim (cat (3, imdata .* addr_signal_b, imdata, imdata .* (1 - addr_signal_b + addr_signal_a)), 1))
%image (flipdim (cat (3, imdata, imdata, imdata), 1))
axis off

subplot (1, 5, 5)
plot (zeros (1, 2 * blocksize), 1 : 2 * blocksize, 'k')
hold on
plot (smean, 1 : 2 * blocksize)
plot (smean_noise, 1 : 2 * blocksize, 'r')
hold off
xlim ([0, max(smean)])
ylim ([0, 2 * blocksize])
axis off
 

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

fragment.wav
fragment_ed_3_vals_2.8_1.5.wav
fragment_ed_3_w7_vals_3.5_2.1.wav

Sony SoundForge Noise Reduction PlugIn, конечно, по-удобней будет, благо что профессионалы клепали для профессионалов, но я тоже старался как мог...

[Schrodinger's cat]

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

я тут распинался, формулы выводил, рассказывал, интересно, для кого я это делал?

Так:
- Имеем вещественный сигнал. $20 \sin(\frac{2 \pi}{30}+0.5)$ Так же 200 отсчетов.
- Создал Ганкелеву матрицу из сигнала.
- Сделал ее SVD разложение. $USV^T$
- Получилось 2 положительных сингулярных числа.
- Взял эти 2 вектора соответствующие этим двум числам и восстановил аппроксимацию. Вроде совпадает с исходным.

Это шаги чисто АГК. Ну только без центрирования, нормирования и создания корреляционной матрицы.

Несколько раз перечитал что вы писали дальше, но не понимаю.
Дальше то что делать? Тут ничего комплексного не было а у вас потом все в комплексных терминах.

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

А вы пробовали

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

[ilghiz]

Schrodinger's cat как вы выше написали, вы на кошках тренируетесь, так возьмите комплексный сигнал и сделайте, что я написал и получите комплексную экспоненту вашего сигнала.

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

[Schrodinger's cat]

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

Пусть у вас сигнал - $f_i$ одна комплексная экспонента

Хорошо, заданную гармонику $20 \sin(\frac{2 \pi}{30}+0.5)$ можно представить в виде $20 \frac{e^{i\frac{2\pi}{30}+0.5}-e^{i\frac{2\pi}{30}+0.5}}{2}$.
Для простоты примем $20 e^{i\frac{2\pi}{30}+0.5}$. Допустим задана такая комплексная экспонента.
Коэффициентов мы пока не знаем: $ae^{i\varphi }$

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

$a f_i = f_{i+1}$

$f_i = ae^{i\varphi }$ так?
В $a f_i$ коэффициент $a$ как понимать, это амплитуда комплексной экспоненты или это коэффициент уравнения а амплитуда подразумевается спрятана в $f_i$?

Цитата:

это ваша задача

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

[svv]

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

заданную гармонику $20 \sin(\frac{2 \pi}{30}+0.5)$ можно представить в виде $20 \frac{e^{i\frac{2\pi}{30}+0.5}-e^{i\frac{2\pi}{30}+0.5}}{2}$.

Тут надо кое-что исправить:
Гармонику $20 \sin(\frac{2 \pi}{30}{\color{magenta}n}+0.5)$ можно представить в виде $20 \dfrac{e^{i\frac{2\pi}{30}{\color{magenta}n}+0.5}-e^{{\color{magenta}-}i\frac{2\pi}{30}{\color{magenta}n}{\color{magenta}-}0.5}}{2{\color{magenta}i}}$.
Маленький и незаметный множитель $i$ в знаменателе меняет фазу на $\frac{\pi}2$.

[Schrodinger's cat]

svv Спаcибо.

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

Есть еще один метод, посдвигать все на две размерности, получить тензор, и его тензорным разложением покромсать...

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

[svv]

Дальше Вам надо каждое слагаемое представить в виде $ae^{\lambda n}$, чтобы не пропала явная зависимость отсчётов $x_n$ от «дискретного времени» $n$:
$x_n=20 \dfrac{e^{i\frac{2\pi}{30}n+0.5}-e^{-i\frac{2\pi}{30}n-0.5}}{2i}=a_1e^{\lambda_1 n}+a_2e^{\lambda_2 n},$
где
$\begin{array}{ll}a_1=-10ie^{+0.5}&\lambda_1=+i\frac{2\pi}{30}\\[0.8ex] a_2=+10ie^{-0.5}&\lambda_2=-i\frac{2\pi}{30}\end{array}$

Представили? Теперь немного "теоретической практики".

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

2) Пусть $x_n=a_1e^{\lambda_1 n}+a_2e^{\lambda_2 n}$ (как в Вашем примере).
Проверьте, что тогда для любого $n$ выполняется
$x_{n+2}-(e^{\lambda_1}+e^{\lambda_2}) x_{n+1}+e^{\lambda_1}e^{\lambda_2} x_n=0$

3) Пусть $x_n=a_1e^{\lambda_1 n}+a_2e^{\lambda_2 n}+a_3e^{\lambda_3 n}$.
Проверьте, что тогда для любого $n$ выполняется
$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$

Таким образом, выполняются линейные (относительно отсчётов $x_n$) соотношения, причём коэффициенты при отсчётах постоянны (не зависят от $n$), не зависят от амплитуд $a_k$, а зависят только от $\lambda_k$.

На этом и основан метод Прони. А именно, будем действовать в обратную сторону: сначала найдём, какое линейное соотношение с постоянными коэффициентами выполняется для соседних отсчётов (т.е. найдём эти коэффициенты при отсчётах). А потом по коэффициентам восстановим $\lambda_k$.