Есть оцифрованный сигнал в виде постоянной составляющей с шумом и импульсными помехами (по оси х - время в секундах).
Нажмите для просмотра прикрепленного файла
Точно известен интервал (+- 5*10^-6), в котором может находиться величина постоянной составляющей (этот интервал показан на рисунке).
Надо определить величину постоянной составляющей. За минимальное время с начала обработки величина постоянной составляющей должна быть определена с точностью около около +- 1*10^-8, при последующей обработке точность должна возрастать (чем больше - тем лучше).

С линейной фильтрацей пока у меня что-то плохо получается. В ЦОС я не силен, поэтому пробовал каскадное включение 10 БИХ-фильтров 1-го порядка (Yn + Yn*(K-1) = K*Xn при K=1/2^6) - слишом много времени уходит (где-то 1500-2000 с) и импульсы мешают.
Пробовал вырезать участки с импульсами и потом "сшитый" сигнал пропускать через ФНЧ, но что-то не очень помогает.
Медианный фильтр не пробовал, т.к. не знаю какого размера окна взять, а при очень большом окне он постоянку исказит (наверное).

Посоветуйте как это лучше (и правильней) сделать.

К сожалению не все выбросы компенсируются противоположными, поэтому фильтрация фнч будет с медленными выбросами.
Единственное хорошее свойство Вашего сигнала это то, что выбросов в выборке немного.
Главное для каждой точки определить выброс это или нет и в случае выброса заменить ближайшим невыбросным значением.
Для этого достаточно знать приближенное значение измеряемой в настоящее время величины и дисперсию присутствующего в настоящее время шума.
Все значения за пределами трех сигм следует считать выбросами.
Таким образом предлагаемый алгоритм такой: Берем выборку 100-200 точек, по предыдущей грубой оценке измеряемой величины и дисперсии определяем новые уточненные значения. итд. Полученные грубые значения можно усреднить или сгладить скользящим средним.
Остается вопрос о первом измерении и о продолжении измерений после изменения измеряемого значения.
Для первого измерения можно использовать просто среднее выборки или среднее участка выборки размах значений которого не превышает Вашего диапазона, а вместо трех сигм Ваш диапазон. Смену значения можно определить по выкидыванию скажем 30% отсчетов.
Если математика поджимает сигмой можете назвать средневыпрямленное значение.

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

Ничего лучшего в голову не приходит.

http://electronix.ru/forum/index.php?showtopic=54351

Я делал бы так же, как написал GetSmart в последнем сообщении в той теме, ссылку на которую выдал psL - сложил бы все выборки и поделил на их количество. По-моему, лучшей фильтрации в природе не существует. Для простоты можно взять 4096 точек, потом и делить не придется.
Ну, или то, что там обсуждалось вначале. Это называется экспоненциальный фильтр. Ресурсов меньше потребуется. Чтобы выбросы меньше влияли на результат, их как-то выделять и отбрасывать.

Оцените максимальную плотность выбросов. Если на 1 выброс гарантированно приходятся два нормальных отсчёта, будет достаточно медианы с окном 3.

На картинке относительная амплитуда выбросов до 3.5*10^-4. При стремлении к 1*10^-8 один такой выброс 1-ым порядком придётся размазать по ~ 35000 отсчётам. Сомневаюсь, что медианный фильтр сможет нагадить больше, ну и он больше ФЧХ ломает.

SPACUM, спасибо за совет, попробую.

Это скользящее среднее. Пробовал в первую очередь. Результаты при приемлемом (т.е. 1024, например) размере окна очень плохие (а для меня больше 2000 уже совсем не приемлемо).
Кстати частный случай ДПФ в виде алгоритма Герцеля для n=0 именно в скользящее среднее и вырождается.

Это его экономисты и статисты так называют. Я же его обозвал БИХ-фильтр 1-го порядка (см. 1-ый пост) - он рекурсивный, значит БИХ.

