AVR. Учебный Курс. Кусочно-линейная аппроксимация

AVR. Учебный курс 48 комментариев

Часто бывает так, что приходится обрабатывать жутко нелинейные величины, задаваемые каким-нибудь извращенным законом. Простейший пример — датчики расстояния SHARP GP2D12. Только поглядите на его характеристику:

Сам черт ногу сломит, а ведь нам бы неплохо иметь выход в человеческих величинах, ну или, хотя бы, линейно зависящие от расстояния. Что делать?

Вариантов тут, на самом деле, всего два. Первый очень быстрый, но жадный до памяти ПЗУ — табличный.
То есть мы просто берем и эту кривулину расписываем в памяти. Например, у нас с 8ми разрядного АЦП идет значение напряжения от 0 до 256, а мы на каждое значение создаем в памяти значение расстояния. Тогда с АЦП сразу гоним в индекс массива, где эти значения хранятся и получаем расстояние:

L=Curve[ADCH];

Недостаток один — прожорливость до памяти, растущая в геометрической прогрессии с ростом разрядности АЦП.

Вариант второй — написать функцию, переводящую одну величину в другую.
И вот тут можно поизвращаться. Дело в том, что если попытаться первести данные напрямую, как нам это советует даташит, на этот дальномер, то мы получим выражение вида:

L = k*X +b — уравнение прямой.

По даташиту видно, что характеристика становится линейной при выворачивании её наизнанку:

То есть у нас получается:

V = 1/L*k+b

  • V — мы знаем, это показание с АЦП, помноженные на вес разряда: ADC * digit_weight. АЦП в данный момент имеет опорное напряжение в 3.3 вольта и 8 бит. Так что вес будет 3.3/255 = 0,013 вольта
  • k — судя по графику, равно около 13.
  • b — на вскидку, если продолжить график до нуля, выглядит примерно на 0.1

Получили:

0.013*ADC = 1/L*13 + 0.1

Путем применения технологий жуткого матана :) мы получаем чистую формулу расстояния:

L = 13/(0.013ADC-0.1)

Прикинув по исходному графику, убеждаемся что нигде ничего не забыли и не напутали.
Теперь покоцаем все на 0.013 и округлим результаты. Пока мы считаем сами нам доступна эта роскошь. МК же в целочисленных вычислениях округлять не умеет.

L = 1000/(ADC-8)

Вот, совсем красиво получилось. Одна непрятность — деление то у нас целочисленное, а значит будут потери округления. Т.е. при целочисленном вычислении 9.9 и 9.1 округляются до 9. Непорядок.

Наиболее точный результат целочисленного деления выражения А=B/C дает формула A = (B+C/2)/C

L = (1000+(ADC-8)/2)/(ADC-8)

В принципе, так можно и оставить и в таком виде скормить компилятору. Но можно и пооптимизировать. Во первых, что в первую очередь просится под нож, так это операция деления — самая мерзкая и тормозная из арифметических операций.

Как минимум одно деление можно оптимизировать:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

u08 conversion (u08 _ADC)
{
u08 Half_DIV;
u16 Result;
 
if(_ADC<15) return 255;		// Очень далеко. Или ошибка датчика.
 
Half_DIV = _ADC-8;
Half_DIV >>=1;			// Делим пополам сдвигом. 
 
Result = 1000+(u16)Half_DIV;	
 
return (u08)(Result/(_ADC-8));	
}

Зачем в самом начале проверку на минимальное значение? Да на всякий случай, чтобы не было деления на ноль, например. Или каких то космических значений. Видишь же, что характеристика датчика по напряжению не опускается ниже 0.2 вольт. А 0.2 вольта в наших показаниях АЦП это 15. Можно и вообще ограничить зону точкой четкого определения расстояния, скажем, 30 сантиметров. Это 0.4 вольта, тогда ограничивать можно сразу от 30 АЦПшных попугаев.

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

Что делать? А тут на помощь опять придет злой матан в лице аппроксимации кусочно-линейными функциями. За злой терминологией скрывается, на самом деле, примитивнейшая вещь:

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

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

