вопрос про интерполяцию Опции

[jorikdima]

Здравствуйте.
Возник вопрос, который быть может, надуман, но хочется разобраться.
Есть сигнал - результат НЧ фильтрации исходного сигнала с целью определить его постоянную составляющую. То есть полоса сигнала в общем то составляет доли Герца, но среднее значение (или нулевая гармоника) может варьироваться в довольно большом диапазоне. Частоту дискретизации этого сигнала нужно поднять в 25 раз, для чего требуется провести интерполяцию. В книжках прочитал, что делается это с помощью добавления нужного числа нулей между исходными отсчетами и пропусканием через ФНЧ для удаления размноженных гармоник. В принципе все работает, но появляются артефакты. Дело в том, что во всех книжках, где рассказывается про интерполяцию, нарисован сигнал с околонулевой постоянной составляющей, ну либо нулевая гармоника много меньше других гармоник в сигнале. И тут все прекрасно, добавляются нули, затем они сглаживаются ФНЧ. Но в случае когда, нулевая гармоника является преобладающей добавление в сигнал отсчетов, сильно отличающихся от текущего значения сигнала (а именно нулевые отсчеты) каждый раз вызывают всплеск на выходе фильтра, пропорциональный разнице между нулем и реальным значением сигнала и зависящим от ИХ фильтра. Этот всплеск у меня много больше шума исходного сигнала, являющейся собственно постоянным током и недопустим.
Что делают гуру в таких ситуациях?
У самого есть идея следующая. Не добавлять нулевых отсчетов и не ставить ФНЧ, а вписывать какой-нибудь полином, тот же сплайн, в реальные отсчеты сигнала, а затем, найдя коэффициенты сплайна, рассчитывать 24 значения между исходными выборками. Минус такого решения, что это не риал тайм, нужно дождаться следующего отсчета прежде чем получить 24 точки до него, но мне не страшно это.
Справедливо ли решение или есть какие-то более мудрые действия?
Спасибо.

[jorikdima]

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

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

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

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

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

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

[Stanislav]

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

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

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

(1)

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

(2)

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

[jorikdima]

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

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

[Stanislav]

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

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

[Самурай]

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

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

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

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

Прикрепленные файлы

 LP_FIR_Decim_Interp.zip ( 20.01 килобайт ) Кол-во скачиваний: 18

 

[Stanislav]

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

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

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

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