P.S. Да! Забыл уточнить - частота дискретизации 1 Гц. Так что вариант скользящего среднего на 4096 - это больше часа. Надо где-то за 500-1000 отсчетов вычислить постоянку с приемлемой точностью (а вообще чем быстрее тем лучше). У меня только через 1500-2000 кое-как получается.


Нет. Есть места где 2 выброса (1 вверх, 1 вниз) идут друг за другом . Окна на 3 мало.
К тому же моя картинка - это пока частный случай. Есть предположение, что в других ситуациях шум помощнее будет, а выбросы будут чаще...
Мы же говорим о применении медианного фильтра ко всем отсчетам, так? При большем размере окна медианный фильтр постоянку может испортить (вроде бы...)

Основной вопрос - что Вы называете "истинной величиной постоянной составляющей". Если мы ничего не знаем априорно о том, что нам дается в виде отсчетов, то, естественно, это будет математическое ожидание - среднее без весов в данном случае.
Если же Вы знаете, что есть наводки (помехи), то их нужно убрать, но для этого нужно знать, что это такое. Соответственно и методы могут быть разными. Мне кажется, что наиболее разумно в данном случае исследовать функцию распределения и выбрасывать отсчеты, которые... (тут можно придумать несколько способов) искажают ее - отклоняют от Гаусса или Пуассона... Если имеется априорная уверенность. Можно, наверное, даже отбрасывать, используя "узость" распределения как целевую функцию. Смотреть можно на производную от сигмы по ширине оставленного коридора.
Примерно такое уже предлагали выше - отбрасывать все, что не влезает в 3 сигмы. Но тут нужно считать вероятность... При большом числе отсчетов вероятность попадания шума в 3 сигмы уже ненулевая, но какой будет вклад от таких отсчетов в среднее...

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

Р. Лайонс его называет "устройство экспоненциального усреднения". У него импульсная характеристика экспоненциальная. БИХ, конечно.

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

Суть задачи - имеется входной цифровой сигнал с частотой около 1 кГц, подверженный фазовой модуляции шумом. Частота исходного сигнала известна с точностью +- 5ppm. Необходимо восстановить исходную частоту сигнала. Начинать можно где-то с +-0,5 ppm.
Имеющиеся отсчеты - это выраженный в виде числа тактов опорного генератора (около 100 МГц) интервал между 1000 входными "зашумленными" импульсами.
Шум, модулирующий фазу, вообще говоря не очень стационарный, кроме того имеет медленно меняющуюся постоянную составляющую. Исключить или уменьшить влияние шума на исходный сигнал нельзя.

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

Вот.
Может что подскажите?



Работаю над этим

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

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

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

Безусловно, но наличие/отсутствие выброса будет видно невооружённым глазом.

Что ж тут загадочного - БПФ пересчитывается для каждого нового отсчёта.

Да-да, я уже понял вашу мысль. Все-таки, это как-то накладно.

Пересчитать 1 (один) столбец матрицы сложнее сигм, динамических порогов или цепочки из 10 фильтров? Ну тады ой...

Товарищи, ДПФ по полной не нужно. Мне нужен только отсчет спектральной плотности только в одной точке - на нулевой частоте. А этом случае ДПФ (а также в случае преоб. Фурье для непрерывного сигнала) получаем просто среднее значение. Усреднение. У меня усреднение приемлемого результата не дает. Пробовал окно на 1000 (и кажись еще на 2000) отсчетов. На выходе сигнал весьма грязненький.

Суть задачи уже озвучил выше.

Может в самой методе чего-то не так?

Этот же алгоритм можно реализовать немного по-другому.
Складываются не абсолютные значения выборок, а их отклонения от среднего значения. Среднее, в данном случае - это постоянное значение для всего процесса суммирования до некоторого момента времени. Стартовое значение среднего можно выбрать произольно. Затем оно будет скорректировано в процессе.
Поскольку суммируются значения со знаком, то текущее значение суммы будет колебаться вблизи нулевого значения. Количество выборок при таком суммировании не нужно ограничивать, а сами значения выборок не нужно хранить. Пока значение суммы по модулю не выходит за заданные нами пределы, можно считать, что выбранное среднее значение соответствует реальному среднему с заданной точностью. Как только сумма превысит заданный предел (по модулю) - включается счетчик сэмплов. Если в течение заданного времени (заданного количества сэмплов) отклонение не вернулась в прежний диапазон - корректируем текущее среднее значение, обнуляем текущую сумму. И продолжаем суммирование, используя уже новое среднее...

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

