Sinussignal in C erzeugen, ohne die Standardfunktion zu verwenden

Question 1

Ich möchte ein Sinussignal in C erzeugen, ohne die Standardfunktion sin() zu verwenden, um sinusförmige Änderungen der Helligkeit einer LED auszulösen. Meine Grundidee war, eine Lookup-Tabelle mit 40 Punkten und Interpolation zu verwenden.

Hier mein erster Ansatz:

const int sine_table[40] = {0, 5125, 10125, 14876, 19260, 23170, 26509, 29196,
31163, 32364, 32767,  32364, 31163, 29196, 26509, 23170, 19260, 14876, 10125,
5125, 0, -5126, -10126,-14877, -19261, -23171, -26510, -29197, -31164, -32365,
-32768, -32365, -31164, -29197, -26510, -23171, -19261, -14877, -10126, -5126};

int i = 0;
int x1 = 0;
int x2 = 0;
float y = 0;

float sin1(float phase)
{
    x1 = (int) phase % 41;
    x2 = x1 + 1;
    y = (sine_table[x2] - sine_table[x1])*((float) ((int) (40*0.001*i*100) % 4100)/100 - x1) + sine_table[x1];
    return y;
}

int main()
{
    while(1)
    {
    printf("%f      ", sin1(40*0.001*i)/32768);
    i = i + 1;
    }
}

Leider gibt diese Funktion manchmal Werte weit größer als 1 zurück. Außerdem scheint die Interpolation nicht gut zu sein (ich habe damit sinusförmige Helligkeitsänderungen einer LED erzeugt, aber diese sind sehr unstetig).

Hat jemand eine bessere Idee, einen Sinusgenerator in C zu implementieren?

Question 2

Das Hauptproblem von OP besteht darin, den Index für die Tabellensuche zu generieren.

Der Code von OP versucht, auf ein externes Array zuzugreifen sine_table[40] führt zu undefiniertes Verhalten. Repariere das wenigstens.

