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

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

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

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 со знаком минус.

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

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

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

Код:

---

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-ти.
Знаю, что инверсия и отрицание отличаются. Тут это большой роли не сыграло. Основная погрешность, всё-таки, из-за небольшого фазового расхождения.

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

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

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

Цитата:

---

Сообщение от Reobne

Это я теорию функций комплексного переменного вспомнил.
Если без неё, то примерно так:
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) единицу, затухать перестало)
Может быть не хватает разрядности бейсика?)

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

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

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

Цитата:

---

Сообщение от Reobne

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

---

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

Код:

---

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

---

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

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

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