Честно простой цифровой фильтр

Вы работаете с АЦП. Получаете результаты преобразования, один за одним. И замечаете, что эти результаты «скачут». А хотелось бы, чтобы стояли, как… Ну, короче, чтобы стояли!
Есть много причин, почему отсчеты АЦП могут быть нестабильны. В своей заметке я не говорю об этих причинах. Я говорю о том, как успокоить показания, получая их AS IS. И как сделать это максимально просто. При этом, возможно, не имея ни малейшего понятия о науке под названием «цифровая  обработка сигналов».
Эта заметка написана в качестве полной замены предыдущей заметки. Ту лучше не читать :)

ПОСТАНОВКА ЗАДАЧИ

Имеется последовательность кодов, дискретно во времени представляющих физический сигнал. Мы будем говорить о последовательности кодов с АЦП. Физический сигнал и полученная последовательность кодов имеют шумы, выбросы, «болтанку» — назвать можно как угодно. Наша задача — сгладить входную последовательность, то есть, выдать выходную последовательность такую, чтобы влияние шумов было уменьшено.
При этом мы стремимся выполнить задачу максимально простыми программными средствами обычного микроконтроллера.
Более того, у нас поставлено еще одно условие. Допустим, что нам неудобно накапливать данные, а потом обрабатывать их и выдавать результат один раз на N полученных кодов. То ли буфер для данных негде организовать, то ли темп выдачи результатов должен совпадать с темпом получения кодов от АЦП, то ли еще что. Но условие поставлено.
И тут нам на помощь приходит цифровой рекурсивный фильтр, самой простой реализацией которого является фильтр первого порядка. Вот его и будем делать.

ТЕРМИНОЛОГИЯ

Мы хотим реализовать блок программной обработки кодов АЦП, на который поступают входные отсчеты, обрабатываются и преобразуются в выходные отсчеты. Вот этот блок обработки и будем называть цифровым фильтром:


В каждый из моментов времени t1-2Т ,  t1-Tt1t1+T  и так далее на входе появляется очередной отсчет последовательности Х, а цифровой фильтр выдает новый отсчет последовательности Y. Если на каждый входной отсчет фильтр выдает и выходной отсчет, то это и есть классический цифровой фильтр.
Иногда поступают не так: набирают некоторое количество входных отсчетов и затем обрабатывают их. Таким образом, на несколько входных получают один выходной отсчет. Это, как правило, уже задача получения того или иного интегрального значения (например, среднего). Она настолько распространена при фильтрации шумов, что часто воспринимается как «обычная». И поэтому классический цифровой фильтр, выдающий отсчеты с каждым входным, называют "скользящим усреднением", как нечто не совсем обычное для усреднения. Что ж, значит мы рассматриваем скользящее усреднение. Но понимаем, что это обычный цифровой фильтр :)
Важной особенностью цифрового фильтра является то, что входные отсчеты приходят с постоянным темпом, в нашем случае — с периодом Т. Обратная величина называется частотой дискретизации и играет огромную роль во всех параметрах цифрового фильтра. Самое простое, что необходимо себе четко представлять: частотная характеристика цифрового фильтра симметрична относительно Fд/2 и полностью повторяется за пределами . Из чего следует, что фильтр нижних частот, который мы проектируем сейчас, не способен подавлять сигналы с частотой , 2Fд, 3Fд, и т.д. Этого понимания пока достаточно для нашей задачи.
Математически цифровой фильтр 1-го порядка описывается различными способами. Мы будем использовать такое представление:

Y(n) = Alfa*Y(n-1) + Beta*X(n)

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

Alfa + Beta = 1

В рамках данной заметки задача расчета цифрового рекурсивного фильтра 1-го порядка состоит в нахождении коэффициентов  Alfa иBeta с учетом удобства их использования в микроконтроллере (МК) для цифровой фильтрации отсчетов.

РАСЧЕТ ФИЛЬТРА

Цифровой фильтр (равно как и не цифровой) имеет как бы 2 лица: частотную характеристику и переходную (временнУю) характеристику. Задачу расчета можно ставить как получение требуемой частотной либо требуемой переходной характеристики. Я люблю исходить из заданного поведения во временной области. И объясню почему.
Дело в том, что частотная характеристика фильтра первого порядка… как бы это выразиться… простенькая. Возьмем 2 фильтра с частотой дискретизации 200 Гц. Первый с частотой среза по уровню -3 дБ, равной 5 Гц, второй — с частотой 1 Гц:



Как видим, в большей части частотной характеристики подавление шумов составляет всего 15...30 дБ, а у второго — 20..40 дБ. Вспомним, я отмечал выше, что за пределами Fд/2 (у нас это 100 Гц) характеристика симметрична и в районе 200 Гц подавления шумов снова нет.  Если нам нужно «взрослое зубастое подавление», то необходимо строить более серьезные фильтры.
Но все же, второй фильтр лучше давит помехи! А что, если частоту среза еще понизить?
И тут оказывается, что в погоне за парочкой децибел дополнительного подавления мы совсем забыли о временнОм «лице» фильтра. А давайте глянем, как реагируют эти два фильтра на ступенчатое входное воздействие (сигнала не было, а потом он вдруг стал равен 1 и таким и остался).



Что мы видим? Если говорить о времени переходного процесса, то первый фильтр (который похуже фильтрует) отрабатывает входной скачек примерно за 160 мс, а второй, наш передовик с подавлением 20...40 дБ — почти за 800 мс. Вот так-то. Лучше фильтрация — хуже переходной процесс. Поэтому и нужно выбрать некий оптимум.
Вот, понимая это, я и предлагаю: исходить из требований по быстродействию. Задавшись временем реакции на ступенчатое входное воздействие, мы получим параметры фильтра (вот те самые Alfa иBeta), а параметры фильтрации примем уж какие получатся. Парочка децибел туды-сюды уже мало что изменят, а быстродействие будет известно.

Вот почему в качестве исходного требования я выбираю обеспечение заданного времени переходного процесса при подаче на вход ступенчатого воздействия. Рассмотрим чуть пристальнее, как  же фильтр отрабатывает такой сигнал. Если дать скачек от нуля до некоторого значения, обозначенного 100%, то на выходе фильтра увидим:


Как видим, на 8-м отсчете фильтр уже отработал 90% входного воздействия, на 10-м — 95%, и так далее. Обычно принято говорить о «недоработанном», т.е. о тех 10 или 5 процентах, на которые фильтр еще «врет» к какому-то отсчету. Говорят еще, что погрешность установления выходного сигнала составляет столько-то процентов. Далее я привожу формулы для погрешности установления от 5 до 0,1%. В первом случае переходной процесс закончился «на глазок», а точности 0,1% обычно достаточно для того, чтобы считать процесс полностью законченным.
Разница между этими 5%-ым и 0,1%-ым фильтрами не столь уж велика. Я предполагаю, что все фильтры, которые будут разработаны по описываемой методике, находятся в континиуме между этими двумя крайними точками. В качестве характеристики фильтра по степени «законченности» переходного процесса введем такой параметр: Ntau — и вычислим его крайние значения:

LN(1/5%)  = 2,996
LN(1/0,1%) = 6,908
Ntau = 2,996...6,908

Смысл Ntau примерно таков: чем более жесткие требования мы предъявляем к завершенности переходного процесса, тем больше это число. Так что нижней границе значений Ntau  соответствует «грубый» процесс, когда мы спешим считать отработку законченной, а верхней границе — «точный» переходной процесс.

Итак, разработчик задался временем переходного процесса Тпп, за которое фильтр отработает скачек на 95...99,9%. Что еще нужно знать? Время выборки — тот период, с которым на вход фильтра поступают выборки и с таким же темпом с фильтра уходят результаты. Он у нас обозначен Т.
Ясно, что за время переходного процесса будет обработано много отсчетов. Сколько?

N = Тпп / Т                                                                                                               (1)

И наша задача — выбрать значения  Alfa иBeta так, чтобы уложиться в Тпп.
Оказывается, все очень просто. Должно выполняться условие:

Alfa < EXP(-Ntau / N)

