На самом деле сейчас в качестве временного решения так и сделал. Но минусы этого дела я прочел в этойтеме, об этом тут и Designer56 говорит.


Не совсем понял в какой цифре? Я и так в цифре все это делаю.
Если вдруг остались непонятные моменты, давайте я еще обрисую картину, хотя вроде все важный моменты упомянул.
Итоговая цель - получить ФВЧ, но особенность в том, что частота дискретизации на несколько порядков выше частоты среза. Для этого я исходный сигнал децимировал в 25 раз, не забыв предварительно про антиалиасинговый ФНЧ и далее применил еще один ФНЧ на нужную мне полосу. Так я выделил постоянную составляющую, но частота дискретизации этой составляющей стала в 25 раз ниже исходной. Хочу поднять ее обратно, чтобы эту постоянную составляющую вычесть из исходного сигнала, для этого реализую интерполятор.
Все фильтры КИХ с линейной фазой. ФНЧ в интерполяторе рассчитан на полосу пропускания [0; fd] где fd это доли Герца, и частоту затухания Fs/(2 * 20) где Fs это исходная (большая) частота дискретизации, а 20 коэффициент интерполирования. Порядок фильра получился около 170! Так что полагаю довольно крутой.
Проблема то все в том, что фильтр давит этот переходной процесс (в данном случае переходной процесс появляется всегда, в течении всего времени работы, а не только в начале) настолько насколько рассчитан, но вот ввиду того, что входной скачок (0 минус мой постоянный сигнал) величина много большая, чем допустимые искажения, то этого подавления недостаточно. Можно конечно еще заузить фильтр, но он итак весьма неслабого порядка.

Непонятно только, нафиг это нужно, если стоит задача определить постоянную составляющую...


ЗЫ. Ах, простите, в аккурат сейчас же всё обрисовали. Теперь будем думать.

Думаю, в данном случае можно попытаться обойтись без прореживания.
Соорудив простейший адаптивный фильтр на частоте дискретизации.
Разностное уравнение:
Y(k)=a*X(k)+ (1-a)Y(k-1),
где Y(k) - текущая оценка постоянной составляющей, X(k) - текущий входной отсчёт, Y(k-1) - оценка постоянной составляющей на предыдущем шаге, a - весьма малый коэффициент, определяющий динамику системы.

Инициализировать оцениватель нужно так:
Y(0)=X(0).

Недостатком является необходимость применения повышенной разрядности вычислений, если есть желание получить несмещённую оценку, т.к, в этом случае, усекать разрядность X(k) нежелательно.

Сообщение отредактировал Stanislav - Nov 5 2010, 19:37

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

С разрядностью проблем нет, все работает на ПК с использованием double либо single.
Буду думать. А не будет проблем в применении тут БИХ фильтра? Я вот не соображу требование линейности фазы при удалении постоянной составляющей необходимо?


Да, полезный сигнал действительно может быть децибел на 80 по амплитуде меньше постоянной составляющей.
Насчет грубого вычитания как раз и думал в первую очередь.
С другой стороны, чем плоха вставка не нулей, а предыдущих значений в моем случае... ну и пусть АЧХ сигнала с постоянной составляющей будет с дополнительным подавлением sin(x)/x, от этого постоянная составляющая не перестанет быть ей, единственное что наверное расчетная частота среза фильтра пострадает.
Спасибо всем, буду думать над предложениями!

Так это будет помеха, а не подавление, и она будет маскировать полезный сигнал. Вам нужен фильтр с затуханием 80дБ + допустимое отношение сигнал/шум

в своих модемах, в качестве корректора постоянки, использую фильтр "памяти" канала для Transform Domain эквалайзеров, взятый из книг по адаптивной фильтрации. Суть в слудующем :
Вычисляем постоянку как y = (1/2^a)*x + (1-1/2^a)*y, a - однозначно определяет полосу, затем, без каких либо преобразований частот вычитаем постоянку из сигнала.
Много позже нашел описание такого подхода на freepatents %) Минусы данного подхода большое время переходного процесса, но при работе линейных модемов это не сильно актуально. Зато метод простой и весит копейки.

Простите, но мне не совсем понятен Ваш подход.
Плясать нужно от определений. Зачем Вам какая-то децимация, если нужно найти оценку постоянной составляющей (сиречь, матожидания - МО) сигнала?
Если волнует проблема нелинейности фазы БИХ-фильтра, соорудите КИХ-фильтр на частоте дискретизации. Оптимальным в смысле СКО оценки МО на заданной последовательности длины N является просто усреднение её элементов, которому в случае непрерывного поступающих отсчётов соответствует фильтр скользящего среднего длины N.

(1)

