Прони на данных с локальным шумом или с дырками

[ilghiz]

Добрый день,

навеяно обсуждением метода Прони в соседней теме https://dxdy.ru/topic150785.html.

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

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

Вопросов два:

1. как устойчиво реконструировать сигнал в области, где есть шум?
2. как устойчиво детектировать регион, где есть шум.

В качестве примера - берем синус на какой-то частоте, набрасываем на него немного шума, и дискретизуем пусть на 1000 точек, так, чтобы частота была такой, что в эти 1000 точек уместилось, скажем 100 периодов. А теперь с 450 по 550 отсчета набросим поверх шума.

Цель этого - уменьшить уровень отношения шим/сигнал в таких данных в спектральной области.

[Schrodinger's cat]

В такой постановке вопроса кажется выгодным вариант метода с SVD который по идее должен отделить шум.
Топ n главных компонент - сигнал, остаток - шум.
А там уже можно и сам район по остатку определить.

Марпл так же про это пишет, глава 11.9 "Идентификация экспонент в шуме."

[ilghiz]

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

В такой постановке вопроса кажется выгодным вариант метода с SVD который по идее должен отделить шум.

Не, не отделит, и сильно завязнет, проверено.

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



Если в центральном регионе нет этой помехи, то все идеально работает, и по Прони, и по Фурье.

[Cos(x-pi/2)]

(по Фурье)

Помеху, мало перекрывающуюся по спектру с полезным сигналом, удаётся устранить КИХ-фильтром с П-образной АЧХ; такой фильтр полезный сигнал пропускает, а помеху нет:

графики сигналов во времени

спектр исходного сигнала

АЧХ фильтра

[ilghiz]

Спасибо, Cos(x-pi/2) за ответ!

Не, через КИХ или через Фурье я умею, но мне там не все нравится и у меня есть гипотеза, что через Прони можно получить более стабильные результаты.

[B@R5uk]

Применять Прони ко всему участку с приведёнными сигналами не корректно, на мой взгляд. Метод предполагает, что сигнал состоит из синусоид с фиксированной амплитудой или с амплитудой, меняющейся по экспоненте. Тут же налицо скачкообразное появление и исчезновение второй синусоиды. Ряд синусов с экспоненциальной огибающей (или даже просто синусов как в Фурье), конечно, сможет аппроксимировать П-образную огибающую, но смысл? Затраты не окупают результат, да с Прони результат будет так себе.

Тут, очевидно, логичнее обработать участок с двумя синусоидами отдельно от всего остального. А участки справа и слева — отдельно от особого участка (можно даже оба одним махом, если известно, что синусоида за время помехи не съехала по какой-то причине по фазе, такое тоже может быть).

-- 12.10.2022, 12:49 --

Cos(x-pi/2), а как вы инициализировали свой фильтр, что он вам from the get-go синус на нужной частоте стал выдавать без переходного процесса?

[ilghiz]

Спасибо B@R5uk за комментарий!

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

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

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

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


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

[B@R5uk]

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

Я хочу использовать Прони метод для распознавания участка, где есть помеха.

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

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

[ilghiz]

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

Для того, чтобы распознать помеху (или что-либо ещё), модель должна поддерживать в себе такую особенность. Метод Прони таким не является, следовательно, для распознавания не подходит.

Не соглашусь. Классический способ распознавания таких помех - считаем норму сигнала и норму производной сигнала и усредняем в окне или как-то локально. Если вдруг произошел скачек по таким нормам в каком-то регионе - то в этом регионе имеется помеха. Очевидно, что норма и производная - как раз автоматически в Прони считаются, то есть можно переиспользовать.

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

Фурье, кстати, тоже не очень для этих целей (для такого рода помех).

Не, Фурье как раз для удаления помехи используют - берут регион с помехой, зануляют, делают Фурье, далее находят в спектральной области 1-2-3 самых ярких пика, все остальное удаляют и считают обратное Фурье и им заполняют регион с помехой. Так как два Фурье (вперед и назад) все-таки довольно трудоемкая операция, хочется заменить это на Прони. Вернее оно заменяется вроде бы, но так как тут много про Прони обсуждалось, хочется вместе об этом подумать, возможно придумается что-то еще лучше.

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

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

Не, мне кажется тут надо формулировать основную цель как - улучшение сигнал/шума у такого рода сигналов в спектральной части.

ЗЫ: это не моя задача, я просто предлагаю на эту тему пофилосовствовать и, может, что-то интересное придумать.

[Cos(x-pi/2)]

(B@R5uk)

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

как вы инициализировали свой фильтр, что он вам from the get-go синус на нужной частоте стал выдавать без переходного процесса?

В предположении, что на выходе фильтра сигнал запаздывает примерно на N/2 отсчётов, где N - длина фильтра, отбрасываю первые N/2 отсчётов в профильтрованном сигнале. Т.е. перенумеровываю - сдвигаю результат влево по оси времени на N/2 отсчетов (точнее: здесь сдвинуто на (N+1)/2, так как в приведённом примере N нечётное).

Вот иллюстрация графиками. На верхнем рисунке Sout это профильтрованный сигнал без устранения задержки. На нижнем красный сплошной график Sf - профильтрованный сигнал с устранённой задержкой, график чёрными точками S1- исходный синус до прибавления помехи: с переходным процессом и без.

[Schrodinger's cat]

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

[B@R5uk]

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

берут регион с помехой, зануляют...

Вот вы говорите "берут регион". Эта операция никак к Фурье не относится, информация о расположении региона получается отдельными методами. Такая же проблема и с Прони.

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

далее находят в спектральной области 1-2-3 самых ярких пика, все остальное удаляют

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

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

предлагаю на эту тему пофилосовствовать и, может, что-то интересное придумать.

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

-- 12.10.2022, 17:41 --

Cos(x-pi/2), спасибо за разъяснение. Что-то про задержку фильтра как-то и забыл совсем. Вот только у вас на графике осцилляции фильтрованного сигнала почему-то начинаются как-то совсем не сразу и довольно внезапно (фильтр "раскачивается" целых 3 периода сигнала практически без отклика). С чем это связано? С крутизной характеристики в области перехода от пропускания к заграждению?

[Cos(x-pi/2)]

(Оффтоп)

B@R5uk

Просто проявился тот факт, что при 4095-точечном фильтре задержка примерно равна 4095/2, т.е. 2048 отсчётов. Это много: почти пятая доля всей длины сигнала в данном примере. А количество периодов синусоиды в задержке, т.е. на этих 2048 отсчётах, зависит от частоты синусоиды и могло бы быть другим при другой её частоте.

У меня на графиках по горизонтальной оси отложено время в секундах, частоту дискретизации я выбрал равной 16 Гц, чтобы на 600 секундах поместилось много отсчётов - почти 10 тыс. (9600 отсчётов, если точно говорить). Частоту основной синусоиды я выбрал так, чтобы на 600 секундах поместилось 14 её периодов, - чтобы её график был бы похож на график синусоиды у ilghiz. Вот и вышло, что на почти впятеро меньшем числе отсчётов (2048) помещается примерно 14/5 периодов такой синусоиды.

Длина фильтра 4095 у меня в программке была установлена "по умолчанию", т.е. она такая оказалась "от балды". Программка допускает длину $2^{\text{целое}}-1.$ Сейчас посчитал с более короткими фильтрами. При 511 уже видна некрасивая корявость на месте помехи; при 255 и сама помеха недофильтрованная остаётся. А при 1023 фильтруется примерно так же, как и более длинными фильтрами. Задержка при этом более короткая: 512 отсчётов; вот картинки, аналогичные предыдущим: результат с фильтром N=1023

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

Поэтому, чтобы мне дальше не оффтопить в этой теме, думаю, лучше обсуждение учебных примеров фурье-спектров, спектрограмм, конкретики алгоритмов, если, конечно, это будет интересно Schrodinger's cat, не вызовет у него возражений, и если, конечно, появятся ко мне вопросы упражненческого уровня (не сильно сложные, у меня опыт невелик - всего лишь типа радиолюбительского :-), вести в его теме про "Разложение Фурье".

[Schrodinger's cat]

Cos(x-pi/2) в сообщении #1566586 писал(а):

это будет интересно Schrodinger's cat, не вызовет у него возражений

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

По теме вспоминается детектор аномалий.
https://keras.io/examples/timeseries/timeseries_anomaly_detection/
Еще автоэнкодеры разного толка, есть подтип как раз шумодавов.
Но это слишком для численных методов, я просто упомяну, что это есть.

Можно спектр усреднять за n районов, и его сравнивать с текущим. С не зашумленным разница будет меньше чем с зашумленным.

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

Но мне все равно кажется разумной идея применить, что то из PCA, SVD... если взять достаточный кусок до и после "вырезаемого аппендицита" то он должен восстановить. В чем там проявляется неустойчивость?
ilghiz раз Вы это проверяли, поясните пожалуйста в двух словах.

С Прони конечно не связано, но пока придумал только это.

[ilghiz]

Огромное спасибо всем за интересные обсуждения!

Позвольте на пальцах рассказать что получается. Для этого я приложу спектр уже в логарифмической шкале.



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

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

Конкретно тут это делалось так:

1. зануляли шумовой кусок, делали Фурье, искали самую яркую гармонику, все остальное выбрасывали, делали Фурье назад и втыкали это вместо помехи и делали снова Фурье.
2. Альтернативно (картинка получается один-в-один) по Прони вытаскивали частоту, на ее основе наименьшими квадратами искали фазу, восстанавливали это в тот сегмент, и снова делали Фурье.

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

Хочется понять на сколько тут все стабильно, если сам сигнал будет состоять из многих компонент. Прони-то захочет много компонент и проиграет Фурье по скорости. А первый метод тоже может завязнуть в поиске того, что надо выбрасывать, а что - нет. И вот тут есть подозрение, что Прони может одной-двумя компонентами схватить больше информации, чем Фурье.

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