Sinus lookup tabel

Op 17 oktober 2020 16:50:46 schreef blackdog:
Het zal waarschijnlijk wel, maar in mijn brein rammeld het wat dat betreft, ben geen wiskunde wonder

En spelling blijkbaar ook niet

2 samples per sinus is genoeg, volgens Shannon.

Op 17 oktober 2020 18:01:31 schreef Frederick E. Terman:
Nou ja, als je iedere keer opnieuw de bijbehorende sinuswaarde uitrekent, ja. Maar hier willen we een beperkt aantal (64) waarden gebruiken, dus bij iedere klokpuls kun je niet een exacte waarde berekenen, maar moet je één van die 64 mogelijkheden kiezen.

Je kunt nog steeds linear interpoleren tussen twee tabelwaarden. Dat gaat verbazingwekkend goed, zeker als je je tabelwaarden optimaliseert voor de interpolatie.

Je bedoelt ze stiekem al wat hoger of lager maken zodat je de interpolatie een beetje kunt bijsturen en zo de 'vervorming' verdeelt?
Rekenkundig is het natuurlijk niet meer dan een add en een shift-right.

Op 18 oktober 2020 01:42:38 schreef blurp:
Je kunt nog steeds linear interpoleren tussen twee tabelwaarden. Dat gaat verbazingwekkend goed, zeker als je je tabelwaarden optimaliseert voor de interpolatie.

Maar dat zou je dan iedere keer moeten doen. Is het dan niet doelmatiger, gewoon je tabel groter te maken, dus toch maar 256 bytes gebruiken? Misschien is die uitbreiding nog wel minder dan je anders voor de interpolatie-code nodig zou hebben.

Tja, dat hangt natuurlijk van je applicatie en MCU af

