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

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

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

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

Сообщение отредактировал R.A.K. - Sep 24 2011, 16:13

К сожалению не все выбросы компенсируются противоположными, поэтому фильтрация фнч будет с медленными выбросами.
Единственное хорошее свойство Вашего сигнала это то, что выбросов в выборке немного.
Главное для каждой точки определить выброс это или нет и в случае выброса заменить ближайшим невыбросным значением.
Для этого достаточно знать приближенное значение измеряемой в настоящее время величины и дисперсию присутствующего в настоящее время шума.
Все значения за пределами трех сигм следует считать выбросами.
Таким образом предлагаемый алгоритм такой: Берем выборку 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) отсчетов. На выходе сигнал весьма грязненький.

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

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