Может быть. Я, н-р, не понимаю, как при частоте выборок 1 Гц можно восстановить частоту сигнала 1 кГц. Или Вам достаточно "среднего по больнице" за тысячу периодов?
Если отвлечься от моих непоняток, то задача распадается на исключение из выборки отсчётов с неправильной динамикой (выбросов) и на расчёт собственно среднего значения выборки.
Если делать скользящее БПФ (ДПФ), то попадание выброса в его окно приведёт к появлению бинов в правой части спектра - имхо, вполне наглядный признак выброса. Что делать с таким отсчётом, думайте сами, но, вероятно, его стОит исключить из расчёта среднего значения.
Если Вас не устраивает "просто среднее значение" на выходе БПФ (ДПФ), придётся придумывать какую-нибудь более сложную модель, которая учитывает, н-р, предысторию или температуру, напряжение, фазу Луны... , т.к. для одного набора данных "просто среднее значение" будет совпадать со средним значением, посчитанным более другим образом.

По-моему, лучший результат даст, все-таки, КИХ, только не скользящего среднего, а с нормально подобранным окном. Это снизит влияние и импульсных помех. Окно должно обеспечить затухание 80-90 дБ (это в случае, если предварительно из сигнала вычесть среднее для сброса разрядности до 16 бит, а то 160 дБ это многовато.
***
Можно несколько раз (5-6) повторить КИХ скользящего среднего, если это менее затратно, чтобы опустить уровень затухания до 80-90 дБ.

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

Вернемся к медианному фильтру. Если из 3-х текущих отсчетов выбирать средний по величине, то выброс вверх, а за ним сразу вниз оба будут отброшены. Если выбросов больше, сделать фильтр длиннее, выбирать средний из 5-ти отсчетов. Мало? - из 7-ми...

Меряется время между соседними фронтами. Результаты 1000 измерений суммируются и выдаются в виде очередного "сэмпла".
Получается каждый сэмпл - среднее значение мгновенной частоты на интервале 1000 периодов.
Соответственно частота сэмплов ~ 1 Гц.
Если надо, объясню почему взял 1000, а не 1, 15, 538000 ...

Мысль про поликлинику не понял (но вы продолжайте развивать). Где и что направильно делаю?


Изголяетесь. Вечереет...
Дело не "устраивает/неустраивает". Сейчас я истинное значение частоты знаю. Я пытаюсь найти алгоритм, который бы в короткое время выдал близкий к истинному значению результат. Все методы, что пока перепробовал, не дают результат с нужной точностью.

А почему, действительно, не взять '1'? Так можно было бы раз в мс просчитывать фильтр 1000 раз вместо одного сейчас. И результат был бы быстрее получен.

А Вы отсортируйте отсчеты, а потом фильтруйте.... Автору нужно получить константу ОДНУ на выходе.


Как вам везет - Вам ведома истина. Тогда посчитайте (Вы же и дисперсию знаете), сколько отсчетов (минимум) нужно для получения истины, лежащей внутри заданного интервала с заданной вероятностью.

Зачем? Наоборот, это же здорово, что "палки" выскакивают изредка. Не обращать на них внимания и призван медианный фильтр.


Так потом фильтром НЧ и получить ОДНУ.


+1
Зачем так неэффективно обходиться с данными? Если, как видим, простое среднее значение не дает нужного результата, а нужны дополнительные меры, вроде медианного фильтра, то, конечно, лучше обрабатывать каждый отсчет.

Если Вы отсортируете, то Ваш медианный фильтр не будет фильтровать Ваши "палки" с некоторой вероятностью. Результат должен быть инвариантен к перестановкам. Нет никакой последовательности - есть просто набор (выборка) данных. Так я понимаю... Поэтому - какие фильтры, фурье-преобразования....

Соображения были такие.
Допустим измеряется частота входного сигнала ("чистого", без шумов).
Измеряемая частота f₀ (входной сигнал). Частота опорного генератора f_изм. Интервал между N импульсами входного сигнала равен М импульсов опорного генератора.
Значит f₀/N = f_имз/(М +- 1).
Значит относительная погрешность измерения составит +-1/М.

Фактически это что-то типа выбора разрядности АЦП.

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



Истина мне неведома.


Не понял
Типа сколько кавалеристов в год будет убито копытом лошади?


P.S. Попробовал обработать последовательность:
1) путем отбрасывания пиков
2) путем замены пиков соседним "более достоверным" значением.
Улучшения как бэ... есть в плане "грязи" выходного сигнала, но результат получается ниже ожидаемого. Ошибка короче получается.

