nabbla (nabbla1) wrote,
nabbla
nabbla1

Categories:

Тестируем atan1 на QuatCore

Пора уже перебираться на "железо" потихоньку. Решил начать с самого первого алгоритма, поскольку он уже был написан на ассемблере.

В программу внёс тестовый кусочек, где берётся 30 единичных векторов, распиханных по окружности с шагом 12° (от 0° и до 348°), и от каждого из них последовательно берётся atan1. Вектора просто засунул в неиспользуемую область памяти
[Spoiler (click to open)]
ORG	0x1A0

;тестовые данные для atan1
;набор из 30 векторов, через каждые 12 градусов, проверить работу на всей окружности
TestVec0x	Int16	32767
TestVec0y	Int16	0
TestVec1x	Int16	32052
TestVec1y	Int16	6813
TestVec2x	Int16	29935
TestVec2y	Int16 13328
TestVec3x	Int16	26510
TestVec3y	Int16	19261
TestVec4x	Int16	21926
TestVec4y	Int16	24351
TestVec5x	Int16	16384
TestVec5y	Int16	28378
TestVec6x	Int16	10126
TestVec6y	Int16	31164
TestVec7x	Int16	3425
TestVec7y	Int16	32588
TestVec8x	Int16	-3425
TestVec8y	Int16	32588
TestVec9x	Int16	-10126
TestVec9y	Int16	31164
TestVec10x	Int16	-16384
TestVec10y	Int16	28378
TestVec11x	Int16	-21926
TestVec11y	Int16	24351
TestVec12x	Int16	-26510
TestVec12y	Int16 19261
TestVec13x	Int16 -29935
TestVec13y	Int16 13328
TestVec14x	Int16	-32052
TestVec14y	Int16 6813
TestVec15x	Int16 -32768
TestVec15y	Int16 0
TestVec16x	Int16 -32052
TestVec16y	Int16 -6813
TestVec17x	Int16 -29935
TestVec17y	Int16 -13328
TestVec18x	Int16 -26510
TestVec18y	Int16 -19261
TestVec19x	Int16 -21926
TestVec19y	Int16 -24351
TestVec20x	Int16 -16384
TestVec20y	Int16 -28378
TestVec21x	Int16 -10126
TestVec21y	Int16 -31164
TestVec22x	Int16 -3425
TestVec22y	Int16 -32588
TestVec23x	Int16 3425
TestVec23y	Int16 -32588
TestVec24x	Int16 10126
TestVec24y	Int16 -31164
TestVec25x	Int16 16384
TestVec25y	Int16 -28378
TestVec26x	Int16 21926
TestVec26y	Int16 -24351
TestVec27x	Int16 26510
TestVec27y	Int16 -19261
TestVec28x	Int16 29935
TestVec28y	Int16 -13328
TestVec29x	Int16 32052
TestVec29y	Int16 -6813



и сразу после инициализации стека добавил вот этот код:
;здесь хотим проверить работу atan1
		Y		Matrix	;пока не используем, сойдёт
		k		29		;чтобы замучать 30 векторов
		X		TestVec29x
		Z		RT_PRF_NO	;отсюда начнут углы лежать (заголовок пожалеем)
@@k_loop:	CALL		atan1
		[Z+k]		[Y+1]	;перенесли угол куда хотели
		Acc		X
		SUB		2
		X		Acc
		kLOOP		@@k_loop


Я не питал иллюзий, что оно с первого же раза заработает правильно, поэтому не стал пока прошивать в ПЛИС, а запустил на симуляции. Очень удачно на осциллограмме мне как раз показали участок, где мы возвратились из процедуры atan1. И там я увидел, что регистр X указывает куда-то не туда.

Семён Семёныч: я забыл в atan1 всё-таки раскомментировать код, отвечающий за сохранение регистров в стек и восстановление их значений в конце процедуры, вот регистр X и "затирался".


Львиную долю этого кода я раскомментировал, и на данный момент atan1 выглядит так:

;atan1(y,x)
;берёт atan(y/x), результат от -pi/2 до pi/2, в формате Q2.12 (т.е представимые углы от -2 до +2)
;также находит длину вектора (x,y). 