Задавшись двумя граничными значениями Ntau, зная требуемое время переходного процесса и период выборки, мы находим 2 значения Alfa, соответствующие «грубому» и «точному» решениям. Выбираем значение Alfa  (ниже я покажу как) в диапазоне:

EXP(-6,908 / N) <= Alfa <= EXP(-2,996 / N)                                                      (2)

где <= означает «меньше или равно»
А потом легко вычисляем второй коэффициент:

Beta = 1 — Alfa                                                                                                          (3)

Вот и все. Коэффициенты фильтра вычислены. Поставив их в утилиту, описанную в заметке коллеги mc2, можете полюбоваться частотной характеристикой. И убедиться, что за требуемое вам время Тпп переходной процесс заканчивается.
Если уж говорить о программе для расчета, то можно обойтись и без «моего» метода. Но тут уж вам решать. Расчет по формулам (1)...(3) несложен, но гарантирует вам безытерационное решение для заданного переходного процесса.

ПРОГРАММНАЯ РЕАЛИЗАЦИЯ

Если мы работаем с МК, то очень хорошей практикой является использование целочисленной арифметики. Даже работа с 32-разрядными LONG переменными в 8-битном МК здесь выполняется намного быстрее, чем с FLOAT. Поэтому я рекомендую выбирать Alfa и Beta таким образом:

Alfa = Na / 2^k
Beta = Nb / 2^k
Na + Nb = 2^k

Иначе говоря, вместо дробных Alfa и Beta нужно взять целые числа, такие, чтобы их сумма была равна целой степени двойки: 2, 4, 8, 16, 32, 64, 128, 256, 512 и т.д.
Вот тогда процессору работы поменьше:

Y(n) = (Na*Y(n-1) + Nb*X(n)) >> k

Микроконтроллер умножает 2 раза целое на целое, складывает, производит сдвиг вправо на k разрядов. Несколько микросекунд — и у вас новый отсчет выходной величины.

ПРИМЕР ДЛЯ ЦЕЛОЧИСЛЕННОЙ АРИФМЕТИКИ 

Рассмотрим пример, в котором я задался N=100. Тогда условие (2) будет выглядеть так:

0,9333 <= Alfa <= 0,9705

Для целочисленной арифметики выберем однобайтную базу коэффициентов. Иначе говоря, пусть сумма Na и Nb будет равна 256. Тогда k=8 и получаем такие границы целочисленного эквивалента Alfa:

239 <= Alfa*256 <= 248
239 <=      Na      <= 248

Выбрав любое значение Na из этого диапазона, мы можем быть уверены, что за 100 выборок после подачи сигнала на вход выходной сигнал войдет в зону 0,1...5 % от установившегося значения. Хотим именно 0,1% — берем меньшее значение Na из указанного диапазона, допускаем 5% — берем большее значение.
И уже от выбранного значения Na расчитываем второй коэффициент:

Nb = 256 — Na

Теперь формула для работы микроконтроллера такова:

Y(n) = (Na*Y(n-1) + Nb*X(n)) >> 8                                                                  (4)

На этом расчет фильтра с целочисленными коэффициентами закончен.

Коллега _pv подсказал в обсуждении хорошую запись вместо (4), полностью эквивалентую по результату:

Y(n) = Y(n-1) + [Nb*(X(n) — Y(n-1)) >> 8]

Здесь имеем экономию на одно умножение (при произвольных соотношениях между коэффициентами) и формула приобретает вид приращения старого значения Y(n) на разность между входом и выходом, умноженную на коэффициент Beta. Чем меньше этот коэффициент, тем медленнее следит выход за входом — более сильная фильтрация.

Приведу сишный код, соответствующий формуле от _pv:

uint filter(ucha x, ucha Nb, ucha k){
  static ucha y = 0;
  static uint z = 0;
  z += (x — y);
  return y = (Nb * z) >> k;
};

Здесь учтено, что при делении (сдвиге вправо) в целочисленной арифметике теряются разряды. Поэтому промежуточный результат сохраняется в переменной z в масштабе суммы (до деления), а точнее, в масштабе выходной величины, деленной на Beta. Кстати, здесь и ответ на вопрос о необходимом диапазоне представления z: по максимальному значению выходной величины, деленному на Beta. Например, однобайтные переменные и коэффициент Nb, равный 1 (при последующем сдвиге на 8 бит, как в показанной функции) требуют двухбайтного z.

Еще один фокус — и я откланяюсь :)
Если не стремиться к конкретной цифре погрешности установления (0,1% или 5% или что-то еще), то можно попытаться выбросить одно умножение, выбрав меньшее из чисел Na и Nb равным 1.
В нашем примере меньшее из чисел именно Nb и легко посчитать, что оно может быть в диапазоне 8...17. Вычислим следующие границы:

256 / Nb = 32...16

Здесь числа равны степени двойки случайно. Важно то, что из диапазона допустимых значений 256 / Nb мы выбираем одно из чисел, равное именно степени двойки. Например, в нашем случае, берем 16.
Тогда 16 - это сумма «удобных» целочисленных Na и Nb, а меньшее из них равно 1, то есть микроконтроллеру не нужно умножать:

Y(n) = (15*Y(n-1) + X(n)) >> 4                                                                    (5)

Видим, что расчет (5) имеет на одну операцию умножения меньше, чем расчет (4).
 

РЕЗЮМЕ

Расчет простого цифрового фильтра нижних частот сводится к таким действиям:

Первое. Зададимся временем переходного процесса, периодом выборок и вычислим, за сколько выборок процесс должен закончиться (1): 

N = Тпп / Т

Второе. Вычислим границы возможных значений коэффициента при Y (2):

EXP(-6,908 / N) <= Alfa <= EXP(-2,996 / N)

Третье. Выбираем k (рекомендую 8, но можно и иное) и умножаем границы Alfa на 2^k. Выбираем любое из возможных значений в качестве первого коэффициента Na, а второй коэффициент находим как дополнение:

 Nb = 2^k — Na

Результирующая формула обработки входной последовательности X(n):

Y(n) = (Na*Y(n-1) + Nb*X(n)) >> k

Как показано в примере, иногда можно подобрать Na или Nb, равными 1, что упрощает вычисления.

БЛАГОДАРНОСТИ
Хочу выразить искреннюю благодарность уважаемому коллеге angel5a за то, что он вернулся из детского садика (кто читал дискуссию, тот поймет) и очень помог мне с переделкой исходной заметки. Радикальной переделкой, хочу заметить.
Также в доработке статьи помогли дискуссии с коллегами _pv и известным в наших кругах товарищем avreal, выступившим на ненашем форуме.

  • +18
  • 21 августа 2012, 18:03
  • drvlas

Комментарии (128)

RSS свернуть / развернуть
Первую статью не читал. Эту с первого прочтения понял. Досконально не разобрался, но я это сделаю при первой же необходимости. Ваша цель достигнута, спасибо )
0
drvlas, Спасибо за подробное толкование первого варианта.
Какую литературу посоветуете почитать человеку, несведущему в данной области, но желающему понять основные принципы?
(хотя Юкио Сато я прочел)
0
Так я ж… Эта… Не читатель. Я писатель. Современных книг не знаю. Учился лет 30 назад. ЩАС, погодь, коллеги обратят внимание, подскажут.
А вообще, основные принципы я понял не в институте (нам такого и не рассказывали), а только тогда, когда понадобилось практически. Брал книги, читал, зная, зачем мне это нужно. Тогда и обучение шло быстро. А иначе… Ой, как там всего много!
0
вне в прошлой статье посоветовали это:
www.eltech.spb.ru/techinfo.html?gid=13
0
:) И у меня в обсуждении это было. Я скачал и русский, и аглицкий варианты.
0
годная заметка
0
Не хватает структурной схемы цифрового фильтра для наглядности. А так, вроде, годная статься для начинающих — что может быть проще звена первого порядка :)
Из литературы для серьезной работы с цифрой могу порекомендовать книжки Сергиенко, к примеру, вот эту: http://www.ozon.ru/context/detail/id/5610414/. Можно еще поискать в интернете материалы моего учителя Глинченко А.С., по которым я, собственно, и учился. Его книжку мы в шутку называли библией :)
Ну а для самых начинающих надо, конечно же, брать Баскакова и Гоноровского :)
0
структурная схема (где плючики в кружочках) это уже для углубляющихся. А для начинающиц «черный ящик», где происходит вся магия, самое то.
0
Так и формулы тогда можно было не писать =)
0
А по чем тогда считать? :)
0
Так готовое же можно взять ;) Из какого-нибудь скрипта задать АЧХ сгенерировать сишный код =)
0
Пригодится, спасибо. В дебри не залезли, вроде понятно всё.
Пока у меня ума хватало только на усреднение по 2^n отсчётам.
+1
  • avatar
  • ploop
  • 22 августа 2012, 08:27