Я отрисовал его по точкам и забросил в маткад, взяв только активные области.
Вот только не понравилось мне то, что входное значение меняется от 0 до 255, а выходное от 0 до 30. Да, на выходе сантиметры, но я бы предпочел полный размах от 0 до 255. Зачем? А шумы фильтровать проще будет.

А еще повернул так, чтобы по оси Х были данные из АЦП, а по оси У расстояние в наших попугаях (1 попугай = 1/7см).

Итак, возьмем и отмасштабируем нашу кривулину до 256, просто помножив на 7 с копейками.

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

Я взял три куска, описав их следующими функциями:

  • X = -4Y + 320 при Y от 0 до 56
  • X = -Y + 150 при Y от 56 до 118
  • X = -0.25Y +60 при Y от 118 до 250

Коэффициенты я старался подобрать такие чтобы прямая не только максимально плотно легла на исходную кривую, но и чтобы они были как можно ближе к степени двойки. Чтобы можно было делить и умножать сдвигами. Если не удается подобрать степень двойки, то можно близкие к ней. Например, Y*3 = (Y<<1)+Y Ну и в таком духе.

Вот, значит, моя спрямленная кривая.

А теперь, собственно, код:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

u08 LinearAPPROX(u08 input)
{
// Too far
if (input<21) return 255;
 
// Line1 X=63+255-4Y
if(input<56)
	{
	input <<=2;         		//4*Y
	return (63+(255-input));
	}
 
// Line2 X=150-Y
if (input<118)
	{
	return (150-input);
	}
 
// Line3 X=63-0.25Y
if(input<250)
	{
	input >>=2;         		//0.25*Y
	return (63-input);
	}
 
// Too close
return 0;
}

Заняло около 90 байт флеша, а выполняется тактов за 15-20 по самой длинной ветке. И еще один такой немаловажный момент. Обратите внимание на функцию X = -4Y + 320 и на то как я ее представил в программе: X=65+255-4Y. Зачем я разбил число 320 на два куска? Да дело в том, что моя выходная величина не вылазит за предел беззнакового восьмиразрядного — это 255. А вот число 320, требуемое по уравнению, в этот лимит не влазит. И тут либо компилятор окажется умным и сам догадается, либо что нибудь куда-нибудь переполнится, что вероятней. Хорошо если Warning даст. А если не даст? Поэтому я решил подстелить соломки, сразу указав в какой последовательности и с числами какого размера оперировать, чтобы не было переполнений. Да и лишней движухи на уровне кода меньше будет.

На чистом ассемблере вообще загляденье.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

; Input: 
; OSRG = R17 - input data
; ACC = R16
; Output
; OSRG - output data
LineAPPROX:	MOV	ACC,OSRG	; OSRG = ACC для удобства операций
 
 		CPI	OSRG,20		;if(OSRG<20) слишком далеко
		BRCS	TooFar
 
		CPI	OSRG,56		;if(OSRG<56) далеко
		BRCS	Far
 
		CPI	OSRG,118
		BRCS	Near		;if(OSRG<118) средняя дистанция
 
		CPI	OSRG,250	;if(OSRG<250) близко
		BRCS	Close
 
TooClose:	LDI	OSRG,0		; совсем близко, сразу 0
		RET			; Выходим с результатом в OSRG
 
; X = 63-0.25Y	
Close:		LSR	ACC		;0.25Y сдвигом делим на 4
		LSR	ACC
 
		LDI	OSRG,63		;63 - (0.25Y) вычитаем
		SUB	OSRG,ACC
		RET			; Выходим с результатом в OSRG
 
; X= 150 - Y
Near:		LDI	OSRG,150
		SUB	OSRG,ACC	; 150-Y Не зря мы в АСС скопировали
		RET			; Выходим с результатом в OSRG
 
; X = 63+255-4Y
Far:		LSL	OSRG		; Y*4 сдвигом умножаем на 4
		LSL 	OSRG
 
		COM	OSRG		; 255-(Y*4)  вычли
		SUBI	OSRG,-63	; 63+(255-(Y*4)) сложили
		RET			; Выходим с результатом в OSRG
 
TooFar:	LDI	OSRG,255		; Совсем далеко? Ну и думать нечего!
		RET			; Выходим с результатом в OSRG

Итого, 52 байта на весь экшн. Причем на самую длинную ветвь отводится всего 10 тактов. Красота! Можно даже в прерывании обрабатывать! :)