;в [X], [X+1] лежат значения x,y, их трогать не будем
;в [Y] после выполнения будет лежать длина вектора, делёная на 2, а в [Y+1]: угол. 

;точность нахождения угла: все 16 бит, если (x,y) нормированы (т.е x^2+y^2=1), в противном случае не хуже 5,625 градусов (5..6 бит значащих) 
;точность нахождения длины вектора: 0,36% (8..9 значащих бит)

atan1 proc

		[SP++]	X	;будет изменён в процессе
		[SP++]	Z	;будет указатель на таблицу углов
		[SP++]	ijk	;все регистры в хвост и в гриву
		; [SP++]	C
					; ;а вот Z оставим нетронутым! Ну и Acc считаем "временным", равно как и флаги
			
;assign ijk = {invreg, kreg, ireg, jreg};
		ijk		0x0804		;Inv=0, k=2, i=0, j=4
		ZAcc		RoundZero
		SUB		[X+i]		;то есть просто [X]
		[SP]		0		;угол инициализируем нулём
		Z		AtanTable
@@j_loop:	Inv		S
		C		[X+i]		;то есть просто [X], здесь всегда i=0
		i		1	
@@i_loop:	ZAcc		RoundZero
		FMPM		[Z+3j+k]	;то есть [Z+3j+2], он же "синус" (мнимая часть). т.к k=2 на всём протяжении
		C		[X+i]
		FMA		[Z+3j+1]	;он же "косинус" (действительная часть)
		[Y+i]		Acc
		iLOOP		@@i_loop	
;вышло 8 строчек, но пока фигня с j (когда не равно единице, со знаками песец выходит)
		X		Y	;возможно, X сохраним в стек, чтобы процедура минимально меняла регистры
		Acc		[SP]
		PM		[Z+3j+i]	;то есть [Z+3j], угол. 
		[SP]		Acc		
		jLOOP		@@j_loop
		DIV2A		[Y+i]
		[Y+i]		Acc
		
		; C		[--SP]
		ijk		[--SP]
		Z		[--SP]
		X		[--SP]

		JMP		[--SP]
atan1 endp


Разве что регистр C пока что решил не сохранять. Скомпилировал повторно, отсинтезировал модуль PipelineQuatCore (это "голое ядро", без периферии), снова запустил симуляцию. На этом же самом месте увидел, что теперь регистр X ведёт себя адекватно, но вот полученный результат, 0x7FFF - ОЧЕНЬ ПОДОЗРИТЕЛЬНЫЙ! Дело в том, что у нас ответ должен выдаваться в формате Q2.14, т.е в диапазоне от -2 до +2 радиан, но ответ должен лежать в более узком диапазоне, от -π/2 до +π/2, т.е уж точно не заходить на упор!

Тут, похоже меня подвели устаревшие комментарии. По ним выходило, что в [Y] будет лежать длина вектора, а в [Y+1] - угол. Но глядя на код и вспоминая, как он работает, я думаю, что ВСЁ НАОБОРОТ: в [Y] лежит угол, а в [Y+1] длина, и она действительно должна получаться единичной всё это время.

Кстати, за 500 мкс симуляции он так и не успевает посчитать все значения (да, этот алгоритм МЕДЛЕННЫЙ), и дамп памяти по окончании 500 мкс выглядит так:



Чуть выше синим обведены те значения, которые уже посчитаны. Тут в основном 0x7FFF, но попадаются чуть меньше. Самое маленькое: 0x7F8A = 32650, что соответствует 0,9964. В комментариях указано, что длина меряется с точностью 0,36% - вот здесь эти самые 0,36% и выходят!

Ладно, заменяем [Y+1] на [Y+i] (у нас должно быть i=0) там, где мы берём результат и помещаем его в "массив", и попробуем ещё разок.

Хотя, самое последнее значение в стеке (начало стека также отмечено синим в дампе, в самом низу) не вселяет оптимизма. Там по идее должен лежать угол поворота, набранный из 45°, 22,5°, 11,25° и так далее, с разными знаками, но ещё без добавления x. Мы уже вовсю идём по I квадранту, вроде углы должны быть положительные, но тут мы видим 0xA1C1, это отрицательное значение -24127, что для Q2.14 соответствует -1,472 рад, или -84,37°. Чушь какая-то...