Спасибо за спасибо :) А еще за то, что (внезапно) навеяло мысль: можно было бы расширить в заметке сравнение фильтрации и усреднения. Понимая при этом, что при фильтрации необходимо выдавать данные на выход с частотой дискретизации (частотой поступления входных отсчетов). А усреднение допускает выдачу одного отсчета на N входных. И невзвешенное усреднение представляет собой простейшую реализацию КИХ-фильтра, который нагло пропускает N-1 отсчетов и выдает только N-й.
Так вот, популярность усреднения связана именно с допустимостью во многих задачах получать один результат за N*T. Как только заходит речь о настоящей фильтрации, тут же выбор между рекурсивным (обычно БИХ) и нерекурсивным (КИХ) фильтрами становится уже не столь очевидным. Тут уж нужно считать количество операций, которые МК придется выполнить для обработки очередного отсчета. А вот по этому параметру БИХ-фильтры обычно и выигрывают.
Поэтому, резюмирую, интерес к рекурсивным фильтрам резко повышается, когда темп выдачи отсчетов на выход высок. Правда, должен заметить, что такие задачи уже решаются более подкованными товарищами и мои заметки под девизом «цифровой фильтр — это просто» для них не нужны.
Но сделать как-то заметку о простом рекурсивном фильтре 2-го порядка — хочется. Ведь я глубоко убежден, что инженерам нужны простые подходы к расчету устройств, покрывающие 90% задач. Заодно и сам бы лучше понял. А то в девайсах у меня работают именно БИХи 2-го порядка, а рассчитывал их, стыдно признаться, моделированием в таблицах :)
0
У меня что-то зачесались руки проверить это на практике. Благо на чём — тоже есть. Стоят часы с аналоговым датчиком атмосферного давления. Разумеется шумит. Дабы сделать задержку обновления показаний ~20 секунд и избавиться от шума, сделал по фону усреднение по 65536 отсчётам (мощно, да? :) ). Цель достигнута — показания как вкопанные. Но чтобы проверить работоспособность пришлось дуть в трубку датчика все эти 20 секунд (я аж вспотел).
Так что для опытов девайс как нельзя лучше подходит. Проект на ассемблере, но тут проблем не будет.
0
Мощно да. Но это действительно просится на изменение.
Так 20 секунд и 64К отсчетов — это 300 мкс на отсчет. Так поступают отсчеты? И насколько ты проверял необходимость 64К? Например, есть ли разница при усреднении, скажем, 1К, 4К?
Вообще, желательно представлять, с чем же мы боремся. Если на входе, скажем, есть явно ВЧ-шум (буквально от отсчета к отсчету болтанка), то его можно придавить именно усреднением, причем невзвешенным и желательно кратным 20 мс. У тебя это примерно 67 отсчетов, если я правильно понял. А если полученные «мегаотсчеты» с частотой 50 Гц дают некоторую остаточную нестабильность — можно попробовать увидеть ее спектр. На худой конец — просто запустить мегаотсчеты в БИХ-фильтр. И сделать его переходной процесс в пару секунд.
А вообще, заведи тему на форуме, там поговорим. И там народ топчется, может получиться хороший топик.
0
Тема, думаю, ни к чему. Сейчас объясню, почему так поступил:
В программе используется много автоматов на прерываниях. Отсчёт времени, программный декодер пульта ДУ, опрос 1-wire датчиков и прочее. Помимо этого медленно ползёт основной цикл, в котором выводятся данные на индикатор, забираются данные с АЦП и он пинается дальше. Данные (16 бит) суммируются с буфером (32 бита). После переполнения 16-битной переменной-счётчика из буфера забираем 2 старших байта как результат.

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

По сути 64К выполняет тупо функцию задержки, «естественной» задержки, без дополнительных условий. По сути задача решена.

Что мне даст «скользящее усреднение»:
1. Понимание принципов работы (практика это не теория)
2. Наглядную демонстрацию работы, где можно «покрутить» к-ты и видеть результат в реальном времени.
0
честное скользящее среднее потребует буфер на все 65536 в данном случае отсчётов.
а этот простой бих фильтр — нет, при почти такой же АЧХ (те же 20 дБ/дек, только провалов нет на частоте кратной длине буфера), ну и ФЧХ немного другая, но кого она тут волнует.
если заменить усреднение на
y += ((adc << 16) — y) >> N;
и брать те же старшие два байта от y, то при N=16 получится абсолютно то же самое по шумам, только значение будет обновляться с каждым отсчётом, а не раз в 20 сек.
0
По мне, например, лучше и нагляднее рассчитывать фильтры по АЧХ, а не по гипотетическим отсчетам… А для этого матлаб тема, или WinFilter проще уж некуда для инженера…
0
Ну, пусть расцветают 100 кустов роз. Кому-то подойдет и временнОй подход. Я постарался объяснить, почему это может быть удобно. И именно для фильтра 1-го порядка. Попробуй описать, как ты будешь формулировать ТЗ по АЧХ. Не вообще сферического фильтра в вакууме, а вот такого, «убирающего болтанку». Придешь к выводу о необходимости исследования спектра помехи, спектра полезного сигнала :) Это не подход для новичка. А я даю именно для него, родимого
+1
Присоеденяюсь. Статья проста и понятна. Спасибо.
0
Внес последние правки, уважаемый соавтор! Спасибо!
0
Только мне кажется, что проще на вход АЦП поставить R/C цепочку, а соотношение R и C дадут те же альфу с бэттой. Я даже примерно понимаю, как их рассчитать…
0
  • avatar
  • dekar
  • 22 августа 2012, 19:45
Я даже примерно знаю, что ты не читал длинннннючие дискуссии к первому варианту статьи :) Там я приводил примеры, когда перед цифровым фильтром вовсе не существует аналогового сигнала: если мы запускаем не его вход уже продукт цифровой обработки. И здесь чуть выше в разговоре с ploop мы говорили о варианте обработки предварительно усредненного в цифре сигнала.
Ну и параметры RC в некоторых случаях могут быть ужасны.
Так что ставь цепочку. Но другой вариант имеет право быть.
0
ок, но если бы я писал статью (смешно, сам знаю), то писал бы про линейный фильтр как аналог RC. А тут в статье вообще не слова.
0
когда перед цифровым фильтром вовсе не существует аналогового сигнала: если мы запускаем не его вход уже продукт цифровой обработки.
Да пусть даже и оцифровываем — собственные шумы АЦП, нехорошее питание, неправильную разводку — да кучу факторов можно придумать, когда RC не спасёт.
Один раз по неопытности развёл так, что на выходе получил шум на 4 разряда, который наблюдался при посаженном входе АЦП на землю. Какой уж тут RC…
+1
Бывает и того хуже, например аллегровские датчики тока имеют на выходе (в зависимости от типа) 7-20мВ шумов (при полной амплитуде сигнала 2.5В). Тут уж как не разводи…
0
Для оценки результата можно взять конкретный сигнал (синус, к примеру), наложить на него известный шум (с равномерным распределением, к примеру) и прогнать через этот фильтр. Поскольку модель сигнала известна, то шум на выходе считаем простым вычитанием, учитывая смещение и коэффициент усиления. Саму оценку количественно можно вычислять при помощи соотношения сигнал/шум до и после преобразования.

