Как почистить сигнал от гауссовского шума?

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

Евгений Машеров
Я оказывается способ 3 предлагал, не зная что он называется когерентным накоплением.

[sup]

Мне кажется, что оптимальный как раз №3.
Пусть у нас имеется сигнал $s_k$ ($k = \overline {1,n}$)
$s_k = Af_k +Nn_k$,
где $f_k$ известна, а $\xi_k$ - случайная величина, распределенная по Гауссу с нулевым средним и единичной дисперсией. Амплитуды $A, N$ неизвестны. Задача - по $s_k$ найти амплитуду $A$.
Ясно, что оптимальным будет назначить $A' = \frac {\sum f_k\xi_k}{\sum f_k^2}$
Аналогично этому можно рассмотреть и случай нескольких сигналов. Пусть, например,
$s_k = Af_k + Bg_k + Nn_k$.
Для простоты считаем их ортогональными - $\sum f_kg_k = 0$.
Тогда
$A' = \frac {\sum f_k\xi_k}{\sum f_k^2}$

$B' = \frac {\sum g_k\xi_k}{\sum g_k^2}$

С проблемой "дрейфа частоты" можно бороться следующим образом - искать максимум отклика $A'(\nu)$ по $\nu$, где $\nu$ - область частот около 50 гц.

-- Чт авг 07, 2014 13:03:39 --

Так, в формулax косяк, надо
$A' = \frac {\sum f_ks_k}{\sum f_k^2}$

$B' = \frac {\sum g_ks_k}{\sum g_k^2}$

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

Ну, таким способом из ЭЭГ вытягивают реакцию на стимул (метод вызванных потенциалов). Правда, проблема неточного знания выделяемой частоты отсутствует, поскольку момент подачи стимула фиксируется. Здесь, видимо, стоит брать число периодов усреднения такое, чтобы фаза сигнала сильно не ушла.
Кстати, можно считать не только среднее, но и дисперсию относительно среднего (что пригодится, чтобы понять, действительно ли это сигнал или артефакт).
Тут может оказаться полезным тождество $(x-E(x))^2=E(x^2)-(E(x))^2$
То есть накапливать отдельно сигнал, отдельно его квадрат, а потом оценивать дисперсию относительно усреднения.

[Boss03]

oleg777,
Matlab не пробовали привлекать? Литературы по нему в части удаления шумов предостаточно.

[sup]

Мне приходилось сталкиваться с похожей задачей. Правда в некотором смысле "наоборот". Надо было "чистить" ЭКГ от наводок 50 гц. Там проблемы есть/нет/артефакт не было Все было очень хорошо видно.
Что касается задачи поиска сигнала под шумом, то ее постановка была именно такая, как я сказал. Сигнал есть, но амплитуда не известна (может и равна 0). Дан сигнал, требуется выдать число.
Для гауссова шума ничего лучше чем я указывал предложить скорее всего нельзя (наверняка это уже давным-давно известно спецам). Это очевидно для поиска образца с носителем в одной единственной точке. А общий случай к этому сводится заменой переменных. Куда более интересным оказался этот же вопрос об амплитуде в случае равномерно распределенного шума. Там тоже можно дать хороший способ предсказания. Впрочем, как я и говорил, это все наверняка известно.

[oleg777]

Boss03 в сообщении #893907 писал(а):

oleg777,
Matlab не пробовали привлекать? Литературы по нему в части удаления шумов предостаточно.

Нет, потому что некогда. Все равно мне надо свой код писать.

Может, это и неправильно...

-- Чт авг 07, 2014 16:32:37 --

Евгений Машеров в сообщении #893848 писал(а):

(Оффтоп)

"Вызывает интерес и такой ещё разрез"

Если суммарная мощность белого шума в 81 ($9^2$) раз выше синусоиды - то при разложении по 2048 точкам спектра пик, соответствующий синусоиде, будет в 25 раз выше шумовых. И даже если заложиться на то, что частота не точно 50Гц (ну, или частота оцифровки немного дрейфует) и пик попал между "зубьев", что также привело к появлению боковых лепестков, и пики 40% от максимального, а среди шумовых есть разной, случайно меняющейся амплитуды, и максимальный втрое больше среднего - всё равно получается отчётливый пик.

Разобрался, это у меня шум был не очень случайный. Поправил это дело.

-- Чт авг 07, 2014 16:46:31 --

--------------

Беру тайм-аут для чтения ваших советов.

И все-таки очень хотелось бы ссылку на хорошие книжки по этому поводу, и по первому способу и, особенно, по второму. Чтобы основательно подковаться.

[oleg777]

