Генерация синуса

[Titus]

Ладно, раз всем лениво, раскрою алгоритм, который весьма прост)

В данном случае используется алгоритм генерации синуса посредством рекурсивного цифрового генератора (о чем упоминалось даже в этой теме неоднократно). Многим он больше знаком по алгоритму Герцеля:

D0 = a1 * D1 - D2 + x,
D2 = D1,
D1 = D0,

который является модификацией данного генератора с добавленным возбуждением коэффициентом x. А сам генератор в чистом виде - это:

D0 = a1 * D1 - D2,
D2 = D1,
D1 = D0,
где a1 = 2 * cos (2 * pi * частота / частота_отсчетов),
или же a1 = 2 * cos (2 * pi * k);

Ввиду того, что алгоритм рекурсивный, требуется большая точность (разрядность) чисел. Приемлимой точности даже 16 битами (не говоря уже о 8) на отрезке в пару сот итераций не добиться. Поэтому выбрана разрядность 24 бита.

Далее о коэффициенте a1. Понятно, что ни о каком умножении речи быть не может. Единственный вариант - представить a1 в виде достаточно простого полинома вида 2 * (2^0 + 2^(...) + 2^(...) + ...) и т.д.
Самым простым полиномом, удовлетворяющим периоду близкому к 256, и при этом имеющим близкий к целому числу корень 1/k - это полином 2 * (2^0 - 2^(-11)), что легко представлятся в виде сдвига вправо на 8 + 3 бит, и одного вычетания. Что я и изобразил в своем алгоритме.

Посчитаем 1/k (число отсчетов) для данного полинома.

a1 = 2 * (1 - 2^(-11)) = 1,9990234375