А то, вроде фильтр, вроде фильтрует то, что «скачет», а где сама оценка фильтрации того что «скачет»? И как «скачет»? Рысью или галопом? А что будет на выходе при одиночном выбросе типа «огого какой выброс»?
0
  • avatar
  • uni
  • 22 августа 2012, 23:48
Немного не понял. Это где, куда и почему? Что именно требуется оценить?
Если рассчитанный по этой методике фильтр — то много проще внести его к-ты в ту самую простую утилиту, которую мы теперь все знаем :) И она покажет к-т передачи на частоте полезного сигнала и на частоте любой помехи. Что будет на выходе при импульсной помехе — можно посчитать как утилиткой, так и экселевской моделькой.
Но смысл? Я прошу понять моих «клиентов»: они все эти модели не знают и не умеют строить. Они даже свои реальные шумы на входе не хотят или не могут оценить — ни по спектру, ни по уровню. Куда уж там говорить о верхней границе спектра полезного сигнала… Но у них показания после АЦП «скачут». Да, уважаемый, галопом, рысью или еще каким конем — но скачут. А хочется, чтобы не скакали. Я и говорю: хрен с ним, с ТЗ на проектирование кашерного фильтра! Давай так: можешь прикинуть на глазок, за сколько времени показания должны устаканиться после подачи на вход сигнала? Знаешь частоту входных выборок? ОК, вот тебе флаг формулка и алгоритм обработки. Сделай и посмотрим, мож поменьше будет «скакать». И уже чел не от балды лепит, а фильтр по формулам. Рекурсивный, йитить его налево! Наука, мля…
Вот за этого новичка я тут и воюю. А вы меня все пытаетесь поставить на более строгий путь.
Со всем уважением,
+1
Если рассчитанный по этой методике фильтр — то много проще внести его к-ты в ту самую простую утилиту, которую мы теперь все знаем :) И она покажет к-т передачи на частоте полезного сигнала и на частоте любой помехи. Что будет на выходе при импульсной помехе — можно посчитать как утилиткой, так и экселевской моделькой.

А причём тут «частота любой помехи»? Равномерный шум — это шум на 0 до бесконечности. Как правило, всегда есть таковой. Клиенты, я так понял, хотели бы оценить насколько таковой фильтр является фильтром. Если «отбросить» всю науку вообще (чего я очень не рекомендую делать), то самое простое — не использовать никаких спектров. Я же вроде понятно написал, что оценку делаем во временной области: «Поскольку модель сигнала известна, то шум на выходе считаем простым вычитанием, учитывая смещение и коэффициент усиления».

Что тут сложного-то? Совсем в ясли-то не будем опускаться. Есть сигнал во временной области. Плюсуем к нему шум какого-то уровня и типа и подаём на вход блока цифровой обработки. На выходе по идее должно быть шума меньше. Можно это оценить и на глазок, но, боюсь, это вряд ли кого удовлетворит. Гораздо полезнее был бы вывод: было отношение мощности сигнала к мощности такого-то шума во временной области такое вот (цифра), а после обработки отношение мощности сигнала к мощности шума стало вот такое вот (цифра). Все радуемся от того, что вторая цифра больше первой. Но если вторая меньше первой или равна ей, то думаем почему так.
0
По поводу «выброса». Если после фильтра стоит решающее устройство и оно должно реагировать на какие-то пороги сигнала, то что будет, если сигнал один раз стал равен нулю? Мы ведь также предполагаем, что сигнал может быть чисто цифровой с датчика температуры к примеру. Температура измеряется через SPI (или I2C) и подаётся на фильтр, а дальше пороговое устройство. И вот незадача, периодически связь куда-то теряется или какая ошибка в логике работы по SPI в результате чего не принятое значение считается нулём. Может такое быть? А почему нет? Что будет на выходе?

Фильтрованием чего занимается такой фильтр? Не плохо бы показать на примерах и количественно «пощупать» действие фильтра.
0
Не пойму, почему Вы так упорно пытаетесь навязать автору свой взгляд на то, как нужно излагать материал. Напишите свою статью. Хорошая статья о цифровых фильтрах никогда не будет лишней.
+3
А почему он мой, собственно, и когда стал? Я попросил оценить во временной области саму работу фильтра. Это разве не само собой разумеющееся?

Есть блок цифровой обработки, неважно какой, есть набор каких-то входных воздействий (у каждого свой, но не у всех ступенька). Разве не нужно после проектирования (расчётов) показать работу блока на чём-то реальном? Или сферические кони в вакууме нынче в моде?

Вот кто-то захочет себе такой фильтр в программу и у него вопрос возникнет: вот у меня сигнал на входе, а что ожидать на выходе у такого фильтра (количественно)?

Или у всех сигнал в виде ступеньки и этим вся цифровая область их разумения этим и ограничена? Ну тогда извиняйте :)

Читаем ещё раз:

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

Автор поставил задачу. А где та мера, которой измеряется степень уменьшения? Каких шумов? Шум бывает узкополосный, широкополосный. Во временной области широкополосный шум — это могут быть как импульсные помехи, так и просто «белый шум». Могут быть одиночные выбросы, как я описал выше, связанные с проблемами при передаче цифры. Вот не сгладит, к примеру, такой фильтр одиночный выброс и сработает ваше устройство от помехи при передаче, а управлять оно может чем угодно.
Как узнать на что рассчитывать? На обтекаемую фразу «влияние шумов было уменьшено»? Так каждый это понимает по-своему.
0
Оставим примеры с картошкой, я с утра не евши :)
Если я правильно понял, Вы предлагаете добавить в обсуждаемую статью следующее: оценку качества полученного фильтра. Верно?
Как (со)автор статьи, я не вижу в этом смысла. Что мы предлагаем в статье? Вернеее, что мы НЕ ПРЕДЛАГАЕМ? Мы не предлагаем фильтр с новым, неизвестным ранее способом обработки сигнала. Наоборот, все, что происходит с сигналом, давным-давно расписано до тошноты, до рвоты. Тогда какой пруф еще нужен? Специалист Вашего (и даже много более низкого) уровня леХко, влет посчитает все результатирующие параметры, леХко найдет ту самую меру улучшения сигнала, леХко добавит на выходе этого фильтра дополнительную решалку проблем с выбросом и т.д.
А мой клиент вообще заморачиваться не будет. Он посмотрит на болтанку после фильтра и скажет: «Спасибо сайту EASYelectronics, полегчало». Или напишет мне недовольное письмо, мол, херня все ваши примочки, у меня вход АЦП на ноль замкнут, а болтанка все равно слишком большая.
Так что — не вижу необходимости усложнять статью.
+1
Вот кто-то захочет себе такой фильтр в программу и у него вопрос возникнет: вот у меня сигнал на входе, а что ожидать на выходе у такого фильтра (количественно)?
20дБ/дек начиная с примерно Fs/K, АЧХ как у RC цепочки.
собственно формула Y[i+1] = Y[i] + (Y[i]-X[i])/K, поведение интегрирующей RC цепочки и описывает. Uc = \int{Ui-Uc}/RC
Или у всех сигнал в виде ступеньки и этим вся цифровая область их разумения этим и ограничена? Ну тогда извиняйте :)
насколько помню теорию линейных электронных схем любой фильтр полностью и однозначно характеризуется реакцией на ступеньку либо на одиночный испульс.
0
Простые примеры. Вы когда машину себе покупаете, то подход какой? Чтобы ездила? А квартиру? Чтобы жить можно было? А жену? Чтобы рожать могла? А комп? Чтобы Линукс устанавливался? Даже когда вы покупаете картошку, то вряд ли вас интересует только то, что на бирке написано «картофель», как минимум нужно, чтобы он был определённого размера, не гнилой, не мороженный, не зелёный, мытый (без земли) и одинаковой формы, пригодной для ручной чистки, так ведь?

А теперь сравните: «наша задача — сгладить входную последовательность, то есть, выдать выходную последовательность такую, чтобы влияние шумов было уменьшено».