Да, к вопросу о том, зачем мне понадобилось растягивать диапазон до 255 попугаев. Дело в том, что несколько повышается точность выхода, но правда увеличивается уровень болтанки сигнала. Т.к. ИК порой выдает много мусора, особенно на дальних расстояниях, когда у него вытягивается характеристика и малейшие колебания превращаются в мощный расколбас выходного значения.

Так вот, когда у нас диапазон большой, то можно нажрать этой каши побольше, потом усреднить и сразу же поделить на попугаев (в том же сдвиге). На выходе, по идее, получим очень стабильные и красивые сантиметры. Впрочем, это только первичный экспериментальный прогон. Я толком еще не тестировал, но первые опыты уже дали более менее адекватные результаты.

Спасибо!!! Вы потрясающие! Всего за месяц мы собрали нужную сумму в 500000 на хоккейную коробку для детского дома Аистенок. Из которых 125000+ было от вас, читателей EasyElectronics!!! Были даже переводы на 25000+ и просто поток платежей на 251 рубль. Это невероятно круто!!! Сейчас идет заключение договора и подготовка к строительству!

А я встрял на три года, как минимум, ежемесячной пахоты над статьями :)))))))))))) Спасибо вам за такой мощный пинок!!!

48 thoughts on “AVR. Учебный Курс. Кусочно-линейная аппроксимация”

  1. Спасибо, DI HALT.
    мне такой информации не хватало, т.к. матан я чесно прогулял, когда учился в универе

    Ответить

  4. Был я техником одной заумной системы. Там Была такая штука, как BWO (VCO в микроволновом диапозоне), но не суть. Управлялся BWO линейным напряжением, но выходные характеристики были далеки от линейности.
    Так умные инженеры поставили между управляющим напряжением и BWO вещь, которую назвали Lineaser. С помощью небольшой матрицы из резисторов управляющее напряжение превращалось в обратную функцию характеристики BWO и в итоге если сложить 2 графика, получалась линейная функция.
    Думаю, в данном случае с датчиком, тоже можно использовать Lineaser, если поставить его между датчиком и контроллером. Заморочиться аппаратно. Тогда не надо будет заморачивать программу. Хотя, каждый выбирает для себя, что лучше. Но это тоже один из способов.

    Ответить

  5. Если уж сказали аппроксимация, надо было вспомнить и про оценку.
    По хорошему коэф. надо подбирать не по графику, а хотя бы использовать метод наименьших квадратов.
    А по поводу коэф. являющихся степенями двойки — делал такую аппроксимацию, начал примерно так-же — нарисовал график, нарисовал возможные варианты прямых, у которых коэф. равны степеням двойки, потом написал программу, которая выбрала отрезки и соотв. им уравнения прямых. Правда результат после округления не сильно отличался от выбранных вручную =/

    Ответить

  6. «Получили: 0.013*ADC = 1/L*13 + 0.1
    Путем применения технологий жуткого матана :) мы получаем чистую формулу расстояния:
    L = 13/(0.13ADC+0.1)»

    было 0.013 стало 0.13. шайтан-матан, однако….

    Ответить

  8. Вообще-то неплохо было бы посчитать погрешность аппроксимации, и показать ее отдельным графиком. И объяснить почему аппроксимация лучше, чем прямое вычисление функции (ну там меньше места занимает, выполняется быстрее итд.)
    А матана тут нет нифига. Алгебра, 9-й класс.

    Ответить

    1. Как это нету нифига? А умные термины? Я их впервые услышал только на Численных Методах в курсе матана в универе. ;)

      Хм. Насчет графика мысль хорошая. Ну а о причинах перехода к аппроксимации я вроде уже упомянул — упрощение и облегчение вычислений.

      Ответить

  10. <>
    Должно быть L = 1000/(ADC-8).

    Ещё можно апроксимировать полиномами второго, третьего итд порядков. Там, в принципе, ничего страшного нет — всего лишь умножения и сложения. Полином по экспериментальным данным нам посчитает тот же Octave или Mathcad (функция polyfit()).

    Ещё один вариант — интерполяция по сплайнам, в частности по сплайну Кэтмула-Рома. Достаточно иметь с десяток экспериментальных точек (если апроксимируемая функция не совсем бешенная).

    Ответить

    1. А он еще и не произошел. Точнее мы стоим сейчас там где и будем стоять, но соседей от нас на другой сервер еще не переселили. А потому работаем с жесткими лимитами на динамику (не более 5одновременных коннектов к апачу). Через несколько дней законфигурят сервак и освободят нашу площадку (кроме нас щас тут еще около 50 сайтов %) )

      Ответить

      2. А не проще что-нибудь типа http://fastvps.ru/ взять. Не в рекламных целях.
        Но даже самая простая тачка лихвой справится с задачей.
        Понятно, что денег стоит, но предлагаю сделать пожертвования (ибо думает мя что наберем с миру по нитки :) сначала и накопить нужную сумму. Даже потестить, вроде, можно. И будет нам площадка и места и всего-всего.
        P.S.: Как раз хочу ближайшее время у них заказать аренду. Можем состыковаться.
        P.P.S.: Не имею ничего против текущего хостера, но всё же как-то проект растет и вроде стоит задумываться об егойным будущем

        Ответить

        1. Ну мы с хостером договорились примерно на то же, но дешевле. Плюс хостера и всю его команду я знаю лично. Просто выделить сервер потребовалось ВНЕЗАПНО, а на это надо время.

          Ответить

          1. Ну тогда да. Если дают отдельный сервак — это есть гуд. Даже всегда гуд. Жду новоиспеченного сервака.

            Ответить

  14. Кто может поделится эфиктивными методами вычисления различных математический функций на ассемблере.

    Ответить

    1. DIHALT может поделиться, когда откроет например отдельную рубрику типа: «Математика в AVR». Потому как это дело в одной — двух статейках не опишеш. Да и усилий потребует не мало и времени. Второй вариант самому собирать с миру по нитке и вникать мозгом. Часть какая то уже есть в курсе по AVR. Тут еще можеш зазырить http://www.mymcu.ru/Articles/Articles.htm К тому же вижу что сейчас DI приналег на ARM микроконтроллеры, поэтому математику пока лучше самому ковырять.

      Ответить

  16. DI HALT День Добрый.
    У меня возникла проблема похожая на эту тему.Не подскажете как можно организовать:
    10-бит АЦП при 2.5 вольта примерно 500, мне так и надо,но начинать мерить я хочу
    от примерно 1-го вольта это 000.В общем начало измерений сдвинуть на 1 вольт.
    Программно это возможно или просто ADC GND приподнять на 1 вольт и с помощью
    AVREF подогнать потом значение при 2.5 вольта.В общем изменить линейность
    измерений как у потенциометра существует три линейности A B C.Проблема в том
    у меня датчик реостат при максимуме я имею к.з. пр минимуме примерно 400 Ом
    перевернуть невозможно поэтому появился приподнятый 0. Пишу в ASM

    Ответить

    1. Хочу извениться за свою невнимательность.Плохо разобрался о чом идёт речь.
      С нолевым положением разобрался.С преобразованием папугаев проблема.Пытаюсь
      собрать с само корректирующейся линейностью .Должно указываться что это 100
      дальше 200 итд.
      r26,r27 АЦП
      r12,r13 нолевое значение (откуда начинаем измерения)
      r14,r15 максимальное значение
      и всё это делал так:
      r14,r15 минус желаемое показание 600-500=100
      r14,r15 делю на результат 600/100=6
      а дальше просто
      АЦП делю на 6 300/6=50
      АЦП отнимаю полученное значение 300-50=250

      max:
      sub r26,r12 отнимаю от АЦП нолевое значение
      sbc r27,r13

      sub r14,r12 отнимаю от максимального нолевое значение
      sbc r15,r13

      mov r24,r26 делаю копию
      mov r25,r27

      mov r18,r14
      mov r19,r15

      mov r4,r14
      mov r5,r15

      subi r18,244 отнимаю от максимального значения желаемое показание
      sbci r19,1

      mov r10,r14
      mov r11,r15

      max1: cp r10,r18
      cpc r11,r19
      brlo max2
      sub r10,r18
      sbc r11,r19
      inc r8
      brne max1
      inc r9
      rjmp max1

      max2: cp r24,r8
      cpc r25,r9
      brlo max3
      sub r24,r8
      sbc r25,r9
      inc r6
      brne max2
      inc r7
      rjmp max2
      max3: sub r26,r6
      sbc r27,r7

      ret

      Это все хорошо для целых чисел .Не подскажете как можно округлять не-то
      погрешность через-чур большая.Всё зависит от того как карта ляжет.

      Ответить

      1. Доброго время суток.Проблему я решил путём разгона АЦП до 16 бит ,а потом
        делить до желаемого результата .Таким образом можно калибровать под любой датчик
        и необязательно знать его таблицу линейности ,просто создаётся желаемое количество ступеней в таблице ,чем больше ступеней тем больше точность.
        Вся калибровка помещается в Эпром ,при включении записывается в СРАМ.
        Премущество в отсутствии каких либо подстроечных резисторов.Скорсть правда не космическая!Для показаний в ЛЦД запаса хватает надо ещё притармаживать.
        При желании могу выложить.

        Ответить

  17. Кстати, не всё так просто тут =).Как определить минимально необходимое число участков КЛА, если допустимая приведённая погрешность аппроксимации равна x?
    http://pr-1.p0.ru/forum/7-43-2
    Автоматизация инженерных расчетов:
    Щепетов А.Г. Автоматизация инженерных расчетов в среде MATHCAD: Учебно — практическое пособие.

    Ответить

  18. Di Halt, помоги пожалста с делением.
    делаю огибающую для пилы типа ADSR
    Атмел предлагает такой алгоритм деления 8ми разрядных чисел

    1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

    ;***** Subroutine Register Variables
 
