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

[Schrodinger's cat]

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

$$Y_n=A_1\sin(\frac{2 \pi n}{F_1}+0.7)+A_2\sin(\frac{2 \pi n}{F_2}+0.3)$$

Где:
$A_1$ и $A_2$ - амплитуды гармоник 20 и 10 соответственно;
$F_1$ и $F_2$ - периоды гармоник, 50 и 6 соответственно;
$n$ - номер дискретного отсчета;
$0.7$ и $0.3$ - начальные фазы гармоник сигнала.
$N$ - всего точек, 198;
$n$ - номер дискретного отсчета $= 0..N-1$

(Оффтоп)

Делаю разложение Фурье. Делаю руками, но не суть, оно совпадает с MathCad'овским.
Строю АЧХ, все хорошо.
По АЧХ нахожу первую гармонику на отсчете 4, вторую на отсчете 33.
Спектр я естественно складываю пополам относительно середины и суммирую. Откуда получаются амплитуды 20 и 10. Отлично!

(Оффтоп)

загрузка фото

Но вот с ФЧХ затык, делаю все тоже самое как с АЧХ получаю нули, оно и не мудрено, ФЧХ симметрична относительно середины.
По разному пробовал, собственно не могу получить начальные фазы гармоник.
Помогите пожалуйста в чем может быть проблема?

Определил формулу гармоники.
Так как этой же формулой можно экстраполировать назовем ее Predict:

$$Predict(a,k,f,x)=a\sin(\frac{2 \pi x}{k}+f)$$

Где:
$a$ - амплитуда;
$k$ - период;
$f$ - начальная фаза;
$a \, k \, f$ - найденные значения.
$x$ - координата, которую можно продлить за диапазон и получить прогноз.

Найденные значения подставляю в формулу. Вид гармоник совпадает с исходником (амплитуды и частоты найдены верно) за исключением фаз $f$.

Что еще будет нужно для ясности вопроса допишу, выложу. Скажите что.

ПС: Я не студент, ничего бежать сдавать не нужно.

[Pphantom]

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

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

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

И заодно еще поправьте текст. "Окружные" частоты обычно называются "циклическими", а те величины, которые вы так называете, вообще являются периодом. Что такое "Pred", надо как минимум пояснить.

[Pphantom]

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

[Schrodinger's cat]

Разложение Фурье дает комплексные числа, чтобы получить амплитуды делаю так:

$$A_{n}=\sqrt{\operatorname{Re}_{n}^{2}+\operatorname{Im}_{n}^{2}}$$

Где:
$n$ - номер комплексного числа в разложении;
$\operatorname{Re}$ и $\operatorname{Im}$ - соответственно действительная и мнимая часть.
Но так как АЧХ получается зеркальной относительно середины, делаю еще так:

$$Amplitude(n)=A_{n}+A_{N-n}$$

Где:
$N$ - всего точек в исследуемом ряде.
Отлично, получаются искомые амплитуды 20 и 10.

Фазы получаю соответственно так:

$$P_{n}=\arctg2(\operatorname{Re}_{n}/\operatorname{Im}_{n})$$

Где:
$\arctg2$ - вариант функции арктангенса учитывающая знаки.
Но формула по аналогии с амплитудами дает нули, из за симметричности диаграммы:

$$Phase(n)=P_{n}+P_{N-n}$$

Но и в $P_{n}$ и $P_{N-n}$ я не нахожу искомых значений 0.7 и 0.3!

Как правильно нужно находить начальные фазы гармоник из разложения?

[ozheredov]

Schrodinger's cat, почему после нахождения периодов не поискать аппроксимацию в линейном виде:

$u_n = C + A_1 \cos (\frac {2 \pi}{T_1}) + B_1 \sin (\frac {2 \pi}{T_1}) + A_2 \cos (\frac {2 \pi}{T_2}) + B_2 \sin (\frac {2 \pi}{T_2})$

[Schrodinger's cat]

ozheredov сперва вообще не понял о чем речь, но потом с трудом догадался. Спасибо, для прогноза это работает.

Но как отсюда найти начальные фазы?

[ozheredov]

Schrodinger's cat
Арктангенс отношения коэффициента разложения к норме вектора коэффициентов разложения по гармоникам одной частоты -- у Вас есть похожая формула после фразы про фазы )

[Schrodinger's cat]

ozheredov вы про эту формулу

$$P_{n}=\arctg2(\operatorname{Re}_{n}/\operatorname{Im}_{n})$$

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

$$P_{n}=\arctg2(\operatorname{?}_{n}/\operatorname{Amplitude}_{n})$$

[ozheredov]

Schrodinger's cat, я дебил -- меня переглючило Вот так надо:

$Phase_1 = \arctg {\frac {A_1}{B_1}}$
$Phase_2 = \arctg {\frac {A_2}{B_2}}$

[Schrodinger's cat]

ozheredov теоретически да, коэффициенты $A_1$ и $B_1$ находятся разложением, фактически это $\operatorname{Re}_1$ и $\operatorname{Im}_1$ но формула $Phase_1 = \arctg {\frac {A_1}{B_1}}$ (в программах нужно учитывать знак, поэтому используется arctg2) не дает искомого, в том то и дело.
При этом обратный Фурье дает правильный ряд. Следовательно коэффициенты найдены верно.

Провел эксперимент. Использовал разложение mathcad'а чтобы наверняка исключить ошибки в моих формулах.
Амплитуды получаются, фазы нет. В чем прикол...

[Александрович]

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

Амплитуды получаются, фазы нет.

Фаза "зашита" в $A_1$ и $B_1$. Откладываете косинус $\cos \varphi_1=\frac{A_1}{\sqrt{A_1^2+B_1^2}}$ и синус $\sin \varphi_1=\frac{B_1}{\sqrt{A_1^2+B_1^2}}$на единичном круге и находите угол $\varphi_1$.

[Schrodinger's cat]

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

$$Y_n=20 \sin(\frac{2 \pi n}{50}+0.7)$$

$N$ - всего точек сделаем 198 (просто считает чуть лучше, не суть те же 200).
$n=0..N-1$ - интервал ряда.
Делаю преобразование Фурье, встроенной функцией в mathcad, чтобы исключить возможные ошибки.
Я уже наперед знаю что искомая гармоника сидит на 4ом отсчете по частоте (по графику, да и собственно $200/50=4$ ).
$\operatorname{Re}_4=1.07 \cdot 10^3$
$\operatorname{Im}_4=-1.666 \cdot 10^3$
зеркальные значения с другого конца преобразования
$\operatorname{Re}_{N-4}=1.07 \cdot 10^3$
$\operatorname{Im}_{N-4}=1.666 \cdot 10^3$
Вычисляю амплитуду. Так как встроенное преобразование не нормирует коэффициенты сразу, придется разделить еще на N:

$$A_{n}=\frac{\sqrt{1070^{2}+(-1666)^{2}}}{198}=10$$

Тоже самое проделаем с зеркальными значениями $N-4$, так же получим 10 и сложим их в итоге получается искомая амплитуда 20.
Прекрасно, работает.
Попробуем считать фазы:

$$P_{n}=\arctg2(\frac{1070}{-1666})=-1$$

Тоже самое с зеркальными значениями $N-4$, получим 1.
Александрович по Вашим формулам получается 0,54 и -0,842.
С зеркальными значениями 0,54 и 0,842.
Арксинусы и арккосинусы этих значений дают 1 и -1.

Как видите ничего похожего на искомое 0,7 не получается!

[svv]

Schrodinger's cat

В этой области используются разные соглашения о знаках, нормировках и т.п. Я использую определения из справочника Корна по математике, стр.150, пункт 4.11-4(b) «Гармонический анализ периодических функций» и стр.687, пункт 20.6.6(а) «Численный гармонический анализ и тригонометрическая интерполяция».

В соответствии с этим стандартом (которому, думаю, следует и Mathcad), фаза $k$-й гармоники — это слагаемое $\varphi_k$ под косинусом в правой части:
$f(t)=\frac 1 2 a_0+\sum\limits_{k=1}^\infty(a_k\cos\omega_k t+b_k\sin\omega_k t)=A_0+\sum\limits_{k=1}^\infty A_k\cos(\omega_k t+\varphi_k),$
где
$f(t)$ — периодическая функция с периодом $T$
$\begin{array}{l}\omega_k = \frac{2\pi}T k \\ [0.8ex]A_0=\frac 1 2 a_0 \\ [0.8ex]A_k=\sqrt{a_k^2+b_k^2}, \; k=1,2,... \\ \varphi_k=\arg(a_k-ib_k)=-\operatorname{atan2}\;(b_k,a_k)\end{array}$
$\operatorname{atan2}$ — имя стандартной функции.

Значит, чтобы найти фазу, надо в Вашей функции выразить синус через косинус:
$Y_n=20 \sin(\frac{2 \pi n}{50}+0.7)=20 \cos(\frac{2 \pi n}{50}+0.7-\frac{\pi}{2})$
Так что фаза была бы $0.7-\frac{\pi}{2}\approx -0.870796$, если бы не один момент.

Дело в том, что, когда Вы изменили число точек $N=200$ на $N=198$, Ваш сигнал $Y_n$ в дискретном времени потерял свойство $Y_{N}=Y_0$:
$20\sin\left(\dfrac{2\pi\cdot 198}{50}+0.7\right)\neq 20\sin 0.7\right$
Теперь период сигнала не соответствует отрезку, на котором выполняется гармонический анализ. Это приводит, в частности, к тому, что ненулевой будет не только четвёртая гармоника, но и другие. Эффект аналогичен скачку функции при непрерывном изменении аргумента. Особенно чувствительны к этому фазы.

Вот первые несколько гармоник, сравните.
$N=200$
$\begin{tabular}{|l|r|r|r|r|} \hline $k$&$a_k$&$b_k$&$A_k$&$\varphi_k$\\ \hline 0& 0.00000000& 0.00000000& 0.00000000& \\ \hline 1& 0.00000000& 0.00000000& 0.00000000& \\ \hline 2& 0.00000000& 0.00000000& 0.00000000& \\ \hline 3& 0.00000000& 0.00000000& 0.00000000& \\ \hline 4& 12.88435374& 15.29684375& 20.00000000& $\color{magenta}-0.87079633$\\ \hline 5& 0.00000000& 0.00000000& 0.00000000& \\ \hline 6& 0.00000000& 0.00000000& 0.00000000& \\ \hline 7& 0.00000000& 0.00000000& 0.00000000& \\ \hline\end{tabular}$

$N=198$
$\begin{tabular}{|l|r|r|r|r|} \hline $k$&$a_k$&$b_k$&$A_k$&$\varphi_k$\\ \hline 0& -0.19738335& 0.00000000& 0.19738335& -3.14159265\\ \hline 1& -0.21229512& -0.09136635& 0.23112124& 2.73517879\\ \hline 2& -0.27234867& -0.22960327& 0.35621827& 2.44114787\\ \hline 3& -0.49227664& -0.60187869& 0.77755659& 2.25635707\\ \hline 4& 10.80406803& 16.83060897& 19.99993211& $\color{magenta}-1.00011535$\\ \hline 5& 0.38997342& 0.71831183& 0.81734397& -1.07341844\\ \hline 6& 0.19051384& 0.39495180& 0.43850022& -1.12135000\\ \hline 7& 0.12458706& 0.28068582& 0.30709358& -1.15305464\\ \hline\end{tabular}$

Вывод. Верните $N=200$ и замените синус на косинус в определении $Y_n$. Получится:
$\begin{tabular}{|l|r|r|r|r|} \hline $k$&$a_k$&$b_k$&$A_k$&$\varphi_k$\\ \hline 4& 15.29684375& -12.88435374& 20.00000000& $\color{magenta}0.70000000$\\ \hline\end{tabular}$

[Schrodinger's cat]

svv осмыслил все, что вы написали, большое спасибо за развернутый ответ.
Собственно эффект называется растекание спектра.
Поэкспериментировал. Работает.
Хочется это переложить в программу и максимально автоматизировать.

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

Попробовал задачу решать несколько иначе. Очевидно:
$A \cos(f)-Y_{0} = 0$
Где:
$A$ - амплитуда гармоники;
$f$ - начальная фаза гармоники;
$Y_0$ - начальное значение дискретного ряда гармоники;
Решая данное уравнение относительно $f$, можно найти начальную фазу.
Но очевидно тут те же проблемы, из за растекания спектра не получить точное значение амплитуды, что сильно влияет на результат. Да и вычислительно очень затратно, нужно выделить гармонику из спектра для получения начального отсчета.

Попробовал гетеродирование сигнала.
Умножаю сигнал на $\cos(dn)$
Где:
$d$ - малая величина, на которую мы хотим сместить гармонику по частоте в спектре АЧХ.
$n$ - номер дискретного отсчета.
Гетеродирование применяю в первый раз.
Да, сигнал разделяется на 2 рядом стоящих, можно подобрать такое значение $d$ чтобы один из них совпал с дискретным отсчетом спектра и не растекался. Но второй все равно стоит рядом и его растекание все портит.
По теории если $d$ выбрать достаточно большим, то второй образовавшийся сигнал должен уйти за пределы спектра. Но такого у меня не происходит, они вместе ездят по спектру на достаточно близком расстоянии.
В общем данным методом точно амплитуду я так и не смог выделить.

Кто нибудь знает как грамотно применить гетеродирование и таки сместить гармонику в спектре АЧХ для ее четкой идентификации?
Нет ли еще каких либо способов выделить параметры гармоник, желательно автоматизированных без ручного анализа?

[ilghiz]

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