Ничему ничто не должно. Есть конкретная задача, которую надо решить.
А если я синус оцифрую, а потом перемешаю отсчеты, тоже должно потом фильтроваться, и БПФ частоту должно находить?


Называется "рандомизация". Был, есть и будет есть.

Это-то я понял. Если по усреднению 1000 периодов получаются выбросы на 10^-4, значит из этой тысячи один крив на 10^-1. Или два - на 5*10-2. Или ...

У Вас на входе уже нечто среднее из шумов, выбросов и полезного сигнала. Оно и было названо "средней температурой по больнице". Может быть есть возможность устранять выбросы ещё в первой тысяче?

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

Пустые слова какие-то. Задача, которую.... Нет там синуса, нет никакой функции ни от чего. Есть набор чисел. Нужно его статистически обработать. И результат должен быть инвариантен по отношению к перестановкам. Иначе мы не решим конкретную задачу, которую нужно решить.


Рандомизация, говорите... Если добавить к "сигналу" шум, уменьшим ли мы дисперсию? Эта ария из другой оперы. Информационный шум.

Пустые слова?
Ну почему же статистически? И почему результат должен ... обоснуйте!

Не дает мне подредактировать свое предыдущее сообщение. По-новой напишу.


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

Имеется входной сигнал примерно 1 kHz, и опорная частота 100 MHz. Формируйте интервал по 10 периодов или больше (но зачем по 1000?) и заполняйте 100 мегагерцами. А дальше фильтруйте все подряд. Похоже, что частота у вас меняется быстро (часто), может, поэтому трудно усредняется. Уже наложение спектра произошло.


А мы не дисперсию ищем. Математическое ожидание.

Это Вы обоснуйте. Автор пишет (писал), что он точно знает частоту. Ну... это он так думает. Но он близок. Это означает, что он считает (, что знает) априорно, что от времени (номера) ничего из того, что его интересует не зависит. Поэтому от перестановки тоже. Что еще нужно пояснять?


Так Вы считаете, что увеличивая дисперсию, мы точнее получаем математическое ожидание? Очень свежо! Обоснуйте.
Ведь когда мы вычисляем дисперсию (для чего!?), то мы потом гордо говорим - среднее значение равно... с такой-то вероятностью..

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

А почему по 10?

А как можно восстановить несущую ФМ-сигнала, когда модулирующий сигнал - шум с медленно меняющимся средним?

Пусть пишет... Шум фазы... А измеряет период... не один. Что такое шум фазы для непериодического сигнала... или периодического. И даже если он (кто, что?) имеет постоянную составляющую (Вы это понимаете?), которая еще медленно (относительно чего?) меняется (каким же образом?). Вот это и есть информационный шум. С огромной дисперсией.

Точно я частоту не знаю, никогда не знал и знать не буду.

Я хочу, чтобы мой алгоритм выдавал результат, находящийся в интервале (a;b). Этот интервал для конкретного устройства, определен экспериментально.
Когда результат будет стабилен, а проверю что и с другими устройствами результат тоже стабилен. Для каждого экземпляра ручками определю допустимый интервал, потом убедюсь что алгоритм выдает результат из этого интервала.
Вот как-то так.
Пусть не научно. Надо чтобы просто кое-как работало.