Это всё равно что ребёнку мелкому дать указание сходить в магазин и купить картошку 1 кг. Он возвращается и приносит, а вы ему: «определённого размера, не гнилой, не мороженный, не зелёный, мытый (без земли) и одинаковой формы, пригодной для ручной чистки». Так ведь? А он и отвечает:

Мама, не пойму, почему ты так упорно пытаешься навязать мне свой взгляд на то, как нужно покупать картошку? Там было написано: картошка, 1 кг — вот и получи.
0
прошу простить, что вмешиваюсь, но спор давно зашел в тупик, как мне кажется.
Вы «проповедуете» правильный инженерный подход (или зачатки оного) и никто, вроде, и не пытается оспаривать всю правильность изложенных соображений по поводу того, каков путь истинный как все должно быть…
Но ведь ув. автор уже не раз ответил, что ставил перед собой иные цели, состоящие в том, чтобы показать неподготовленным новичкам наипростейший способ решения некоторых насущных проблем, при пом. цифрового фильтра. И часть из этих новичков захочет, естественно, разобраться более подробно, начнет читать литературу или другие, более близкие к теоретическим основам статьи, чтобы понять общий подход и научиться использовать более правильные методы проектирования. А кто-то просто опробует эту методу и решит что ему ее достаточно для своих конкретных задач. И зачем человеку, который хочет сделать лично для себя «умную табуретку» или регулятор оборотов вентилятора в сортире изучать матапарат, читать литературу по ЦОС и углубляться в теорию ТАУ?
Или вот такая «аллегория»:
При условии, что вы хотите помочь утопающему человеку, что рациональнее: вытащить его из воды, а потом обьяснять ему как важно научиться плавать самостоятельно, или же, давая советы и показывая движения учить его плавать, чтобы он самостоятельно выбрался на сушу? (а если у человека «горит» проект?).
+1
Спасибо за понимание, коллега! Думаю, что мы с уважаемым uni уже поняли друг друга.
В развитие темы — мечтаю сделать еще парочку заметок с таким же подходом. Думаю, что буду концентрироваться на фильтрации и усреднении. Когда-то у меня были очень неплохие практические результаты. Ну, вспомнить все… непросто. Однако это больше вопрос времени и вдохновения.
0
спасибо за статью.
В развитие темы неплохо было бы про медианные фильтры рассказать.
Весьма эффективные от разных «выбросов»
0
Ну, медианный фильтр штука хорошая. Но я лично его не применял. ЩАС специально полез в Википедию, чтобы уточнить, как там работает. Да, это очень похоже на КИХ-фильтр. И должно намертво отметать определенную долю шумов (короткие выбросы, искажающие входной сигнал на длительность, более короткую, чем импульсная характеристика). Занимательно!
Но я люблю рассказывать только о том, что прошло через собственные ручки :)
Можешь написать, а мы почитаем :)
0
Медианный фильтр очень хорошо справляется с «выбросами» (особенно выигрышно но смотрится если выбросы случайны и с большой амплитудой). Но его реализация требует очень много ОЗУ (притом затраты растут нелинейно в зависимости от разрядности АЦП). Я, честно говоря, очень редко встречал медианный фильтр в проектах на МК.
0
Это похоже подразумевает использование radix sort с одним разрядом, ради эффективности. И лучше ничего не придумаешь.
0
Гм. Нинай про какой ты, но вроде медианный фильтр для изображений особо много ОЗУ не требует — по сути только под несколько отсчетов.
Правда расчет веселый — эти отсчеты нужно отсортировать и выдать средний.
0
Но его реализация требует очень много ОЗУ (притом затраты растут нелинейно в зависимости от разрядности АЦП).
Извиняюсь, коллеги. Вспылил, был не прав :) Я спутал медианный фильтр с другим алгоритмом фильтрации.
0
Да неплохо было бы расскзать. Мне его работа что-то не кажется очевидной. И ещё более неочевидна эффективная численная реализация.
0
Принцип-то, практически очевиден (в процессе сортировки выбросы в обе стороны уходят в начало и конец списка, таким образом в средине остаются только «хорошие» элементы), а вот эффективная реализация это да, не просто…
0
То, что выбросы уходят это да, ясно. А вот, что остается, не очень. Достаточно представить случай когда выбросов нет, чтобы начали возникать вопросы.
0
Ну, не так все сложно, ИМХО. Этот фильтр срезает вершины. То есть, пока в пределах его ИХ входной сигнал монотонен (в т.ч. и постоянен) — фильтр тупо вносит задержку на половину своей ИХ. Как только начинаем переаолзать через вершину — срезает ее, в зависимости от соотношения длины ИХ и остроты вершины.
Вот простая моделька в таблицах. Покзано просто синус — и его фильтрация. Видно, что FILTR_SIN имеет маленькие площадочки на вершинах. Потом добавил 5-ю гармонику в синус. И все равно видно, что монотонные участки FILTR_GARM простоотстают от входного сигнала GARM. А вершинки срезает.
Ну, нелинейный, что возмешь...
0
Я думаю, что тут в игру вступает распределение отсчетов с учетом шума. Наиболее вероятные значения будут чаще других попадать в средину списка, откуда потом они попадут на выход фильтра. Впрочем, это чисто умозрительное заключение, возможно я ошибаюсь.
0
«на выходе получил шум на 4 разряда, который наблюдался при посаженном входе АЦП на землю»
Делать на это цифровой фильтр-как говорит мой знакомый- лепить из г… пулю. Не обижайтесь! Какая же там разводка?(при посаженном входе АЦП на землю).
0
  • avatar
  • svs39
  • 22 августа 2012, 23:53
Да ладно, получить такую фигню — запросто. Аналоговую землю с цифровой по цепочке соединить и все…
0
Давно это было, не помню точно, но помню, что АЦП (или опора) сидел на сильноточной дороге, где ШИМ молотил.
0
Не, я конечно понял в чём дело, и распаял эту хрень сразу, переразвёл, кондёров добавил, всё по правилам. Но в памяти засело намертво. Ещё у меня там МК ресетился сам по себе :)
0
спасибо за статью! проверил на аналоговом датчике, хорошие результаты и легко настраиваемые результаты! до этого использовал RC фильтр и усреднение по 8 отчетам.
0
еще хочу немного докопаться — хотелось бы формулы видеть не «курсив + италик» а в нормально, «формульном» виде (пусть и в виде картинок)
0
Понятия не имею, как вставлять сюда формулы! Про картинки не думал, но немного вломно :) Может кто-то из товарищей подскажет, можно ли формулы более красиво средствами HTML изобразить.
А за доброе слово — спасибо! Ради этого все делалось — чтобы «хорошие и легко настраиваемые результаты».
0
Ну, смотрите, одну формулу (2) переписал! Чуть попозже переделаю и другие. ЩАС дела…
0
Что-то показатель степени шрифтом крупнее основания. Со школы ещё привычно наоборот.
0
Да и знаки неравенства великоваты. Думаю, нужно как-то увеличить сами Алфьа и Е. Но не могу сейчас.
0
«курсив + италик»
Гм. Это, как бы, одно и то же)
Вообще, формулы такого уровня можно прилично форматировать при помощи HTML. Но в новой версии LiveStreet почему-то такие полезные вещи, как <sup> и <sub> отключены в css.
Алсо, alfa и beta вполне можно заменить на α и β (&аlpha; и &bеta;).
0
sup — это верхний индекс, sub — нижний? Так Бог с ним, с тем LiveStreet, я не знаю о чем ты. Значит ли это, что сейчас я могу войти в редактирование файла, выложенного в Сообществе, и ручками сделать формулу?
0
sup — это верхний индекс, sub — нижний?
Да. Стандартные HTML-тэги. Тока тут они не работают — для них в css выставлено «отображать как обычный текст».
LiveStreet — движок, на котором это сообщество работает. До обновления его sub и sup показывались нормально… Зато сам движок регулярно падал.
Значит ли это, что сейчас я могу войти в редактирование файла, выложенного в Сообществе, и ручками сделать формулу?
Да, в тех пределах, в каких позволят HTML и местный движок.
0
ой, курсив + полужирный :)
0
Теория хорошо, я думаю если бы добавить код на С чтобы можного скопилировать и посмотреть на практике — было бы супер!!!
0
  • avatar
  • Nemo
  • 24 августа 2012, 11:22