Я проанализировал ситуацию и понял: мой фильтр идеологически близок к варианту №2. Отличается он тем, что вместо искусственного добавления гармоник он превращает шум в гармоники (с известным спектром).

Теперь мне есть чем заняться: сравнивать второй метод с моим.

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

Как-то задавался вопросом, есть ли разница между двумя подходами избавления от помех? Период полезного сигнала $Т$ известен и имеется последовательный набор оцифрованных значений длинной $nT$.
1. Каждый $i$-тый период раскладываем в ряд Фурье и определяем $1-$ую гармонику.
Далее находим среднее значение.
2. $nT$-считаем периодом и разложением находим $n$-ю гармонику.
Результаты должны совпасть, но проверить это как-то не получилось.

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

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

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

Я понимаю что окна нужны когда не известна частота гармоники.

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

Окна нужны для уменьшения влияния разрывов на краях.
Если можно согласовать частоту опроса с частотой интересующего сигнала - можно от разрывов избавиться.
Но именно в данном случае я бы просто усреднял. Без Фурье.

[profrotter]

sup дело пишет: тут задача оценки параметра сигнала. Про оценку амплитуды сигнала в белом гауссовом шуме можно найти, например, тут http://strts-online.narod.ru/files/lec9.pdf в п.9.6.1. Только следует учесть один важный ньюанс: в этой теме, насколько я понял, следует рассматривать в качестве полезного сигнал со случайной начальной фазой и отношение правдоподобия брать для него. После чего, как я предполагаю навскидку, вылезет коварная нелинейнщина связанная с необходимостью обращения функции Бесселя нулевого порядка.

[sup]

Мне кажется, что проблем с фазой нет.
Имеется весьма простой подход к этой задаче. Удобно перейти к формулировке в $R^n$. Тогда задача выглядит так
$\bar s = A \bar f + N \bar n$,
где
$\bar s$ - зашумленный сигнал
$\bar f$ - образец
$\bar n$ - вектор с независимыми компонентами, случайно распределенными по Гауссу.
Совместное распределение этих компонент зависит только от "радиуса", а посему инвариантно относительно поворотов. Значит если $Q$ - ортогональная матрица, то
$\tilde s = Q\bar s = A Q\bar f + N \tilde n$
Подбираем $Q$ так, чтобы $\bar g = Q\bar f$ имело только одну ненулевую компоненту, а именно $g_1 \neq 0$. Теперь уже ясно, что у вектора $\tilde s$ информативна только первая компонента. Оставшиеся компоненты содержат только независимый гауссов шум. Их величина никак не может нам помочь в определении $A$. Таким образом задача свелась к анализу одного единственного равенства
$\tilde s_1 = A g_1 + N \tilde n_1$
Ну а теперь уже ясно, что следует положить $A = \frac {\tilde s_1}{g_1}$, поскольку матожидание шума равно 0.
Разумеется, вся матрица $Q$ нам не нужна. Всего лишь одна строка. С точностью до множителя это в точности вектор $\bar f$. Отсюда и возникает формула
$A'= \frac {\sum f_k s_k}{\sum f_k^2}$
А что делать, когда образец состоит из двух компонент, и, следовательно, требуется определить две амплитуды?
Да то же самое. Пусть у нас образцы $\bar f$, $\bar h$. Без потери общности можно считать, что $\sum f_k h_k =0$. В противном случае можно рассмотреть $\bar h_1 = \bar h - \lambda\bar f$. Итак, пусть образцы "ортогональны". Тогда опять таки рассматривая подходящий поворот, можно получить два независимых соотношения на искомые амплитуды. Как следствие получаем формулы
$A'= \frac {\sum f_k s_k}{\sum f_k^2}$

$B'= \frac {\sum h_k s_k}{\sum h_k^2}$
Теперь уже можно решить вопрос с фазой. Достаточно рассмотреть два образца $\sin k\nu$ и $\cos k\nu$. Находим их амплитуды и дело в шляпе.

(Оффтоп)

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

[levtsn]

если БПФ то фильтр не нужен. бпф уже дает ответ.

[oleg777]

Мне вот подсказали еще такую новую штуку: метод эмпирической модовой декомпозиции. По этому поводу кто-нибудь хочет высказаться?

P.S. Пардон, что никак не реагирую на ваши посты. Это потому, что я врубаюсь очень медленно, заодно теорию вспоминаю.

-- Сб авг 09, 2014 22:29:02 --

Мне вот подсказали еще такую новую штуку: метод эмпирической модовой декомпозиции. По этому поводу кто-нибудь хочет высказаться?

P.S. Пардон, что никак не реагирую на ваши посты. Это потому, что я врубаюсь очень медленно, заодно теорию вспоминаю.