Рекурсивный фильтр скользящего среднего



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

    Как ты уже наверняка знаешь, дорогой читатель, существует два способа построения цифровых фильтров. Это рекурсивные фильтры, они же фильтры с бесконечной импульсной характеристикой (БИХ), и трансверсальные фильтры, они же фильтры с конечной импульсной характеристикой (КИХ). Самым простым и широко используемым фильтром КИХ является «фильтр скользящего среднего». Результат фильтрации такого фильтра, есть среднее арифметическое последних N отсчетов входного сигнала.



    Или, в развернутом виде, для N=4:



    Функция на языке С реализующая фильтр скользящего среднего:

    #define N (4)
    int filter(int a)
      {
      static int m[N];
      static int n;
      m[n]=a;
      n=(n+1)%N;
      a=0;
      for(int j=0;j<N;j++){a=a+m[j];}
      return a/N;
      }
     

    Комплексный коэффициент передачи фильтра скользящего среднего, нормированный относительно частоты дискретизации, определится как преобразование Фурье от импульсной характеристики:



    График амплитудно-частотной характеристики (АЧХ), нормированной относительно частоты дискретизации, при различных значениях длинны фильтра (N=4;8;16), приведен на рисунке:



    Соответственно, график фазо-частотной характеристики (ФЧХ):



    Данный фильтр нашел широкое применение в обработке сигналов, отчасти благодаря своей простоте, но самое главное его свойство — линейная фазо-частотная характеристика, и, соответственно, постоянное во всей полосе частот время запаздывания сигнала. Этот фильтр трансформирует амплитудный спектр сигнала, не затрагивая фазовый, что делает удобным его использование в системах регулирования. Фильтр скользящего среднего, благодаря своей линейной переходной характеристике, широко применяется при линейной интерполяции, передискретизации сигнала и т.п.

    Главный недостаток фильтра скользящего среднего — вычислительная сложность, пропорциональная длине фильтра N. Для решения этой проблемы существует рекурсивный фильтр скользящего среднего. То есть, фильтр, имеющий те же характеристики, что и классический фильтр скользящего среднего, но реализованный по рекурсивной схеме. Такие типы фильтров широко известны в узких кругах, и называются: рекурсивные фильтры с линейной ФЧХ [Введение в цифровую фильтрацию. Под. ред. Богнера Р. М: 1976] или CIC-фильтры [DspLib]. Существует научная школа проф. Турулина И.И. [РГБ], занимающаяся исследованием подобных фильтров.

    Покажем способ построения рекурсивного фильтра скользящего среднего на примере фильтра длинной N=4. Впоследствии нетрудно будет обобщить результаты на произвольную длину фильтра.

    Как было отмечено выше, значение n-ного отсчета сигнала на выходе фильтра можно определить как:



    А значение предыдущего, ((n-1)-го) отсчета:



    Вычтем из первого выражения второе, в результате получим:



    Нетрудно сообразить, что при произвольной длине фильтра N, уравнение запишется в следующем виде:

    (1)

    На основании приведенного уравнения можно записать код фильтра на языке С, но сначала мы выполним проверку наших выкладок. Найдем частотные характеристики рекурсивного фильтра скользящего среднего, для чего выполним Z-преобразование уравнения фильтра. Хочу напомнить дорогому читателю, что для выполнения Z-преобразования необходимо заменить переменные (xn,yn) их Z-отображениями (X,Y), каждое понижение индекса переменной на единицу соответствует умножению Z-1.



    Решим полученное уравнение относительно Y/X –коэффициента передачи фильтра



    Перейдем из Z-области в частотную область, заменив на , и получим комплексный коэффициент передачи:



    Построим график модуля комплексного коэффициента передачи (АЧХ фильтра), при различных значениях N=4,8,16:



    А также аргумент комплексного коэффициента передачи (ФЧХ фильтра):



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

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



    В полученном выражении выполним подстановку:




    В результате уравнение фильтра (1) преобразуется в систему уравнений:




    На основании системы уравнений фильтра запишем код, реализующий фильтр на языке С.

    #define N (4)
    int filter(int x)
      {
      static int n;
      static int m[N];
      static int y;
      y=y+(x-m[n]);
      m[n]=x;
      n=(n+1)%N;
      return y/N;
      }

    Стоит отметить, что любители «совершенного кода» могут наложить ограничивающие требования на длину фильтра N. При длине фильтра, равной степеням двойки (2,4,8,16…), можно заменить операцию деления (/) арифметическим сдвигом (>>), операцию остаток от деления (%) — побитовой конъюнкцией (&).

    Заключение

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

    Успехов тебе, дорогой читатель!
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More
    Ads

    Comments 45

      +2
      Ещё одна очевидная проблема рекурсивного фильтра — малейшее однократное повреждение данных приводит к перманентному искажению результата, а оно в реальных системах может произойти по разным причинам — космическое излучение, несанкционированная запись в область памяти массива со значениями и т.д. поэтому для надёжности иногда приходится пересчитывать по полной весь массив. Или должна быть какая-то обратная связь которая устранит расхождение.
        0
        Селяви, как говорят французы
          0
          В рекурсивном варианте некоторое отклоение d на шаге k даст вклад в выход фильтра +d/N, но потом на шаге k+N это же отклонение вычтется.
            0
            Если значение в массиве изменится несанкционировано, то вычтется уже совсем другая величина…
              0
              (удалено, дубликат)
              0
              В системах, у которых есть эти самые «малейшие однократные повреждения данных», проблем куда больше.
              Потому что вероятность, что слетит стек, какой-нибудь индекс или прочий указатель ровно такая же.
              Нет, я знаю про ватчдог, но это есть костыль — система перед его срабатыванием имела полное право делать чёрт-знает-что.
                0
                А какие еще есть костыли для систем с возможным повреждениями данных?
                Из программных средств для контроля целостности стека вызовов можно каждой функции передавать параметр CalledById и внутри функции его проверять. Желательно иметь поддержку такой фичи на уровне компилятора.
                  0
                  Вариант, в котором можно быть уверенным — дублировать (троировать) системы управления.

                  А «иногда пересчитывать» этот несчастный фильтр — это из анекдота о поиске кошелька под фонарём.
                    0
                    Троирование системы управления не поможет, т.к. программа на них будет одна и программная ошибка сделает плохо на всех экземплярах. Она в любом случае сделает плохо, даже если всё работать будет идеально — как только используешь арифметику с плавающей запятой — будет накапливаться неустранимая ошибка из-за неизбежного округления.
                      0
                      У меня тут вопрос возник. У Вас какая проблема?
                      Космические излучения? Дублируйте. Иначе за то время, что пройдёт между прилётом частицы и очередным пересчётом фильтра Ваш космолёт включит самый-главный-двигатель и улетит чёрт знает куда.
                      Программные ошибки? Пересчитывать каждые 5 секунд одни и те же данные с одной и той же ошибкой… Не лучшее решение.
                      Накопление ошибки во флоатах? Всё просто — не делайте так никогда. Разговор не с этого начался.

                      Или просто поговорить надо на отвлечённые темы?.. Тогда, пожалуйста, не ко мне.
                        0
                        Дело в том что однократная ошибка в случае полного пересчета уйдёт со временем и больше беспокоить не будет, а для рекурсивного алгоритма ошибка засядет надолго и если это девайс с аптаймом в 15 лет, за такой срок ошибок накопится просто вагон. А если его перезапускать каждые сутки никто ничего не заметит, конечно. Вот для этого и нужны пересчёта по полной, чтобы периодически сбрасывать возможные ошибки и не допускать их накопления. Независимо от их природы.
                  0
                  По большей части, повреждения данных не приводят к перманентным искажениям данных, даже повреждение стека — поток рано или поздно будет завершен аварийно, но система продолжит работать. А тут данные, обнаружить изменение которых нет никакой возможности а последствия — перманентные искажения, на основе которых будут приниматься какие-то важные решения.
                  Вач-дог спасает совсем от других проблем.
                  Система должна обладать свойством самовосстановления после сбоя, даже если часть данных будет потеряна она должна продолжить работать.
                0

                "Главный недостаток фильтра скользящего среднего — вычислительная сложность, пропорциональная длине фильтра N." — ну надо ж было ляпнуть такую ерунду. Вычислительная сложность — одно сложение, одно вычитание, одно деление на константу (заменяется на умножение, а то и на сдвиг). O(1).
                Вот памяти O(N) надо, это да. Поэтому на совсем мелких девайсах непопулярен. И других недостатков хватает.

                  +1

                  Автор имел в виду нерекурсивный фильтр (т.е. функция от предыдущих N входных отсчетов) — а он реализуется как сумматор-делитель всех отсчетов, со сложностью O(N).

                    0
                    Сложить, вычесть — это и есть «открытый» автором рекурсивный фильтр.
                    А изначальный фильтр у него — все сложить.
                    Причем у автора в программе ошибка.

                      0
                      Простите… где?
                        0
                        полагаю в этой строке
                        m[i]=x;
                        так как до этого i не объявлялась.
                          0
                          Прошу прощения, поправил. Подгонял код по формулы и пропустил.
                        +2

                        В следующей статье автор расскажет, как можно 2+2+2+2+2 заменить на 2*5?
                        При этом приведя кучу формул и наукообразной терминологии?

                      0
                      Зачем переменные static? это не позволяет делать в программе 2 параллельных фильтра.
                        0
                        Потому что это С11, переменные с модификатором static объявленные внутри блока, имеют глобальное время жизни и локальную область действия
                          0
                          Прошу прощения описался C99
                            0
                            Это я знаю.

                            мне параллельно надо делать 2 фильтра. Но вызывать вашу функцию я не могу. Она не реентерабельна.
                              0
                              В чем сложность с двумя фильтрами?
                              опишите 2 функции или объявите void filter(int y, int x, int* m, int n);
                          0
                          Чем же рекурсивный вариант лучше прямого? Все равно N значений в памяти хранить
                            0
                            но вычислять то одно
                              0

                              тогда зачем хранить n отсчетов?

                                0
                                Нда, без буффера никак, по сути надо реализовывать задержку в N отсчетов. Чего кстати в последней реализации на си не видно, думаю так будет корректнее:

                                #define N (4)
                                int filter(int x)
                                {
                                static int m[N];
                                static int y;
                                y=y+(x-m[0]);
                                for (int i=1;i<N;i++) m[i-1]=m[i];
                                m[N-1]=x;
                                return y/N;
                                }
                                  +1
                                  Ужасно-ужасно! Двигать массив поэлементно? Это будет дольше чем считать «по полной программе».
                                  Вместо того чтобы двигать «окно» лучше организовать классический кольцевой буфер — двигается указатель, а массив остаётся на месте. Число по указателю извлекаем — это будет последний элемент буфера, новое на его место записываем — это будет начало буфера, и сдвигаем указатель по кругу на 1 позицию. Никаких перекидываний целых кусков памяти!
                                  Представьте себе крайнюю ситуацию, когда N = 2000000 и каждую итерацию вам надо будет 2млн элементов сдвинуть на одну позицию… да быстрее их будет сложить и не парится с рекурсивным вычислением.
                                    0
                                    оно може и понятнее, но работает дольше

                                    «Чего кстати в последней реализации на си не видно»
                                    по моему отлично видно «n=(n+1)%N;»

                                    я так понимаю большинство собравшихся тупит по патерну КИХ фильтра
                                    Поясню. Входные данные складываются в буфер, указатель буфера закольцовывается n=(n+1)%N;
                                    Значение указателя и есть нулевой отсчет, или согласно математическим формулам n-ный отсчет
                                    Минус первый (n-1) отсчет это (n+N-1)%N, или в условиях N=2^L (n-1)&(N-1)
                                    Минус k-тый (n-k) это (n+N-k)%N при условии что k<N, отсчеты k>N не сохранены
                                      0
                                      > указатель буфера закольцовывается n=(n+1)%N;

                                      Да точно, сначала не разглядел.
                                0
                                Прямой — надо складывать все N значений, потом делить…
                                А рекурсивный — сложить текущее значение(начало окна) и вычесть последнее(конец окна), и одно деление.
                                  0
                                  Вероятно этот алгоритм медленнее так как предусматривает дополнительную запись в память.
                                  Думаю сложение 4х будет быстрее
                                    0
                                    Это какую дополнительную? Обычно в таких алгоритмах используют кольцевой буфер — по указателю считали ячейку и это будет наш последний элемент, и на его место вписали свежее значение, передвинули указатель на следующий элемент или перешли на начало если это конец выделенной области под массив.
                                      0
                                      N может быть куда больше 4. Просто пример приведен для N = 4.
                                        0
                                        Во-первых, тут ОДНА дополнительная запись. В нерекурсивном оригинале будет 4 (а в большинстве случаев куда больше, скажем 64) чтения из той же памяти.
                                        А во-вторых, это микроконтроллер. Если ОЗУ набортная (есть «толстые» модели с внешней памятью, но во всяких DIY это редкость), доступ к ней — ровно один такт.
                                          0
                                          Заметьте, ОДНА дополнительная ЗАПИСЬ.
                                          Так как есть требование когерентности кеша то данная операция отразиться на замедлении всех ядер.
                                          Сейчас основные тормоза это совсем не сложение или умножение — это чтение и запись в память.

                                          доступ к памяти десятки и сотни тактов.
                                            0
                                            Когерентность кэша нужна, если один фильтр считается из разных потоков. Так сделать можно, но зачем?
                                            ПС. Если считать строго, то записей в память больше:
                                            не рекурсивный) чтение и запись индекса, запись в массив, N ЧТЕНИЙ из массива
                                            рекурсивный) чтение и запись индекса, ЧТЕНИЕ и запись в массив, ЧТЕНИЕ и ЗАПИСЬ суммы.
                                            В первом случае мы задействуем N+1 адрес памяти за раз, во втором 3
                                            Чем больше N тем больше промахов кэша будет в первом случае. Во втором количество промахов почти не зависит от N.
                                    0
                                    А как же CIC фильтр и Hogenauer pruning для пущей экономии?
                                      0
                                      так этож оно и есть

                                      Такие типы фильтров широко известны в узких кругах, и называются: рекурсивные фильтры с линейной ФЧХ [Введение в цифровую фильтрацию. Под. ред. Богнера Р. М: 1976] или CIC-фильтры [DspLib]. Существует научная школа проф. Турулина И.И. [РГБ], занимающаяся исследованием подобных фильтров.
                                        0
                                        Я бы поостерегся называть CIC фильтр рекурсивным. В нем отсутствует обратная связь по выходу — отличительный признак.
                                      0
                                      Можно пойти дальше и таким же образом реализовать треугольный фильтр.
                                        0
                                        Это здесь [Введение в цифровую фильтрацию. Под. ред. Богнера Р. М: 1976]
                                        +1
                                        Ещё такой фильтр неустойчивый — при реализации на числах с плавающей точкой он через некоторое время станет накапливать постоянную составляющую.
                                          0
                                          Вот как раз в предложенной реализации такого не происходит

                                        Only users with full accounts can post comments. Log in, please.