Мне не надо "среднее равно х с вероятностью y".
Мне надо "среднее находится в интервале от a до b c вероятностью 1".

Вам нужно 0.5 ppm. Если заполнить 10 периодов, то числа будут около 1000000. Если 20 - 2000000.
Вот и получается, что единица младшего разряда соответствует требуемой точности. Примерно, так.
Здесь от скорости изменения вашей частоты все зависит. Нужно успевать ее оцифровывать, чтобы удовлетворить теореме Котельникова. Иначе получается наложение спектра. Т.е. частота меняется быстро, а мы этого не видим, а видим какое-то медленное переползание. Или того хуже.

Метод этот для чего существует?!
Вы сами-то читали, что там написано?


Так так же оно и есть. Говорят, что истинное значение лежит с вероятностью ... в интервале одна сигма, с вероятностью ... в интервале 2 сигма. И даже можем сказать про вероятность того, что оно лежит в интервале полсигма. С некоторой вероятностью. А интервалы отсчитываются от абсолютно точно посчитанного среднего значения выборки. Мы его всегда с бесконечно хорошей точностью вычисляем. А с вероятностью 1... нельзя даже сказать, что... некто умрет через миллион лет.

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

Тогда проблему надо брать осадой.
1.Хотелось бы видеть гистограмму сигнала без выбросов хотя бы из 1000 точек, 20 градаций внутри Ваших границ.
2.Хотелось бы видеть спектр выборки из 2048 точек.
Это позволит оценить максимально возможную точность.

Сделайте так: все выборки, превышающие пределы возможного изменения постоянной составляющей + 3 сигмы шума (о нём позже) исключите из расчёта. От оставшихся найдите среднее арифметическое. Сигму шума вычислите по формуле выборочной дисперсии от отобранных отсчётов.

Автору: отсортируйте точки по амплитуде, отбросьте 3% самых больших и 3% самых маленьких точек, а для оставшихся точек усредните амплитуды. Это и будет ваша постоянная составляющая. Такая операция называется ранговым фильтром — это один из тех нелинейных фильтров, про которые вам говорит Таня.

P.S. Рандомизация (дитеринг) нужна в процессе оцифровки. А раз данные уже оцифрованы, никакой рандомизации не нужно.

Еще раз предлагаю автору темы перепроверить данные, если брать не среднее значение на интервале 1000 периодов, а средневзвешенное значение на подобранном окне непрямоугольной формы. Этим должны значительно уменьшиться выбросы.
Второе то, что каждый следующий результат должен получаться сдвигом окна не на 1000 исходных сэмплов, а на меньшую величину. Форма окна дает нам частотную характеристику и мы можем выбрать оптимальную децимацию. Возвращаясь к исходному окну на 1000 сэмплов, которое имеет частотку sin(x)/x, понимаем, что сглаженное окно значительно лучше подавляет шумы (перемалывая их в мелкодисперсионную кашу), чем прямоугольное. Нас не интересует сильно узкая полоса, а именно подавление шума в большей части полосы. Еще другими словами, у автора в исходной схеме слишком большой коэффициент децимации. Можно даже просто взять то же прямоугольное окно в 1000 сэмплов с шагом скольжения 100 (а не 1000), например, сэмплов, и уже будет улучшение или в уровне шума или в скорости.
**
т.е. коротко:
1. сменить окно предварительного фильтра с прямоугольного на сглаженное (предлагаю применить гауссовское)
2. уменьшить децимацию с выхода предварительного фильтра

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

Почти оффтопик. Вы меряете фазу + фазовый шум. Передо мной листовка "Фазовый шум кварцевых генераторов компании Мэджик Кристал". на каждом графике есть белый шум + низкочастотный шум. Граница около 100Гц. Причем низкочастотный шум бывает 30 раз на декаду и бывает более 100 раз на декаду. Как получается так много не знаю. А у Вас шум белый или низкочастотный?