Код-то здесь простенький получается:
#define Na 15
#define Nb 1
#define k 4

int filter(int x){
  static int y1 = 0;
  return y1 = (Na * y1 + Nb * x) >> k;
};

Вот только какие данные ему скармливать и что с выходными данными делать? И под какую платформу вообще делать?

Хотя можно в маткаде просимулировать. Сгенерить «данные» (синус скажем), зашумить их рандомом, прогнать через фильтр и нарисовать графики до и после.
0
Если интересно просто посмотреть работу такого фильтра, то лучше Экселя ничего и не придумаешь. Могу, конечно, и сам сделать — специально для коллеги Nemo. Но думаю, что это лучше делать под свои задачи. Например, каким шумом задаться? Какой полезный сигнал интересует? Ступенька, импульс, гармонический? Масса вариантов. Я сделаю то, что мне кажется наглядным и актуальным, а у Nemo задачи в другой области…
0
как только будет время, хочу залить код в микроконтроллер AVR или Cypress. Сделать цепочку ADC — uC — DAC. На вход подать сигнал с помехами а на выходе DAC посмотреть что получилось. Спасибо!
0
Если интересно просто посмотреть работу такого фильтра, то лучше Экселя ничего и не придумаешь.
Потратьте немного времени на изучение Matlab (или его свободного аналога SciLab)
Уверен, не пожалеете!
После него анализ сигналов в экселе — как разводка плат в паинте.
Тем более, что ЦОС там представлен в виде отдельного раздела с соответствующими библиотеками.
Можно перепробовать все известные науке велосипеды, прежде чем изобрести свой.
+1
тогда уж:
int filter(int x){
  static int y = 0;
  return y += (x - y) >> k;
};
0
Побойся Бога! Формула неверна. И потерял Na и Nb. Правильно у Vga — см выше
0
и правда, чего это я тут ересь несу какую-то: 15/16 это ведь совсем не тоже самое что 1 — 1/16
0
Все правильно. просто это оптимизированный вариант. Подходит для случаев, когда Nb=1.
0
Ребята, вы правы. Пардон!
0
по смыслу такая запись куда более похожа на описание RC цепочки.
Uc = \integral{(Ui — Uc) / RC}
приращение значения на выходе (зарядный ток) пропорционально разнице значений между входом и выходом.
для лучшего понимания я бы так описывал:
Y[i] = Y[i-1]*(1-K) + X[i-1]*K <=> Y[i] = Y[i-1] + (X[i-1] — Y[i-1]) * K;
ну и пару тактов на умножениях экономит.
0
Разница напряжений — это вход минус выход. И она, деленная на R, дает ток, который интегрируется с к-том 1/С и ПРИБАВЛЯЕТСЯ к имеющемуся. Все верно. И коэффициент, с которым интеграл от разности напряжений прибавляется к выходному — как раз 1/RC. Так что аналогия есть.
Но насчет экономии — да, и она есть. Покажу, как это выглядит применительно к моим расчетным значениям Alfa и Beta.
Перепишем формулу, по которой ты предлагаешь работать:
Y(n) = Y(n-1) + Beta*X(n) — (1-Alfa)*Y(n-1)
или, с учетом того, что (1-Alfa) = Beta:
Y(n) = Y(n-1) + Beta*(X(n) — Y(n-1))
Таким образом, только при значениях Beta, равных 1/2^k, мы получим вообще без умножения — сдвигом (X(n) — Y(n-1)) на k разрядов. А если Beta «кривое», то умножаем на Nb и делим на 1/2^k сдвигом.
Это нужно внести в статью. Спасибо!
0
Забавный вариант.
0
При x — y < 2^k у нас на выходе не будет изменений. т.е. будет не ассимптота к x уже. Не забываем про особенности целочисленной арифметики.
+1
Мое почтение, уважаемый соавтор! Но ты не прав, ИМХО. По разрядности я промоделировал (потому что тоже показалось стремно) и увидел, что разницы между записью
Y(n) = (Na*Y(n-1) + Nb*X(n) >> k
и записью
Y(n) = Y(n-1) + Nb*(X(n) — Y(n-1)) >> k
нет. Я взял эксель
ROUNDUP((1000+С1*15)/16)
и
ROUNDUP(C1+(1000-C1)/16)
где 1000 — входной сигнал Х, С1 — предыдущее значение Y.
Оба варианта сходились к 1000 на 73-м шаге — и ни одного шага не было с разницей хоть в 1.
0
округлять не результат надо, а промежуточные вычисления
A2: =A1+ROUND((100-A1)/16;0) и заполнить вниз
0
возражение принимается:
int filter(int x){
  static int y = 0;
  return (y += ((x << k) - y) >> k) >> k;
}

правда теперь надо за максимальными значениями x следить чтоб x<<k в int влезало.
0
Теперь сдвигов дофига. А на некоторых МК это не такая уж дешевая операция (например x51 и AVR — сдвиг доступен только на один бит, причем у первого есть аппаратное умножение и деление).
Да и менее понятен такой код.
0
return (y += ((x << l) — y) >> k) >> l;
l можно кратным 8 сделать, тогда компилятор (хороший компилятор) впринципе должен разрулить это просто подвинув адреса на единичку у команд ld, st.

а с кодом
int filter(int x){
  static int y1 = 0;
  return y1 = (Na * y1 + Nb * x) >> k;
};

та же самая проблема, x округляется до кратного (1 << k) / Nb.
поэтому перед фильтром его надо бы в старшие разряды подвинуть на 8/16бит, чтобы сам фильтр работал с большими числами и не боялся округлений, а потом выход фильтра на те же 8/16 бит подвинуть обратно (взять только старшие байты)
0
Я выше ошибся с функцией округления (механически щелкнул по выбору ROUNDUP вместо ROUNDDOWN). Поставив округление ROUNDDOWN, которое автоматически получается в целочисленной арифметике, получил, что в обеих записях процесс не доходит до конца. Застревает на 985…
Но решение есть. Если на каждом шаге просто сохранять новое значение Y до его деления на 2^k:
int filter(int x){
  static int y = 0, Y=0;
  Y += x + y;
  return y = (Nb * y) >> k;
};

Тогда уже на 117-м шаге процесс доходит до 1000.
То есть, мы не выплескиваем с водой ребенка: выход приводим к масштабу Х (как и положено), пусть и с потерей дробной части но неделенный выход помним и используем в следующем цикле
0
там пара опечаток в коде.

int filter(int x){
static int y = 0, z = 0;
z += (x — y);
return y = z >> k;
};
0
Ну да, x-y. Я тут в экселе проверяю, а не в Си. Но умножение на Nb надо добавить. Это более общий случай. ЩАС подумаю
И почему обязательно z? Си различает регистр?
0
человек на мониторе плохо различает однобуквенные переменные которые отличаются только регистром.
y = (Nb * Y) >> k;
0
Не, надеяться на кратность сдвигов восьми — не наш метод. Смотри ниже, есть выход универсальный
0
ОК, проверил. Да, если формула

Y(n) = Y(n-1) + Beta*(X(n) — Y(n-1))

как более общая, то функция такая:

int filter(int x, int Na, int k){
static int y = 0, z = 0;
z += (x — y);
return y = Na * z >> k;
};


Кстати, хорошо видно при моделировании, что процесс ускоряется с ростом Na, то есть к-то при разности. Если он 1, то успокаивается до 0,1 % (вот та самая 1000 на выходе при 1000 на входе) за 117 циклов. А если взять Na=3, то тот же результат за 37 тактов. Понятно, мы ведь уменьшили ТАУ в 3 раза — вот и переходной процесс ускорился в 3 раза.
Это еще одно преимущество твоей записи!
0
return y = Na * z >> k;
так смело без скобочек можно писать только если наизусть знаешь все приоритеты всех операций в С.
я не знаю и оно скорее всего у сдвига меньше чем у умножения, но вот так вот можно на ровном месте обратно получить округление.
поэтому лучше return y = (Na * z) >> k;
0
Да я всегда леплю скобок побольше. Тут спешил. ОК, учту. Переделаю статью на красивые формулы, поставлю твою формулу вместо своей, ибо лучше же, что и говорить! И покажу реализацию на Си. Но! Все это после поездки. ЩАС сборы, тудам-сюдам, не дадут сосредоточиться
0
Какой своей? Изначальной y = (Na * y + Nb * x) >> k? ИМХО, ее стоило бы оставить, но дополнить формулой _pv с комментариями на тему различий.
0
Найди мне хоть одно преимущество этой моей формулы — и я ее оставлю. Видно же, что умножение в принципе одно (или ноль, но это частный случай), и лишь на 1 вычитание больше в формуле от _pv. Или я опять не все понимаю?
0
Она очевидна и понятна. И на ней можно показать типичные проблемы — причину перехода к менее понятной формуле. Если ее выкинуть — повествование станет непоследовательным.
0
Может вместо Na все-таки Nb?

int filter(int x, int Nb, int k){
static int y = 0, z = 0;
z += (x — y);
return y = (Nb * z) >> k;
};
0
Дык… Исправил в статье еще вчера. Там новый вариант функции лежит. Но все равно мерси, во входных переменных оставалось Na. Поправил :)
0
Там у тебя еще какой-то «ucha». Никогда не слышал о таком типа денных.
0
  • avatar
  • Vga
  • 28 сентября 2012, 16:29