.def	drem8u	=r15		;remainder
.def	dres8u	=r16		;result
.def	dd8u	=r16		;dividend
.def	dv8u	=r17		;divisor
.def	dcnt8u	=r18		;loop counter
 
;***** Code
 
div8u:	sub	drem8u,drem8u	;clear remainder and carry
	ldi	dcnt8u,9	;init loop counter
d8u_1:	rol	dd8u		;shift left dividend
	dec	dcnt8u		;decrement counter
	brne	d8u_2		;if done
	ret			;    return
d8u_2:	rol	drem8u		;shift dividend into remainder
	sub	drem8u,dv8u	;remainder = remainder - divisor
	brcc	d8u_3		;if result negative
	add	drem8u,dv8u	;    restore remainder
	clc			;    clear carry to be shifted into result
	rjmp	d8u_1		;else
d8u_3:	sec			;    set carry to be shifted into result
	rjmp	d8u_1

    что то тут не то, ибо он должен помещать значение в регистр dres8u, но в коде даже не вспоминает про него. что не так?

    Ответить

  19. «Наиболее точный результат целочисленного деления выражения А=B/C дает формула A = (B+C/2)/C»
    Это как так может получиться??? Чего-то тут с теорией не того: (B+C/2)/C = B/C + 0.5 != B/C

    Ответить

  20. L = 1000/(ADC-8)
    “Наиболее точный результат целочисленного деления выражения А=B/C дает формула A = (B+C/2)/C”

    А может в этом конкретном примере, при B=1000 такое будет еще точнее?
    L = (64000/(ADC-8))>>6

    Ответить

  21. Здравсвуй DI HALT!

    Какие меры можно предпринять для защиты от покосячинга EEPROMа при снижении напряжения питания? Я так понимаю в первую очередь это монитор питания. А программно? И дадут ли эти меры 100% гарантию?

    И вот ещё что… наткнулся я на просторах интернета на сайт allelectronics.3dn.ru
    Дык вот его дизайн слизан с вашего + в разделе КНИГИ коменты к книгам тоже ващи… Короче воруют…..

    Ответить

    1. 100% гарантии не даст никто. Если нужна полная гарантия, то лучше использовать внешнюю ЕЕПРОМ. Т.к. при ее цикле записи запортачить содержимое случайным сбоем очень сложно.

      Да кто только у меня не тащит статьи. Что делать — популярность. А с сайта этого мы всем сообществом оборжались.

      Ответить

Добавить комментарий

Ваш e-mail не будет опубликован.

Этот сайт использует Akismet для борьбы со спамом. Узнайте как обрабатываются ваши данные комментариев.