const int sine_table[40] = {0, 5125, 10125, ...
    ...
    x1 = (int) phase % 41;                     // -40 <= x1 <= 40
    x2 = x1 + 1;                               // -39 <= x2 <= 41  
    y = (sine_table[x2] - sine_table[x1])*...  // bad code, consider x1 = 40 or x2 = 40,41

Vorgeschlagene Änderung

    x1 = (int) phase % 40;   // mod 40, not 41
    if (x1 < 0) x1 += 40;    // Handle negative values
    x2 = (x1 + 1) % 40;      // Handle wrap-around 
    y = (sine_table[x2] - sine_table[x1])*...

Es gibt viel bessere Ansätze, aber um sich auf die Methode von OP zu konzentrieren, siehe unten.

#include <math.h>
#include <stdio.h>

const int sine_table[40] = { 0, 5125, 10125, 14876, 19260, 23170, 26509, 29196,
31163, 32364, 32767, 32364, 31163, 29196, 26509, 23170, 19260, 14876, 10125,
5125, 0, -5126, -10126, -14877, -19261, -23171, -26510, -29197, -31164, -32365,
-32768, -32365, -31164, -29197, -26510, -23171, -19261, -14877, -10126, -5126 };

int i = 0;
int x1 = 0;
int x2 = 0;
float y = 0;

float sin1(float phase) {
  x1 = (int) phase % 40;
  if (x1 < 0) x1 += 40;
  x2 = (x1 + 1) % 40;
  y = (sine_table[x2] - sine_table[x1])
      * ((float) ((int) (40 * 0.001 * i * 100) % 4100) / 100 - x1)
      + sine_table[x1];
  return y;
}

int main(void) {
  double pi = 3.1415926535897932384626433832795;
  for (int j = 0; j < 1000; j++) {
    float x = 40 * 0.001 * i;
    float radians = x * 2 * pi / 40;
    printf("%f %f %f\n", x, sin1(x) / 32768, sin(radians));
    i = i + 1;
  }
}

Ausgabe

         OP's     Reference sin()
0.000000 0.000000 0.000000
0.040000 0.006256 0.006283
0.080000 0.012512 0.012566
...
1.960000 0.301361 0.303035
2.000000 0.308990 0.309017
2.040000 0.314790 0.314987
...
39.880001 -0.020336 -0.018848
39.919998 -0.014079 -0.012567
39.959999 -0.006257 -0.006283

Besserer Code würde die Werte nicht weitergeben i, x1, x2, y als globale Variablen, sondern als Funktionsparameter oder Funktionsvariablen. Vielleicht ist das ein Artefakt des Debugging von OP.

Hat jemand eine bessere Idee, einen Sinusgenerator in C zu implementieren?

Das ist ziemlich breit. Besser in Bezug auf Geschwindigkeit, Präzision, Coderaum, Portabilität oder Wartbarkeit? sine() Funktionen sind einfach zu erstellen. Hochwertigere erfordern mehr Aufwand.

Obwohl unscharf, ist die Verwendung einer kleinen Nachschlagetabelle durch OP ein guter Anfang – obwohl ich sehe, dass dies ohne Gleitkomma-Mathematik möglich ist. Ich empfehle OP, eine getestete und funktionierende Lösung zu erstellen und zu veröffentlichen Code-Review für Verbesserungsideen.

Question 3

…eine bessere Idee, einen Sinusgenerator in C zu implementieren?

Edit: Erste Lektüre vorschlagen Dieser Artikel um ein Verständnis dafür zu bekommen, was OP verlangt.

Aus dem in Ihrer Frage angegebenen Kontext “besser” wahrscheinlich bedeutet Größe und/oder Geschwindigkeit von kompiliertem Code, vielleicht um einen kleinen Mikroprozessor zu unterstützen.

Das CORDIC ( COOrdinate RDrehung DIGital CComputer )-Algorithmus eignet sich sehr gut für die Verwendung auf kleineren uP- und FPGA-Implementierungen, die über begrenzte mathematische Berechnungsmöglichkeiten verfügen berechnet den Sinus und Cosinus eines Wertes nur mit Grundrechenarten (Addition, Subtraktion und Verschiebungen). Mehr über CORDIC und wie man damit Sinus/Cosinus eines Winkels erzeugt werden hier bereitgestellt.

Es gibt auch mehrere Websites, die Beispiele für die Implementierung von Algorithmen bereitstellen. Einfach CORDIC enthält detaillierte Erläuterungen zum Generieren einer Tabelle, die dann für die Verwendung auf Ihrem Zielgerät vorkompiliert werden kann, sowie Code zum Testen der Ausgabe der folgenden Funktion (die Festkomma-Mathematik verwendet):

(Siehe Dokumentation der folgenden und andere Funktionen im Link)

#define cordic_1K 0x26DD3B6A
#define half_pi 0x6487ED51
#define MUL 1073741824.000000
#define CORDIC_NTAB 32
int cordic_ctab [] = {0x3243F6A8, 0x1DAC6705, 0x0FADBAFC, 0x07F56EA6, 0x03FEAB76, 0x01FFD55B, 
0x00FFFAAA, 0x007FFF55, 0x003FFFEA, 0x001FFFFD, 0x000FFFFF, 0x0007FFFF, 0x0003FFFF, 
0x0001FFFF, 0x0000FFFF, 0x00007FFF, 0x00003FFF, 0x00001FFF, 0x00000FFF, 0x000007FF, 
0x000003FF, 0x000001FF, 0x000000FF, 0x0000007F, 0x0000003F, 0x0000001F, 0x0000000F, 
0x00000008, 0x00000004, 0x00000002, 0x00000001, 0x00000000 };

void cordic(int theta, int *s, int *c, int n)
{
  int k, d, tx, ty, tz;
  int x=cordic_1K,y=0,z=theta;
  n = (n>CORDIC_NTAB) ? CORDIC_NTAB : n;
  for (k=0; k<n; ++k)
  {
    d = z>>31;
    //get sign. for other architectures, you might want to use the more portable version
    //d = z>=0 ? 0 : -1;
    tx = x - (((y>>k) ^ d) - d);
    ty = y + (((x>>k) ^ d) - d);
    tz = z - ((cordic_ctab[k] ^ d) - d);
    x = tx; y = ty; z = tz;
  }  
 *c = x; *s = y;
}

Bearbeiten:

Ich habe die Dokumentation zur Verwendung der Beispiele unter gefunden Einfach CORDIC Website sehr einfach zu folgen. Eine kleine Sache, auf die ich gestoßen bin, war jedoch das Kompilieren der Datei cordic-test.c der Fehler: Verwendung des nicht deklarierten Bezeichners ‘M_PI’ aufgetreten. Es scheint, dass beim Ausführen des kompilierten gentable.c Datei (die die cordic-test.c Datei) die Zeile:

#define M_PI 3.1415926535897932384626

Obwohl es in seinen eigenen Deklarationen enthalten war, war es nicht in den printf-Anweisungen enthalten, die zum Erstellen der Datei verwendet wurden cordic-test.c. Nachdem dies behoben war, funktionierte alles wie angekündigt.

Wie dokumentiert, erzeugt der Bereich der erzeugten Daten 1/4 eines vollständigen Sinuszyklus (-π/2 – π/2). Die folgende Abbildung enthält eine Darstellung der tatsächlichen Daten, die zwischen den hellblauen Punkten erzeugt werden. Der Rest des Sinussignals wird durch Spiegeln und Transponieren des ursprünglichen Datenabschnitts hergestellt.

Question 4

Das Erzeugen einer genauen Sinusfunktion erfordert eine Menge an Ressourcen (CPU-Zyklen und Speicher), die in dieser Anwendung nicht gerechtfertigt ist. Ihr Ziel, eine “glatte” Sinuskurve zu erzeugen, geht an den Anforderungen der Anwendung vorbei.

Während Sie beim Zeichnen der Kurve möglicherweise Unvollkommenheiten bemerken, nimmt das menschliche Auge diese Unvollkommenheiten überhaupt nicht wahr, wenn Sie diese Kurve auf einen LED-PWM-Antrieb anwenden.
Auch bei einer 40-stufigen Kurve wird das menschliche Auge den Helligkeitsunterschied zwischen benachbarten Werten wahrscheinlich nicht wahrnehmen, sodass eine Interpolation nicht erforderlich ist.

Im Allgemeinen ist es effizienter, wenn Sie eine Sinusfunktion generieren, die die entsprechenden PWM-Antriebswerte direkt ohne Gleitkomma erzeugt. Tatsächlich eher als eine Sinusfunktion a skalierter erhöhter Kosinus sinnvoller, so dass eine Eingabe von Null zu einer Ausgabe von Null führt und eine Eingabe der halben Anzahl von Werten im Zyklus den Maximalwert für Ihren PWM-Antrieb ergibt.

Die folgende Funktion generiert eine Raised-Cosinus-Kurve für eine 8-Bit-FSD-PWM aus einer Suche mit 16 Werten (und 16 Bytes), die einen Zyklus mit 59 Schritten erzeugt. Es ist also sowohl speicher- als auch leistungseffizient im Vergleich zu Ihrer 40-Schritt-Gleitkommaimplementierung.

#include <stdint.h>

#define LOOKUP_SIZE 16
#define PWM_THETA_MAX (LOOKUP_SIZE * 4 - 4)

uint8_t RaisedCosine8bit( unsigned n )
{
    static const uint8_t lookup[LOOKUP_SIZE] = { 0, 1, 5, 9,
                                                 14, 21, 28, 36,
                                                 46, 56, 67, 78,
                                                 90, 102, 114, 127} ;
    uint8_t s = 0 ;
    n = n % PWM_THETA_MAX ;

    if( n < LOOKUP_SIZE )
    {
        s = lookup[n] ;
    }
    else if( n < LOOKUP_SIZE * 2 - 1 )
    {
        s = 255 - lookup[LOOKUP_SIZE * 2 - n - 2] ;
    }
    else if( n < LOOKUP_SIZE * 3 - 2 )
    {
        s = 255 - lookup[n - LOOKUP_SIZE * 2 + 2] ;
    }
    else
    {
        s = lookup[LOOKUP_SIZE * 4 - n - 4] ;
    }

    return s ;
}

Für eine Eingabe von 0 <= theta < PWM_THETA_MAX die kurve sieht so aus:

Was ich vorschlage, ist viel glatt genug für die Beleuchtung.

In der Praxis könnten Sie es so verwenden:

for(;;)
{
    for( unsigned i = 0; i < PWM_THETA_MAX; i++ )
    {
        LedPwmDrive( RaisedCosine8bit( i ) ) ;
        Delay( LED_UPDATE_DLEAY ) ;
    }
}

Wenn Ihr PWM-Bereich nicht 0 bis 255 beträgt, skalieren Sie einfach die Ausgabe der Funktion; 8-Bit-Auflösung ist mehr als genug für die Aufgabe.

Question 5

Der klassische Trick, einen Kreis zu zeichnen (und damit auch eine Sinuswelle zu erzeugen), ist Hakmem #149 von Marvin Minsky. Z.B,:

#include <stdio.h>

int main(void)
{
    float x = 1, y = 0;

    const float e = .04;

    for (int i = 0; i < 100; ++i)
    {
        x -= e*y;
        y += e*x;
        printf("%g\n", y);
    }
}

Es wird leicht exzentrisch sein, kein perfekter Kreis, und Sie können einige Werte leicht über 1 erhalten, aber Sie können dies anpassen, indem Sie durch das Maximum dividieren oder runden. Außerdem kann ganzzahlige Arithmetik verwendet werden, und Sie können Multiplikation/Division eliminieren, indem Sie eine negative Zweierpotenz für verwenden ealso kann stattdessen shift verwendet werden.

Question 6

Haben Sie sich überlegt, den Anteil der Sinuskurve aus zu modellieren [0..PI] als Parabel? Wenn die Helligkeit der LED nur von einem menschlichen Auge beobachtet werden soll, sollten die Formen der Kurven ähnlich genug sein, so dass ein geringer Unterschied erkannt wird.

Sie müssten nur die entsprechende Gleichung finden, um sie zu beschreiben.

Hmmm, …

Scheitelpunkt bei (PI/2, 1)

Schnittpunkte der X-Achse bei (0, 0) und (PI, 0)

f(x) = 1 - K * (x - PI/2) * (x - PI/2)

Wo K wäre …

K = 4 / (PI * PI)

Question 7

Für eine LED könnten Sie wahrscheinlich nur mit etwa 16 Schritten auskommen, ohne überhaupt zu interpolieren. Das heißt, ich kann mindestens zwei merkwürdige Dinge an dir erkennen sin1() Funktion:

1) Sie haben 40 Datenpunkte drin sine_tableaber Sie nehmen den Index x1 Modulo 41 der Eingabe. Das scheint nicht der richtige Weg zu sein, um mit der Periodizität umzugehen, und lasst uns x1 zeigen hinter den letzten Index des Arrays.
2) Sie addieren dann +1, also x2 kann noch weiter über die Grenzen des Arrays hinausgehen.
3) Sie verwenden i in der Funktion, wird aber nur im Hauptprogramm gesetzt. Ich kann nicht sagen, was es tun soll, aber die Verwendung eines solchen Global in einer einfachen Berechnungsfunktion scheint zumindest schmutzig zu sein. Vielleicht soll es den Bruchteil für die Interpolation liefern, aber sollten Sie nicht verwenden phase dafür.