Тьфу, черт. «типе данных».
0
  • avatar
  • Vga
  • 28 сентября 2012, 16:30
Тьфу, черт :)

typedef unsigned char ucha;
0
Хоть бы uchar его назвал. Но лучше поменять на uint8_t и int16_t. int и char могут быть совершенно произвольного размера.
0
  • avatar
  • Vga
  • 28 сентября 2012, 16:39
Дело привычки. У меня обычно не публичные проги, так что за 20 лет что-то выработалось (пусть убогое) — и меня не жмет.
Но вообще, ты прав. Теперь, когда приходится людЯм показывать, надо думать о едином алфавите :)
Кстати, верно ли, что
typedef unsigned char ucha;

лучше, чем
#define ucha unsigned char

Мне кажется, что первое дает компилятору больше инфо, чем просто подстановка вместо ucha. И это полезно. Или нет?
0
Мне кажется, что первое дает компилятору больше инфо, чем просто подстановка вместо ucha. И это полезно. Или нет?
Понятия не имею.

Но типы в примере все же лучше поменять на имеющие фиксированную разрядность.
0
  • avatar
  • Vga
  • 28 сентября 2012, 16:51
Удивлен. Мне казалось, что ты просекаешь многие моменты, недоступные мне :)

Ну, ладно. А типы в примере… Потом, как-нибудь. Когда табличку сделаю и сяду править.
0
Погодь, я что-то не понял. Расчет Na и Nb, которому, СОБСНО, была посвящена статья, делается на абакусе калькуляторе. А исполняемый алгоритм описан вот такой формулой:

Y(n) = (Na*Y(n-1) + Nb*X(n)) >> k

Она, практически, на Си :)
Мне кажется, что пример программы на МК только затемнил бы изложение.
0
Вот так примерно оно выглядит в маткаде:
+2
Спасибо! Именно это я имел в виду: результаты предсказуемы, увидеть их можно в любом пакете, который позволяет строить графики.
То есть, мне с Экселем можно не напрягаться :)
0
Хочу поделиться небольшим опытом. Я для фильтрации взял себе на заметку «Practical Analog and digital filter design». Автор Les Thede. Тупо из книжки вбил в максиму 8 алгоритмов для Баттерворда и Чебышева 1 рода соответственно для низких частот, для высоких, полосовой и режекторный. И от плавающего усреднения в чистом виде отказался.
0
Книжка хороша! Спасибо, скачал.
Однако же, что значит «тупо вбил из книжки»? Написал программу расчета? Так — просим, просим! Или воспользовался прилагаемым диском с прогой WFilter?
То есть, я голосую за максимально простые подходы к расчету, дающие решения новичкам в подавляющем числе случаев (во всех, фактически), а обычным инженерам — просто в большинстве. Если из этой книги, скорее всего дающей именно такие инструменты, ты можешь предложить альтернативу скользючке, то было бы очень хорошо увидеть это в виде статьи. Или хотя бы здесь, в комментариях. Дело в том, что самое интересное — увидеть альтернативный метод именно для такой задачи, как я сформулировал, а не для сферического коня в вакууме.
+1
Давайте, я пример приведу.
В максиму вбил 5 формул из соответствующего раздела для нужного фильтра. Пример ниже для режекторного фильтра для вырезания частоты розетки. Скрипт даёт целочисленные коэффициенты для задаваемой разрядности (32, 16, 8 бит) ну чтобы и на AVR и на STM32 работало с естественной для МК разрядностью. Ну и графики рисует. Собственно, 70% кода для рисования графиков амплитуды и фазы.
Скрипт, приведённый ниже даёт результат для коэффициентов при X
0.97861561896753
-1.864775369303321
0.97861561896753;
а при Y
1.864775369303321
-0.95723123793506
23 битные целочисленные аналоги для X:
8209223
-15642870
8209223
и для Y
15642870
-8029838.

Параметры фильтра задаются в самом верху скрипта: частоты пропускания, частоты подавления, ну и во сколько раз задавить. Скрипт подбирает минимальный порядок фильтра и для заданной разрядности выдаёт коэффициенты (и с плавающей точкой и целочисленный для заданной разрядности). Наверное, мог бы и статью оформить, если люди посчитают это полезным.

<em>-->	 numer: true;
-->	 fpprec : 64;
Filter parameters.

-->	 aPass: -1;
-->	 aStop: -60;
-->	 wPass: 1;
-->	 wStop: 2;
Unnormalized paremeters.

-->	 fPass1: 43;
-->	 fPass2: 57;
-->	 fStop1: 49.8;
-->	 fStop2: 50.2;
-->	 Wcs: ( fPass2 - fPass1 ) / ( fStop2 - fStop1 );
Calculatable constants.

-->	 eps: sqrt( 10^( -0.1*aPass ) - 1 );
-->	 n: ceiling( acosh( sqrt( ( 10^(-0.1*aStop) - 1 )/( 10^(-0.1*aPass) - 1 ) ) ) / acosh( Wcs ) );
-->	 D: asinh( 1/eps )/n;
-->	 if ceiling( n/2 ) = n/2 then
  ( cnt: n/2, sigmaN: 1 )
else
  ( cnt: (n-1)/2, sigmaN: -sinh(D) );
-->	 Hs: 1;
-->	 for i from 0 thru (cnt-1) do
    ( fi: %pi*(2*i+1)/(2*n),
      sigma: -sinh(D) * sin(fi),
      w: cosh(D) * cos(fi), 
      B2m: sigma^2 + w^2, 
      B1m: -2 * sigma, 
      Hs: ratsimp( Hs * B2m / ( s^2 + B1m*s + B2m ) ) );
-->	 if ( ceiling( n/2 ) = n/2 ) then
    ( Hs: ratsimp( Hs * 10^(aPass/20) ) )
else
    ( Hs: ratsimp( Hs * sinh(D)/(s+sinh(D)) ) );
1 Normalized bandstop.

Parameters for normalized to bandstop convertion.

-->	 BW: (fPass2 - fPass1) * 2 * %pi;
-->	 w0: sqrt( fPass1 * fPass2 * (2 * %pi)^2 );
-->	 Ss: BW * s / (s^2 + w0^2);
-->	 notchHs: ratsimp( subst( Ss, s, Hs ) );
2 Frequency responce.

Frequency response.

-->	 complexHs: subst( 2*%pi*%i*s, s, notchHs );
-->	 absHs: cabs( complexHs );
-->	 phHs: carg( complexHs );
-->	 wxplot2d([absHs], [s,0,120]);
-->	 wxplot2d([phHs], [s,0,120]);
3 Discretizing.