Ну ладно, тем временем проект отсинтезировался по-новой, симуляция выполнилась, и смотрим ещё один дамп:


Какая-то логика во всём этом проглядывает. Начнём с конца. 0x0D67 = 3431, что в Q2.15 означает 3431/16384 = 0,209 радиан или 11,9984° Правильный ответ: -12°. Где-то перепутали знак (в нечётном количестве раз). А вот в остальном, ошибка составляет 5,75 угловой секунды, это БОЛЕЕ ЧЕМ ВДВОЕ меньше цены младшего разряда (ЦМР), т.е по сути точнее дать ответ мы и не могли!

Следующее значение: 0x1ACE = 6862, или 0,4188 радиан, или 23,9968°. В этот раз ошибка 11,505 угловых секунд, т.е тут всё-таки ошиблись на единичку: лучше было поставить 6863. Ну и знак, разумеется, не тот.

Давайте дальше в таблицу сведём, что получилось и что должно было получиться:


Жёлтым я подсветил, где ошибка превысила 0,5 ЦМР, т.е ответ можно было представить точнее в 16-битном представлении. А красным - где ошибка превысила 1,5 ЦМР, т.е мы аж НА ДВЕ ЕДИНИЦЫ ошиблись. Такое было ровно один раз.

После нуля видим полное повторение значений, "пошли на второй круг". Оно и понятно, atan(y/x) проходит все свои значения за половину окружности, а затем повторяет их. (в отличие от atan2(y,x))

Как же так мы со знаком напортачили...

Давайте самый первый вызов atan1 посмотрим на симуляции:


Для понимания происходящего, листинг:
atan1 proc
127  F3CD              [SP++]  X   ;будет изменён в процессе
128  F3ED              [SP++]  Z   ;будет указатель на таблицу углов
129  F3A4              [SP++]  ijk ;все регистры в хвост и в гриву
12A  A757              ijk     0x0804      ;Inv=0, k=2, i=0, j=4
12B  8801              ZAcc    RoundZero
12C  83C4              SUB     [X+i]       ;то есть просто [X]
12D  FC0E              [SP]    0       ;угол инициализируем нулём
12E  ED58              Z       AtanTable
12F  A381  @@j_loop:   Inv     S
130  8AC4              C       [X+i]       ;то есть просто [X], здесь всегда i=0
131  A007              i       1   
132  8801  @@i_loop:   ZAcc    RoundZero
133  91EB              FMPM    [Z+3j+k]    ;то есть [Z+3j+2], он же "синус" (мнимая часть). т.к k=2 на всём протяжении
134  8AC4              C       [X+i]
135  92E3              FMA     [Z+3j+1]    ;он же "косинус" (действительная часть)
136  D480              [Y+i]   Acc
137  A859              iLOOP   @@i_loop    
138  CDDD              X       Y   ;возможно, X сохраним в стек, чтобы процедура минимально меняла регистры
139  80FC              Acc     [SP]
13A  81E7              PM      [Z+3j+i]    ;то есть [Z+3j], угол. 
13B  FC80              [SP]    Acc     
13C  A95A              jLOOP   @@j_loop
13D  8ED4              DIV2A   [Y+i]
13E  D480              [Y+i]   Acc
13F  8900      NOP  0 ;AUTOMATICALLY INSERTED BY TRANSLATOR TO PREVENT HAZARD
140  A7FF              ijk     [--SP]
141  EDFF              Z       [--SP]
142  CDFF              X       [--SP]
143  B0FF              JMP     [--SP]
atan1 endp


И у нас [X] = 32052 = 0x7D34, [X+1] = -6813 = 0xE563.

Поехали!

Сохраняем регистры в стек, инициализируем индексные регистры (все скопом), обнуляем аккумулятор, вычитаем из него [X] = 0x7D34. Результат [SP] инициализируем нулём, в регистр Z заносим адрес AtanTable, а в регистр Inv заносим знак: 1. Всё правильно, мы так и ожидали: здесь знак X будет обратным, но и число у нас лежит "-i", специально для этого.