Hier ist ein einfacher Interpolator, der zu funktionieren scheint. Nach Geschmack anpassen.

#include <assert.h>

int A[4] = {100, 200, 400, 800};    
int interpolate(float x)
{
    if (x == 3.00) {
        return A[3];
    }
    if (x > 3) {
        return interpolate(6 - x);
    }
    assert(x >= 0 && x < 3);
    int i = x;
    float frac = x - i;
    return A[i] + frac * (A[i+1] - A[i]);
}

Einige willkürliche Beispielausgaben:

interpolate(0.000000) = 100
interpolate(0.250000) = 125
interpolate(0.500000) = 150
interpolate(1.000000) = 200
interpolate(1.500000) = 300
interpolate(2.250000) = 500
interpolate(2.999900) = 799
interpolate(3.000000) = 800
interpolate(3.750000) = 500

_{(Ich überlasse es dem interessierten Leser, alle Vorkommen von zu ersetzen 3 mit einer richtig definierten symbolischen Konstante, um die Funktion weiter zu verallgemeinern und auch die Berechnung der negativen Phase zu implementieren.)}

Question 8

Ich würde mit der Annäherung von Bhaskara I an eine Sinusfunktion gehen. Mit Grad von 0 bis 180 können Sie den Wert so annähern

float Sine0to180(float phase)
{
    return (4.0f * phase) * (180.0f - phase) / (40500.0f - phase * (180.0f - phase));
}

Wenn Sie einen Winkel berücksichtigen möchten, würden Sie hinzufügen

float sine(float phase)
{
    float FactorFor180to360 = -1 * (((int) phase / 180) % 2 );
    float AbsoluteSineValue = Sine0to180(phase - (float)(180 * (int)(phase/180)));
    return AbsoluteSineValue * FactorFor180to360;
}

Wenn Sie es im Bogenmaß tun möchten, würden Sie hinzufügen

float SineRads(float phase)
{
    return Sine(phase * 180.0f / 3.1416);
}

Hier ist ein Diagramm, das die mit dieser Annäherung berechneten Punkte und auch die mit der Sinusfunktion berechneten Punkte zeigt. Sie können die Annäherungspunkte kaum sehen, die unter den tatsächlichen Sinuspunkten hervorschauen.