-->	 f: 1000;
-->	 Hz: subst( 2*f*(1-z)/(1+z), s, notchHs );
-->	 Hz: ratsimp( expand( Hz ) );
-->	 numHz: num( Hz );
-->	 denHz: denom( Hz );
-->	 coeffX0: coeff( numHz, z, 0 );
-->	 coeffY0: coeff( denHz, z, 0 );
-->	 load( powers );
-->	 res: powers( numHz, z );
-->	 szX: res[1]+1;
-->	 res: powers( denHz, z );
-->	 szY: res[1]+1;
Divide all coefs on denominator's coef at power 0.

-->	 array( coefX, szX );
-->	 array( coefY, szY );
-->	 for i from 0 thru (szX-1) do
    ( coefX[i]: coeff( numHz, z, i )/coeffY0, 
      print( coefX[i] ) );
-->	 for i from 0 thru (szY-1) do
    ( coefY[i]: -coeff( denHz, z, i )/coeffY0, 
      print( coefY[i] ) );
Impulse responce.

-->	 kill( w );
-->	 absHz: cabs( subst( exp(-%i*w), z, Hz ) );
-->	 wxplot2d([absHz], [w,0,%pi]);
Created with wxMaxima.</em>


Амплитуда аналогового фильтра

Фаза аналогового фильтра.

Амплитуда цифрового фильтра для частоты дискретизации 1кГц
0
Да, извините, у меня код как-то криво вставился с лишними символами.
0
Ну, статья была бы интересной. А в приведенном примере мне мало что понятно. И самое главное: какова была бы процедура, которую наш читатель применил бы для решения задачи в формулировке, изложенной выше в разделе «Постановка задачи»? Оставь на секундочку те 8 решений, что ты реализовал по книге. Пусть они пойдут в статью с соответствующей темой. Но приведи то из них, которое успокаивало бы «болтанку» с неизвестным спектром шумов минимальными ресурсами МК и простейшей программой.
Вот это решение я бы и хотел сравнить с предложенным мной «честно простым фильтром». Не для определения победителя, а для того, чтобы люди могли использовать на практике лучшее. Или просто имели выбор на свой вкус.
0
Простите что я встрял, пожалуйста, не ищите во мне конкурента! Я совершенно не с целью побеждать или проигрывать! Просто несколько раз сталкивался с необходимостью усреднять, вырезать и выделять диапазон цифровым образом. Хотел поделиться соображениями, что и сделал.

Навскидку мне думается, для решения задачи подошёл бы IIR фильтр нижних частот Чебышева — в нём максимально быстро падает частотная характеристика. Для первого порядка он и будет выглядеть практически так же, как скользящее среднее. Но если нужен более быстрый спад, то порядок будет чуть больше (задействовано не только последнее значение с АЦП, но и предпоследнее и, возможно, второе с конца), зато усреднённый сигнал будет, чётко отделять целевые изменения среднего от шума, так как частота среза задаётся явно, см. картинку.

Фильтр нижних частот Чебышева

Но если предполагается, что значение фелевого сигнала почти постоянно или динамика процесса не важна, скользящего среднего больше чем достаточно! Извините что встрял!
0
Да что ты, дорогой, я же не наезжаю! Наоборот, интересно. Ведь мы же здесь не показываем свои методы — а лишь приложение известных. Тут авторское эго не ущемляется.
Да-с, теперь к твоей заметке о возможном фильтре. Я уже отметил — спектр помехи никто толком не знает и получить не умеет, такая у нас ситуация. Поэтому и считаю, что плясать от заданной АЧХ — непонятно от чего же плясать. Так что старика Чебышева сюда вроде как и не притянуть. Кстати, первый порядок всегда одинаков, он не может быть ни чебышевским, ни бесселевским, ни каким иным.
А более высокие порядки — хороши, но как их считать, если есть минимальная информация о сигнале и шуме… Вот я и применяю скользячку: что-то уж точно усреднит-сгладит, места не требует, ресурсов минимум, да и методика расчета фильтра — вот она, простая до слез.
Вот, скажем, попросил бы я по твоей методике посчитать усреднитель, про который известно только одно: переходной процесс с точностью до (0,1...5)% пусть закончится через заданное время. И чтоб попроще. Какой будет ход расчетов?
0
Вы, безусловно правы, помеха на то и помеха, что про неё мало известно (ну кроме шума от розетки, конечно). Я бы решал задачу в несколько шагов.

1) Если примерно известно время изменения целевого сигнала, то меня бы интересовали частоты меньшие 2*ПИ/<это время>.

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

Безусловно, минус в том, что нельзя рассчитать коэффициенты прямо в микроконтроллере. Точнее можно, но нецелесообразно.

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

Понятно, контретно применение бегушего среднего или более высокого порядка зависит от конкретной ситуации — вида целевого сигнала и помехи.

Вот, вроде получилось свою мысль словами изложить :) Извините на многобуквенность.
0
Книжка хороша! Спасибо, скачал.
А ссылкой поделиться?)
0
Дык, гугель все знает. Порылся в своей истории, вот это, кажется:
iitkgp.vlab.co.in/userfiles/9/file/PADFD.pdf
0
делаю проще
value -= value>>6; /// чем больше сдвиг
value += new_adc_value<<6; /// тем сильнее сглаживает.
0
Так… А чем это проще? Это так же по вычислениям, только у тебя БЕТА из ряда отрицательных степеней двойки.
И как-то не наглядно с порядком операций. Хорошо, если ты уже проверил и мин нет :)
0
проверено
0
Спасибо.
Кстати, на другом форуме один товарищ тоже предлагал ограничить ряд значений БЕТА только степенями двойки. Ну ясно, что у меня предлагается более богатый вариант. Но! Я задумался: а какова будет ситуация с достижимыми параметрами фильтра, если вот такую епитимью на себя возложить? То есть, какое множество фильтров получится, если перебрать все БЕТА от 1/2 до 1/256? Это же легко посчитать, всего 8 фильтров. А может для простых применений этого достаточно? Может удобна будет такая табличка:
Номер фильтра (степень двойки) Тактов до окончания перех.процесса
1 3
2 10
3 50

(числа от балды!)
Так что спасибо за то, что укрепил меня в этой мысли. Доделаю. Не сегодня :)
0
Коэффициенты фильтра вычислены. Поставив их в утилиту, описанную в заметке коллеги mc2, можете полюбоваться частотной характеристикой.
Подскажите, а как подставить коэффициенты в эту утилиту? Я там вижу только поля для ввода частоты сэмплирования и частоты среза…
0
М-дя… Ошибочка вышла. Не та утилита. Там только подсчитать к-ты можно. Ну, значит другой программой для расчета АЧХ, Матлаб и прочие.
0
спасибо! очень помогло. в рантайме персчитываю константы, т.к. должно меняться время «усреднения», на stm32 тормозов не замечено. всё считается во флоате.
0
а про более высокие порядки фильтрации где мона почитать?
0
Запрос «рабинер л. гоулд б. теория и применение цифровой обработки сигналов» дает ссылки на ту книгу, по которой я лет сто назад учился. Были и другие, не вспомню теперь.
Толстая такая… :)
0
Литературы много, очень много.
Можно начать с: Богнер, Консиантинидис, «Введение в цифровую фильтрацию»
digteh.ru/LIB/books.php
или Жан Макс, «Методы и техника обработки сигналов».
Для самых «маленьких»: Юкио Сато, «Обработка сигналов: первое знакомство».
easyelectronics.ru/yukio-sato-obrabotka-signalov-pervoe-znakomstvo.html
0
Жан Макс: depositfiles, letitbit.
0
Автору, нельзя ли привести статью в соответствии со временем?
Я о картинке с подписью "(проба пера: изображение формул через www.texify.com)".
Изображения нет, при попытке зайти на www.texify.com выдает это:

Спасибо.
0
В соответствии со временем, говоришь? Да, времена изменились.

Там весь фокус был в дублирующем представлении формулы. Убрал. Осталось все, что необходимо.
0
Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.