Теперь заносим [X]=0x7D34 в регистр C, инициализируем i=1, опять обнуляем аккумулятор. И выполняем FMPM [Z+3j+k], причём загружается значение 0x8000. Результат пока "не виден" - умножение с накоплением занимает 18 тактов. А тем временем начинает выполняться C [X+i], значение 0xE563 поступает на шину данных, и затем мы ждём, пока освободится АЛУ - нельзя в него загружать C, пока оно предыдущую команду не завершило, иначе запорет!

Наконец, это дело завершено. Далее, FMA [Z+3j+1], и на шину поступает значение 0x0000. Тоже логично. И тут же мы останавливаемся, поскольку не можем вытащить из АЛУ результат, Acc (SrcAddr = 0x80), пока не завершатся вычисления.

Следующий слайд:


Всё верно: результатом вычислений стал 0x7D34, т.е действительно Inv=1 отработал как надо, превратив "-1" в "+1" (причём даже -32768 вполне себе превратилось в 32768, без "насыщения" на 32767, т.к аккумулятор у нас шире на 1 бит), именно так мы и ожидали.

И по-моему, я знаю, в чём подлянка! Программа тут не виновата, в самом алгоритме косяк :) Когда мы этот "вектор", (32052; -6813), повернули на 90° против часовой стрелки, у нас вышло (6813;32052), вот тут знак и сменился, про что я немножко забыл. В общем: комплексные числа у нас введены абсолютно верно, они явно ведут вектор в нужном направлении, а вот углы в таблице AtanTable все нужно взять со знаком "минус", и будет нам счастье. И последнее значение тоже добавлять надо со знаком "-", т.е заменить DIV2A на DIV2S.

Попробуем ещё разок:


И сведём в таблицу:


Да, примерно то же самое, но теперь знак правильный.

Ну и для очистки совести прошьём это безобразие на макет. Фиттеру в очередной раз "поплохело" и он завалил тайминги, всего 23,26 МГц. Но это же "наихудший случай", когда самый косячный кристалл, минимально-допустимое напряжение и максимальная рабочая температура. При комнатной температуре и номинальном напряжении должно по идее работать. Запустим:



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

РАБОТАЕТ!

Краткие итоги: 29 слов кода (58 байт) и 15 слов данных (30 байт). Точность: не хуже 2 единиц младшего разряда, преимущественно 0..1 единица. Время выполнения: около 21,6 мкс (540 тактов). Дополнительная функция: нахождение длины вектора с точностью не хуже 0,36%.

Да вообще-то, меня такое устраивает во всех отношениях :) Хотя есть некий спортивный интерес в получении полной точности, возможной в 16-битном представлении. В теории, алгоритм "atan на ЧЕТЫРЁХ умножениях" на это способен, ещё и уложится тактов в 100, но я не проверял ещё.

PS. Я неправильно ошибку считал: брал относительно исходных 12°, 24° и так далее, но ведь и представить их входными векторами точно я не мог. Сейчас это суммарная точность после двух операций. Возможно, если сравнивать с честно посчитанном atan(y/x) для двух целочисленных значений y,x, результат немного изменится в лучшую сторону. Ну и выборка пока маленькая.

PPS. Погоды сильно не делает. Там, где ошибка свыше 1,5 ЦМР, попробовал честно посчитать. Получил вместо 22 угл секунд ошибку в 21 угл секунду, ну да.
Tags: ПЛИС, математика, программки, работа, странные девайсы
Subscribe

  • Ремонт лыжных мостиков

    Вернулся с сегодняшнего субботника. Очень продуктивно: отремонтировали все ТРИ мостика! Правда, для этого надо было разделиться, благо народу…

  • Гетто-байк

    В субботу во время Великой Октябрьской резни бензопилой умудрился петуха сломать в велосипеде. По счастью, уже на следующий день удалось купить…

  • А всё-таки есть польза от ковариаций

    Вчера опробовал "сценарий", когда варьируем дальность от 1 метра до 11 метров. Получилось, что грамотное усреднение - это взять с огромными весами…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 0 comments