Инициализировать память такого фильтра можно значением X(0) - переходный процесс при этом будет существенно уменьшен. Можно придумать и более "продвинутые" способы, обеспечивающие лучшую сходимость к оптимальному значению.
Такой фильтр не содержит умножений при накоплении, а если выбрать N равным степени двойки, операция умножения на коэффициент 1/N реализуется при помощи сдвига с округлением.
Если требуемая длина фильтра получается большой, количество операций в нём можно радикально изменить, преобразовав его в рекурсивный вид. Для этого, положим

(2)

Убедиться в эквивалентности (1) и (2) могут все желающие.
Фильтр содержит только 3 операции сложения/вычитания и одну операцию сдвига (умножения на к-т) на отсчёт, и это количество не зависит от его длины! Промежуточный результат накопления Y(n) должен содержаться в регистре повышенной разрядности, дабы получить несмещённую оценку.
Способ инициализации такого фильтра попробуйте предложить сами.

Если рассчитать КИХ фильтр на частоте дискретизации (это, понятно, первая мысль, которая приходит в голову), то его порядок получается что-то вроде 4000-5000, если память не изменяет. Это потому что отношение частоты среза к частоте дискретизации 10^-4

Ваши сомнения мне непонятны...
Если память позволяет, пуркуа бы и не па?
Рассчитывать там особенно ничего не надо: все коэффициенты равны 1 или -1. Проще пареной репы, да ещё и оптимально к тому же...

Фильтр скользящего усреднения со всеми единичными коэффициентами имеет довольно паршивую АЧХ в области пропускания и, особенно, в области подавления: первый лепесток АЧХ всего лишь 1/5 от основного, да и остальные весьма хреново спадают. Ситуация с подавлением улучшается последовательным включением нескольких таких фильтров, но опять таки, это очень далеко по характеристикам от оптимального фильтра, хотя и требует гораздо меньше вычислений.

Что касается основного вопроса автора, то, как уже было замечено, Вы, скорее всего где-то ошибаетесь. Возможно у Вас слишком короткий фильтр, или же неправильно выбрана задержка сигнала. Есть еще один момент: если Вы хотите обеспечить подавление низких частот на уровне 80дБ, то все ФНЧ в дециматоре и интерполяторе должны иметь пульсации в области пропускания такого же порядка! К подавлению требования наверно еще более жесткие.

В прикрепленном файле простейшая модель для Simulinkа, только я немного уменьшил Ваши требования к соотношению частот среза и дискретизации: Fd = 1MHz, Fpass = 500Hz, Fstop = 1500Hz. На картинке результат моделирования:

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

Свят-свят, братцы... Ей-богу, до несварения желудка довести способны...
При чём здесь АЧХ, если нужно произвести оценку постоянной составляющей сигнала? ДУМАЙТЕ прежде, чем отвечать и приводить графики.

Проторенных дорог не существует. Но рано или поздно всё получится.

Это не такая уж и 'постоянная составляющая', хотя так называется для разговора. И так как она меняется во времени, у нее есть спектр, компоненты которого на частоте сигнала должны быть на (скажем) 100 дБ ниже нулевой компоненты. А сделать, чтобы она не менялась по времени, это надо взять всю выборку целиком, сложить и вычесть.

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

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

3. Если вам надо вырезать какую-то часть спектра вблизи нуля, то для этого необходимо ставить фвч. При этом частота среза вашего цифрового ФВЧ будет зависеть от полосы вблизи нуля и частоты дискретизации. Если скажем вы хотите на частоте дискретизации 1 МГц отфильтровать все что ниже 1 Гц, то вам потребуется FIR фильтр ОЧЕНЬ большого порядка. IIR в этом случае вообще может быть неустойчив даже при 32 битном представлении с плавующей точкой. Это обусловленно самым фундаментальным принципом - принципом неопределенности и ни один смертный ничего не может с этим поделать.

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

4. Не понятрно, почему вы боитесь исказить фазу IIR фильтром? если у вас нулевая пч то ставьте один и тотже IIR на I и на Q и наслаждайтесь. Никакие характеристики вы не потеряете, потому что разность фаз квадратурных каналов в останется неизменной.

Все правильно.
Только хочется заметить, что никаких магических универсальных формул не бывает.
Без априорной информации об "истинном" сигнале любые фильтры не имеют смысла.
Вот имеются "измерения" в какие-то моменты времени. Без дополнительных знаний или предположений они и являются сигналом. А если я предполагаю, что должна быть линейная функция, то должна провести прямую (линейная аппроксимация) максимально близко к этим точкам. Или другой функцией приблизить. Тут еще нужно знать, какова ошибка "измерений".
Вот пусть есть 4 точки - (0, 0), (1, 1), (100, 1) (101, 0). Можно приблизить квадратичной функцией... А, может, это горизонтальная прямая?
Извините за банальности.