Zoals ik al zij, ik heb ooit een DTMF generator gemaakt op basis van een 4bit AD converter, maar daar had ik wel een behoorlijk nauwkeurige samplefrequentie. (dat was in 1995, toen FPGA's nog nauwelijks bestonden en een ASIC nog "betaalbaar" was).

Dat werkte zonder analoog filter prima op het toenmalige PTT net.

Met een 4e orde-filter zou het me niets verbazen als je bruikbare DTMF kunt maken door gewoon de juiste frequentie als blok naar buiten te sturen.

@EricP: Ja, tabelwaarden wat hoger (lager onder nul) kiezen zodat je vervorming minimaliseert. Precieze waarden zou ik in een spreadsheetje met >10 keer oversampling uitproberen.

Update: Ik heb even met een spreadsheetje gespeeld. Met interpolatie tussen 4 waarden (0, 30, 60 en 90 graden) krijg ik een RMS-afwijking van 0.008. Ruimschoots geschikt voor DTMF dus.

[Bericht gewijzigd door blurp op maandag 19 oktober 2020 04:28:36 (12%)

edit: Ik zie FET me voor was met het noemen van dit algorithme. Het verschil: Ik optimaliseer door direct een macht van 2 voor de constante te gebruiken: een bitshift. Maar doe je het algemeen dan kan je inderdaad de periode mikken.... En: ik noem de stabiliteits-truuk die FET vergeet.

Op 17 oktober 2020 12:58:39 schreef meander:
Ik weet niet meer precies hoe t ging, de wiskunde erachter is lastig, maar je krijgt dan iets in de vorm van: y[n] = a . y[n-1] + b . y[n-2].

   x[n] = x[n-1] - c * y[n-1];
   y[n] = y[n-1] + c * x[n];

met c bijvoorbeeld 2^-8 .

Let op dat je in de tweede regel x[n-1] verwacht. Maar het is makkelijker te programmeren als je daar x[n] gebruikt (de zojuist veranderde x waarde!), bovendien (belangrijker) zorgt dit er voor dat de berekening stabiel is. Anders loopt ie weg.

Hierdoor wordt de code lekker simpel: Alle indexen kan je weglaten en je hoeft alleen x en y te onthouden.

De wiskunde is redelijk simpel: De afgeleide van sinus is cosinus, de afgeleide van cosinus is min sinus. In mijn voorbeeld moet je x[0] op 1 initializeren en y[0] = 0. Dan wordt x de cosinus en y de sinus.

Maar het probleem hiermee is dat je in iets van 2 * pi / c iteraties door de hele cyclus gaat. Dat is geen mooi getal. In het geval van Arco z'n probleem zou ie dus voor enige nauwkeurigheid inderdaad c=2^-8 moeten gebruiken, maar voor z'n 1463 Hz sinus moet ie dan voor ongeveer 100 iteraties door het algorithme heen.

Op 17 oktober 2020 19:30:48 schreef flipflop:
2 samples per sinus is genoeg, volgens Shannon.

Nee: BIJNA genoeg. Er is een < en geen <= ....

Als ik 1 * sin (pi*t) sample op iedere hele seconde krijg ik altijd nul. Dat is niet te onderscheiden van 0 * sin (maaktnietuit). Maar als ik weet dat 0.5 Hz de max freq in mijn signaal is, dan kan ik uit de samples van 1 * sin (0.999 * pi*t) het origineel reconstrueren.

als je daar x[n] gebruikt (de zojuist veranderde x waarde!),

Goed dat je dit noemt. Daar was ik inmiddels ook achter.. door het per ongeluk in een BASIC-loopje precies zo te programmeren, omdat ik vergat x even in een temp te stoppen.

Achteraf had ik het me moeten herinneren van toen ik eens een programmaatje had gemaakt om de bewegingen van hemellichamen te simuleren. Paar arrays voor positie, snelheid, versnelling, massa etc. van alle lichamen en simuleren maar. Dat komt in feite op hetzelfde neer, en inderdaad bleken banen stabiel te blijven met die truc (waarvan ik vermoed dat ik daar toen óók gewoon door slordigheid achter was gekomen).

(Ik maakte daarin trouwens nog een fraaie denkfout. Als eerste test: zon en aarde. Tik tik tik, save, RUN. Plaatje. Niets beweegt. Wat gaat er fout? Eerst koffie. Daarna vergeten.
Volgende dag: hé, aarde toch een stukje opgeschoven.
Tja, logisch: ik had voor ALLE input de werkelijke waarden genomen. En het zou dus een vol jaar duren voordat mijn aardetje rond het zonnetje gegaan zou zijn. )

64 staps tabel valt prima mee te leven...

Trapjes komen door scoop, polarisatievlek door fototoestel (bij normaal gebruik zie je geen vlekken)

Ziet er netjes uit zo.

Grote kans dat de scoop zélf maar 8 bits is . Nee, dat is het probleem niet. De samplerate van 10 kHz levert op het oog veel meer afwijking van de sinusvorm op (zie mijn plaatje hierboven), maar zelfs dat maakt, na filtering en 'gewoon even afwachten' niet zoveel uit.
In de sim (ideaal om in te knoeien en proberen) kom ik in dat geval op waarden voor de spurious van evengoed nog ruim beter dan 40 dB down.

/OT hoe werkt dat met die polarisatievlek? Ik weet hoe het display werkt, maar wat doet het fototoestel daar dan mee?

[Bericht gewijzigd door Frederick E. Terman op maandag 19 oktober 2020 17:33:18 (13%)

Het LCD heeft een filter, en de camera heeft ook een coating. Waarschijnlijk kunnen die twee niet samen door één deur...
Mijn oude Philips heeft geen last van polarisatie...

[Bericht gewijzigd door Arco op maandag 19 oktober 2020 20:37:19 (23%)

D'r gaat niks boven een klassiek elektronenbombardement op een ordentelijke fosfor.

Op 17 oktober 2020 12:58:39 schreef meander:
Ik weet niet meer precies hoe t ging, de wiskunde erachter is lastig, maar je krijgt dan iets in de vorm van: y[n] = a . y[n-1] + b . y[n-2].

Op 19 oktober 2020 08:59:57 schreef rew:

code:

   x[n] = x[n-1] - c * y[n-1];
   y[n] = y[n-1] + c * x[n];

met c bijvoorbeeld 2^-8 .

M.i. bedoelt meander iets met uitsluitend historische y-waarden.
Het recept hierboven is dan equivalent met
y[n] = (2 - c^2)*y[n-1] - y[n - 2] voor n >= 2, y[0] = 0, y[1] = c.
Dit geeft exact een reeks waarden van sin(t), maar de tijdstap is iets groter dan c en de amplitude iets groter dan 1. De afwijkingen worden voor kleine c (meer stappen per periode) snel verwaarloosbaar. Bijv. voor c = 0.1 (zowat 64 stappen per periode) resp. 0.042% en 0.125% .

Op 19 oktober 2020 21:44:42 schreef Frederick E. Terman:
D'r gaat niks boven een klassiek elektronenbombardement op een ordentelijke fosfor.

Inderdaad. .

@arco.

Zwaar jaloers op je. De pm331X vindt ik de gaafste scoop ooit.

Op 20 oktober 2020 19:57:30 schreef aobp11:
[...]
[...]
M.i. bedoelt meander iets met uitsluitend historische y-waarden.
Het recept hierboven is dan equivalent met
y[n] = (2 - c^2)*y[n-1] - y[n - 2] voor n >= 2, y[0] = 0, y[1] = c.
Dit geeft exact een reeks waarden van sin(t), maar de tijdstap is iets groter dan c en de amplitude iets groter dan 1. De afwijkingen worden voor kleine c (meer stappen per periode) snel verwaarloosbaar. Bijv. voor c = 0.1 (zowat 64 stappen per periode) resp. 0.042% en 0.125% .

Ik heb het nog eens opgezocht, zie de afbeelding hieronder geknipt uit een dictaat wat ik nog ergens had liggen. Ik heb de hele wiskundige afleiding weggelaten, maar het komt er op neer dat je een sinus signaal kan omschrijven als een complexe vector en dat helemaal kan uitwerken naar die notatie.

y[n-1] is dan de laatste waarde die je uitgestuurd hebt, y[n-2] de waarde daarvoor. bij het opstarten heb je deze getallen natuurlijk nog niet, dus wordt gesteld dat y[n-1] = 1 en y[n-2] = 0; als je dit vervolgens door laat lopen komt er automatisch een sinus uit.

de getallen voor a1 bepalen de frequentie. om dit uit te rekenen heb je ook een cosinus nodig, maar als je slechts een beperkte hoeveelheid tonen wilt uitsturen, kun je deze met de hand al uitrekenen en als constante programmeren.

Op 21 oktober 2020 11:06:01 schreef High met Henk:
@arco.

Zwaar jaloers op je. De pm331X vindt ik de gaafste scoop ooit.

Ik kan zelfs nog kiezen. Er zit nu een bruin front op, maar ik heb ook nog een grijze... (was in beide kleuren leverbaar)

@meander,
De afleiding in het dictaat komt waarschijnlijk neer op het nog eens afleiden van gonioformule
sin(a + b) + sin(a - b) = 2 sin(a)cos(b).
Wellicht bekend van de zijbanden bij amplitudemodulatie.
Met de hoekstap α is y[n] = sin(nα) ( bij schaalfactor 1).
Met a = nα en b = α komt er
y[n - 1] + y[n + 1] = 2 y[n] cos(α)
en heb je de recurrente betrekking al.
Ik zie net dat in het diktaat ten onrechte een - teken voor 2 cos(ω*Tsample). Later a1 = 1.96 wel in orde.
Om piekwaarde 1 te krijgen moet je y[1] = sin(α) nemen, niet y[1] = 1.

Wie snapt de afrondingsregels van C goed genoeg om te verklaren dat de resultaten verschillend zijn tussen de twee versies?

#include <stdio.h>
int main (int argc, char **argv)
{
  int x,y;
  double c;

  c = 1/256.0;

  x = 65000;
  y = 0;

  while (1) {
#if 0
    x -= y * c;
    y += x * c; 
#else
    x -= y >> 8;
    y += x >> 8; 
#endif
    printf ("%d %d\n", x, y);
  }
}

Edit: Als y negatief is, verschilt y >> 8 en y*c bijna altijd precies 1.

OK. Versie voor 1444 Hz met samples op 10kHz. (die 1444 had 1477 moeten zijn. oefening voor de lezer om dat te fixen).

#include <stdio.h>
#include <math.h>
int main (int argc, char **argv)
{
  int x,y;
  double c;

  //  c = 1/256.0;
  c = 3.1415926535 * 2 *1444/10000;

  x = 65000;
  y = 0;

  while (1) {
    x -= floor (y * c);
    y += floor (x * c); 
    printf ("%d %d\n", x, y);
  }
}

Ik krijg dan dit:

Resteert om c te schrijven als 29730 / 32768 . Dan staat er

    x -= floor (y * c);

en maak je er van:

    x -= (y * 29730) >> 15;

Dan moet je de berekening rond y hier wel in 32 bits doen. Ohja, Kennelijk krijg ik wat overshoot, Die 65000 om bijna 16 hele bitjes te benutten is zowiezo al 1 bitje te veel, maar met zo weinig stappen schijn ik ook tussenresultaten te krijgen die groter zijn dan de initiele waarde. Dus kennelijk wordt ie 12% groter dus de initializatie moet dan op 2^15 / 1.14 o.i.d. worden.

Op 21 oktober 2020 16:51:08 schreef meander:
[...] geknipt uit een dictaat wat ik nog ergens had liggen.

Dat is een leuke! Met y(0) bepaal je de fase, en y(1) geef je de golfvorm de juiste begin-helling om op de gewenste amplitude uit te komen. Wat aardig.
Hoe heet het dictaat? Misschien vind ik het ook nog ergens.

Op 21 oktober 2020 20:16:04 schreef Frederick E. Terman:Hoe heet het dictaat? Misschien vind ik het ook nog ergens.

Ik denk het niet. Het was ooit geschreven door een docent op het HBO. het is een stukje theorie verpakt in een programmeeropdracht om een microcontroller für Elise te laten spelen met sinusvormige signalen.
Ik weet niet hoe leuk ze het vinden als ik dat hele ding hier online gooi (met name die docent was nogal gehecht aan intellectueel eigendom volgens mij), maar wellicht kan ik het je toesturen.

Op 21 oktober 2020 20:16:04 schreef Frederick E. Terman:
[...]Dat is een leuke! Met y(0) bepaal je de fase, en y(1) geef je de golfvorm de juiste begin-helling om op de gewenste amplitude uit te komen. Wat aardig.

Zo simpel ligt het niet. y[0] en y[1] bepalen samen welke sinus eruit komt. y[0] alleen legt de fase nog niet vast. Je zult van de gewenste sinuslijn zowel y[0] als y[1] moeten aflezen.

Ik heb trouwens zo'n bruin vermoeden dat het schema met (x[n],y[n]) minder last heeft van rekenonauwkeurigheid dan het schema met alleen y[n]. Behalve misschien als je
y[n + 1] = (2 - c^2)*y[n] - y[n-1]
uitrekent als
y[n + 1] = y[n] + { (y[n] - y[n-1]) - c^2*y[n]}.

Ik zei het inderdaad onnauwkeurig; de fase en amplitude bepalen y(0). Maar natuurlijk ook y(1).
In elk geval: als je een bepaalde amplitude en fase wilt hebben, hoef je maar op die twee plaatsen de daarbij horende waarden in te vullen, en dan loopt het vanzelf goed door. Dat vond ik wel geinig.

Als ik in Excel een miljoen punten laat berekenen, vallen de perioden nog steeds precies over elkaar - de onnauwkeurigheden blijven dus gelijk en op- dezelfde plaats. Beslist niet slechter dan de x,y-machine lijkt het dus, althans de paar waarden die ik geprobeerd heb.

@meander: Nou, graag! Al zal ik zonder ook beslist wel overleven, het is wél erg fraai.

Nog even: Mijn arm-compiler maak van de loop:

.L2:
        mul     r3, r6, r5
        sub     r4, r4, r3, asr #15
        mul     r3, r6, r4
        add     r5, r5, r3, asr #15

maar dit is dus met de frequentie-constante preloaded in het R6 register. Toen ik met -O2 compileerde (meer op snelheid) ging ie slim doen om de mul instructie te omzeilen: drie shift-en-add dingen om de "langzame" mul te voorkomen.

In mijn cortex-m0 handleiding staat echter:

providing high-end processing hardware including a single-cycle multiplier.

Dus dit zal zo wel sneller zijn. Nu met -Os gecompileerd.

Op 17 oktober 2020 19:30:48 schreef flipflop:
2 samples per sinus is genoeg, volgens Shannon.

Klopt, dat noemen we een blokgolf... en ja, als je die perfect filtert zodat de eerste harmonische afdoende onderdrukt is dan is het een prachtige sinus. Blijft wel dat als je de frequentie een klein beetje aanpast je extreem veel last van aliasing krijgt. Idealiter wil je dus samples hebben zodanig dat het een veelvoud is van ieder van je frequenties, dus in het slechste geval een verdubbeling van het aantal samples per frequentie die je wilt kunnen genereren.
Daarom is het redelijker om een sinus te creeëren die voor alle gewenstie frequenties min of meer in orde is, zodanig dat het filter de eerste harmonische van de hoogste frequentie ie iedere situatie afdoende onderdrukt.