1 / k = 1 / (arccos(a1 / 2) / (2 * pi) = 201.0536,

что является достаточно точным, близким к целочисленному периодом.

Если кого интересует полином для периода 256, то я тоже прикидывал его, у меня получился 2 * (2^0 - 2^(-12) - 2^(-14) + 2^(-18)), что потребует уже целых 3 сдвига и 3 сложения, усложняя код, ведь это все надо в 24-х битной арифметике сдвигать и складывать (ну, младшие части полинома можно и в 16, и даже в 12-битной).

Посчитаем 1/k (число отсчетов) для данного полинома:

a1 = 2 * (1 - 2^(-12) - 2^(-14) + 2^(-18)) = 1,99939727783203125

1 / k = 1 / (arccos(a1 / 2) / (2 * pi) = 255.92,

что уже менее точно, но все же достаточно приемлимо.

Вот вышеприведенный уже код, но с комментариями:

	org	$8000

	LD	B,201		; Длина таблицы синуса 201 байт
	LD	HL,$8100	; Адрес таблицы синуса
	EXX			;
	
	XOR	A		; D0/D1 = $00.00.00
	LD	H,A		;
	LD	L,A		;
	
	PUSH	HL		; D2 = -$02.00.00
        LD      IXL,-2          ;
	
	JR	LoopIn		; --> LoopIn
	
;------------------------------- Главный цикл
MLoop:	
	EXX			;
	
	SBC	HL,BC		; A.H.L = A.H.L - LX.B.C
	SBC	A,IXL		; (D0 = a1*D1 - D2)

        LD      IXL,D            ; Сохранить LX для D2 будущей итерации

;-------------------------------
LoopIn:				
				; D0/D1 = A.H.L
				
	LD	D,A		; D = A (и заодно сохранить A)
	ADD	A,A		;
	SBC	A,A		;
	LD	IXH,A		; A = $00/$FF, в зависимости от знака
	LD	A,D		; Восстановить A
	LD	E,H		;
        LD      B,3             ; DE = (A.H.L >> 8) >> 3 = A.H.L >> 11
RLoop:	SRA	D		;     
	RR	E		;
	DJNZ	RLoop		;
	
	POP	BC		; Восстановить D2 = LX.B.C
	PUSH	HL		; Сохранить D1 = (A).H.L
	
	SBC	HL,DE		; A.H.L = A.H.L - HX.D.E 
        LD      D,A             ;
	SBC	A,IXH		;
	
	ADD	HL,HL		; A.H.L = A.H.L << 1
	ADC	A,A		; (D0 = a1 * D1)
	
	EXX			; Сохранить A в (HL)+
	LD	(HL),A		; и цикл
	INC	L		;
	DJNZ	MLoop		;
	
        POP     AF              ; Восстановить стек
       
	RET			; Выход

---------- Post added at 15:21 ---------- Previous post was at 15:19 ----------

Да, начальная амплитуда задается через коэффициент в D2 со знаком минус.

[Titus]

Если хотите, программа на бейсике, чтобы побаловаться с коэффициентом a1 и периодом. Специально сделан 10-кратный прогон, чтобы было видно, как периоды ложаться друг на друга, т.е. чтобы сразу была видна погрешность.

[Titus]

Кстати, можно написать программу поиска близких к целочисленным корней 1/k для коротких полиномов (ну, скажем до 4-й степени), которая простым перебором найдет решения. Если они есть)

[Reobne]

Вот что получилось.
Алгоритм, немного переделал, но смысл тот-же:

Real x=0,y=1
const d=2*PI/N/M, si=sin(d), co=cos(d)

for i=0 to N-1
  for j=1 to M
    x = x*co+y*si
    y = -x*si+y*co
  next M
  SinTab[i]=x
next N

Разрядность вычислений принял 16.
N я принял за 256, а M подбирал, чтобы si и co получились удобные для умножения. В итоге принял М за 5.
При этом константу "co" можно вообще выкинуть, она, для 16 бит неотличима от единицы.
si=#0.014
То есть сдвинуть на 8 разрядов, сумма, потом ещё 2 разряда, и ещё сумма.
В итоге получилась такой фрагмент кода:

sintabgen5
  ld de,289
  ld hl,32769
  xor a
  exx
  ld hl,sintab+#80
  ld d,h
  ld e,l
stb_loop1
  neg
  ld (de),a
  ld b,10
stb_loop2
  ld a,b
  rra
  sbc a,a
  exx
  ex de,hl
  xor d
  ld c,a
  rla
  sbc a,a
  ld b,a
  add hl,bc
  sra c
  sra c
  add hl,bc
  ld a,d
  exx
  djnz stb_loop2
  inc e
  dec l
  ld (hl),a
  jr nz,stb_loop1

Занимает она 43 байта. Погрешность около 2.
Работает, конечно, медленнее, есть внутренний цикл, но это меня совершенно не тревожит.
XOR там, чтобы знак менялся. По формуле надо то плюс, то минус. Заметили, что EX DE,HL только одно? Поэтому внутренний цикл до 10-ти.
Знаю, что инверсия и отрицание отличаются. Тут это большой роли не сыграло. Основная погрешность, всё-таки, из-за небольшого фазового расхождения.

Я конечно старался кодить покороче, до рекорда немного не дотянул. Интересно, совсем другой подход, а размер получился похожий. Прям закон информатики!
Но потенциал оптимизации, как всегда, неизвестен! А может даже точность повысить получится. Неизвестно.

[drbars]

Все эти рекорды ни к чему, если процедура не универсальна. Кто напишет универсальный код, чтобы на входе были параметры периода и амплитуды — тот и герой Коротко, быстро и точно.

А тот, кто демо в 256 байт пишет синус из калькулятора возьмёт. Точно и сердито

[Titus]

Это я теорию функций комплексного переменного вспомнил.
Если без неё, то примерно так:
Real x=0,y=1
for i=0 to 200
x, y = x*cos(d)+y*sin(d), -x*sin(d)+y*cos(d)
SinTab[i]=x
next


d=PI/100

Ты уверен в этой формуле? Написал на бейсике, дает затухающие колебания, каждый следующий период меньше предыдущего)

---------- Post added at 23:08 ---------- Previous post was at 23:06 ----------

Сделал вместо cos(d) единицу, затухать перестало)
Может быть не хватает разрядности бейсика?)

[Reobne]

В формуле уверен. Это поворот вектора на угол.

Шлю свой бейсик-пример, сравни.

[Reobne]

Исправь там 200 на 2000. Увидишь, что 10 кругов пиксель в пиксель легли.

[Titus]

Исправь там 200 на 2000. Увидишь, что 10 кругов пиксель в пиксель легли.

Понятно почему)
Я это:

   x = x*co+y*si
    y = -x*si+y*co

В лоб считал, а надо сперва t = x*co + y*si, потом уже x=t;

Ну не пользуюсь я в повседневной жизни, никаких своих программах, а так же при походе на вершину Джомолунгмы - комплексными числами)

[Reobne]

Чем меньше шаг угла, тем меньше погрешности от такого упрощения. Поэтому я его и применил, у меня угол-то не ПИ/100, а ПИ/(128*5)=ПИ/640.
Надо было сразу об этом нюансе написать, не подумал.