Статистическая обработка медленно изменяющихся сигналов
Два раза ку, дорогие товарищи!
Немало интересного сказано по математической обработке сигналов. Рискну добавить.
Обрабатывая сигнал с АЦП простым процессором нам в общем случае нужно решить три задачи:
- сглаживание сигнала
- увеличение разрядности АЦП
- достижение 1 и 2 пункта с минимальными затратами ресурсов микроконтроллера.
Самым простым, не требующим памяти фильтром является фильтр модифицированного скользящего среднего. Математически он задается так:
(1) Ki+1 = (N * Ki + Xi) / (N + 1)
Ki - текущее значение фильтра
Ki+1 - новое значение фильтра
Xi - текущее значение параметра
N - коэффицент инерции фильтра, он же постоянная числа итераций этого фильтра. Проще говоря, если сигнал Х в какой-то момент изменил уровень, то каждое N раз измерений значение фильтра будет в e (2,718..) раз ближе к сигналу Х (грубо - на 63%). Постоянная числа итераций связана с постоянной времени фильтра τ (тау) через частоту f опроса датчика: τ=N/f
Проще говоря, если, скажем, температура объекта изменилась скачком с 0 до 100°, а N = 8 и измерения мы проводим 1 раз в секунду, то через 8 секунд фильтр даст 63°, еще через 8 будет 88, потом 95, далее 98, 99... акцентирую внимание на том, что в данном случае N - это не постоянная времени, а постоянная числа измерений. Ведь если частота опроса увеличится в 10 раз, то и значения фильтра будут изменяться в 10 раз быстрее.
Однако даже в таком виде эта формула сложновата для процессора, ибо приходится применять долгие операции деления, а в худшем случае - и числа с плавающей точкой. И вот уже код не помещается в память, а время обработки просто зашкаливает.
Операции деления и умножения можно заменить на операции сдвига, если речь идет о степенях двойки. Однако в вышеприведенной формуле, даже если установить N=2ˆR, то все равно придется делить на N + 1. По этому слегка упростим себе задачу, воспользовавшись формулой
(2) K / (N + 1) ≈ (K - (K / N)) / N
как видно, в данном случае мы имеем деление на степени двойки, которое можно заменить сдвигом. Таким образом, мы получаем слегка модифицированную формулу, применимую для расчётов на МК. Приведём пару выкладок:
Ki+1 = (N * Ki + Xi) / (N + 1) = N * Ki / (N + 1)+ Xi / (N + 1) ≈ (N * Ki - (N * Ki / N)) / N + (Xi - (Xi / N)) / N = Ki - Ki / N + Xi / N - (Xi / N2)
При этом считаем, что величина Xi / N2 достаточно маленькая, и поэтому ею можно пренебречь - раз уж мы используем округление по формуле (2).
(Для особых эстетов если N=2R, тогда N2=(2R)2=2(R*2). то есть Xi / N2 = Xi shr (R shl 2))
Введём дополнительную переменную:
Pi+1 = Ki+1 * N
Тогда
Ki+1 * N = Ki * N - (Ki * N) / N + Xi
(3) Pi+1 = Pi - (Pi / N) + Xi
(4) Ki+1 = Pi+1 / N
Дополнительная переменная P введена для того, чтобы при целочисленном делении не терять точность. В общем, если N = 2R, то в виде программы это выглядит вот так:
Си:
Filter = Filter - (Filter >> R) + X; Result = Filter >> R;
Pascal:
Filter := Filter - (Filter shr R) + X; Result := Filter shr R;
Если требования к скорости обработки данных высоки, а вот точность не столь важна, то можно пожертвовать одним измерением, но сократить число сдвигов вдвое.
Си:
Result = Filter >> R; Filter = Filter - Result + X;
Pascal:
Result := Filter shr R; Filter := Filter - Result + X;
Теперь, разобравшись со сглаживанием сигнала и выкинув шумы, зададимся вопросом, а может, в этих шумах было что-то полезное, и мы зря их выкидывали? Возможно. Дело в том, что из шума можно легко вытащить дополнительные разряды АЦП.
Если, к примеру, стоит задача измерять температуру термопарой с достаточно большой точностью, но с малой скоростью, то нет смысла ставить специализированную дорогую микросхему типа AD7799 (мою любимую, кстати) можно вполне обойтись и внутренним 10-битным АЦП.
Используя тот факт, что последние разряды АЦП обычно "шумят", мы можем извлечь из этого шума дополнительную информацию о сигнале: если вычислить уровень этого шума, то можно сделать вывод об уровне сигнала. То есть, взяв среднее арифметическое от значений АЦП за достаточно долгий промежуток времени (при неизменном уровне сигнала) мы можем получить его значение с точностью, превосходящей разрядность АЦП.
Иными словами, если значения АЦП "болтаются" между 3 и 4, причем 3 выпадает в 65% случаев, а 4 в 35% соответственно, то можно утверждать что сигнал приблизительно равен 3,65
И вот тут очень пригодится арифметика с фиксированной точкой. То есть когда веса в ячейках имеют нецелочисленные значения. Обычно в байте веса ячеек - это положительные степени двойки.
S = B7*(2ˆ7)+B6*(2ˆ6)+B5*(2ˆ5)+B4*(2ˆ4)+B3*(2ˆ3)+B2*(2ˆ2)+B1*(2ˆ1)+B0*(2ˆ0)
S = B7*(128)+B6*(64)+B5*(32)+B4*(16)+B3*(8)+B2*(4)+B1*(2)+B0*(1)
Однако если принять что нулевая ячейка имеет показатель степени не 0, а, к примеру -2 то получим следующую картину.
S = B7*(2ˆ5)+B6*(2ˆ4)+B5*(2ˆ3)+B4*(2ˆ2)+B3*(2ˆ1)+B2*(2ˆ0)+B1*(2ˆ-1)+B0*(2ˆ-2)
S = B7*(32)+B6*(16)+B5*(8)+B4*(4)+B3*(2)+B2*(1)+B1*(0,5)+B0*(0,25)
в результате число S представлено уже с точностью до 0,25. Что самое замечательное в данном представлении - это то, что процессор не замечает такой подставы. Все действия с числами с фиксированной запятой полностью идентичны действиям с целыми числами.
В нашем случае представим что нам нужно повысить разрядность АЦП на Q разрядов. Это значит что нам надо повысить точность в M раз. M=2ˆQ. Тогда наша формула усреднения приобретает вид
(5) Pi+1 = Pi - (Pi / (N * M)) + Xi * M
(6) Ki+1 = (Pi+1 / N) / M
Следует однако учесть, что в формуле (6) деление /N можно заменить сдвигом, а вот /М уже придется делать с помощью математики с плавающей точкой, иначе мы просто потеряем точность. Однако утешает тот факт, что формулу (6) нужно применять не каждый раз, когда считываются данные ацп, а лишь при выводе результата измерения.
Если же речь идет не об индикации а об управлении каким-то объектом (к примеру, ПИД-регулятор), или выводе на бинарный индикатор (светодиодную линейку), то приняв фиксированную точку с точностью M для всей программы можно вообще не делить на M.
Теперь немного о разрядности переменных, участвующих в фильтрации. Если разрядность АЦП принять как L, то разрядность переменной P должна быть не менее L+Q+R
Удачи, дерзайте.