En lugar de calcular el seno como una función del tiempo, mantenga un par seno/coseno y avance a través de la multiplicación de números complejos. Esto no requiere funciones trigonométricas ni tablas de búsqueda; solo cuatro multiplicaciones y una renormalización ocasional:

 static const double a = 2 * M_PI * 280 * 30e-6; static const double dx = cos(a); static const double dy = sin(a); double x = 1, y = 0; // complex x + iy int counter = 0; void control_loop() { double xx = dx*x - dy*y; double yy = dx*y + dy*x; x = xx, y = yy; // renormalize once in a while, based on // https://www.gamedev.net/forums/topic.asp?topic_id=278849 if((counter++ & 0xff) == 0) { double d = 1 - (x*x + y*y - 1)/2; x *= d, y *= d; } double sine = y; // this is your sine }

La frecuencia se puede ajustar, si es necesario, volviendo a calcular dx , dy .

Además, todas las operaciones aquí se pueden hacer, bastante fácilmente, en punto fijo.

Racionalidad

Como @ user3386109 señala a continuación (+1) , el 280 * 30e-6 = 21 / 2500 es un número racional, por lo que el seno debería dar la vuelta después de 2500 muestras exactamente . Podemos combinar este método con el de ellos reiniciando nuestro generador ( x=1,y=0 ) cada 2500 iteraciones (o 5000, o 10000, etc...). Esto eliminaría la necesidad de volver a normalizar, así como eliminar cualquier imprecisión de fase a largo plazo.

(Técnicamente, cualquier número de punto flotante es un racional diádico. Sin embargo, 280 * 30e-6 no tiene una representación exacta en binario. Sin embargo, al reiniciar el generador como se sugiere, obtendremos un seno exactamente periódico como se esperaba).

Explicación

Algunos pidieron una explicación en los comentarios de por qué esto funciona. La explicación más simple es usar las identidades trigonométricas de suma de ángulos :

 xx = cos((n+1)*a) = cos(n*a)*cos(a) - sin(n*a)*sin(a) = x*dx - y*dy yy = sin((n+1)*a) = sin(n*a)*cos(a) + cos(n*a)*sin(a) = y*dx + x*dy

y la corrección se sigue por inducción.

Esta es esencialmente la fórmula de De Moivre si consideramos esos pares seno/coseno como números complejos, de acuerdo con la fórmula de Euler .

Una forma más perspicaz podría ser mirarlo geométricamente. La multiplicación compleja por exp(ia) es equivalente a la rotación a radianes. Por lo tanto, al multiplicar repetidamente por dx + idy = exp(ia) , rotamos incrementalmente nuestro punto de inicio 1 + 0i a lo largo del círculo unitario. La coordenada y , según la fórmula de Euler nuevamente, es el seno de la fase actual.

Normalización

Mientras la fase continúa avanzando con cada iteración, la magnitud (también conocida como norma) de x + iy se aleja de 1 debido a errores de redondeo. Sin embargo, estamos interesados en generar un seno de amplitud 1 , por lo que necesitamos normalizar x + iy para compensar la deriva numérica. La forma directa es, por supuesto, dividirlo por su propia norma:

 double d = 1/sqrt(x*x + y*y); x *= d, y *= d;

Esto requiere un cálculo de una raíz cuadrada recíproca. Aunque normalizamos solo una vez cada X iteraciones, sería bueno evitarlo. Afortunadamente |x + iy| ya está cerca de 1 , por lo que solo necesitamos una ligera corrección para mantenerlo a raya. Expandiendo la expresión para d alrededor 1 (aproximación de Taylor de primer orden), obtenemos la fórmula que está en el código:

 d = 1 - (x*x + y*y - 1)/2

TAREAS: para comprender completamente la validez de esta aproximación, es necesario demostrar que compensa los errores de redondeo más rápido de lo que se acumulan y, por lo tanto, obtener un límite de la frecuencia con la que se debe aplicar.

Como se señaló en algunos de los comentarios, el valor del time crece continuamente con el tiempo. Esto plantea dos problemas:

1. Es probable que la función sin tenga que realizar un módulo internamente para obtener el valor interno en un rango admitido.
2. La resolución del tiempo empeorará cada vez más a medida que aumente el valor, debido a la adición de dígitos más altos.

Hacer los siguientes cambios debería mejorar el rendimiento:

 double time; void control_loop() { time += 30.0e-6; if((1.0/280.0) < time) { time -= 1.0/280.0; } double sine = sin(2 * M_PI * 280 * time); ... }

Tenga en cuenta que una vez realizado este cambio, ya no tendrá una variable de tiempo.

La función se puede reescribir como

 double n; void control_loop() { n += 1; double sine = sin(2 * M_PI * 280 * 30e-6 * n); ... }

Eso hace exactamente lo mismo que el código en la pregunta, con exactamente los mismos problemas. Pero ahora se puede simplificar:

 280 * 30e-6 = 280 * 30 / 1000000 = 21 / 2500 = 8.4e-3

Lo que significa que cuando n llega a 2500, ha generado exactamente 21 ciclos de la onda sinusoidal. Lo que significa que puede establecer n de nuevo en 0. El código resultante es:

 int n; void control_loop() { n += 1; if (n == 2500) n = 0; double sine = sin(2 * M_PI * 8.4e-3 * n); ... }

Siempre que su código pueda ejecutarse durante 21 ciclos sin problemas, se ejecutará para siempre sin problemas.

Utilice una tabla de consulta. Su comentario en la discusión con Eugene Sh.:

Una pequeña desviación de la frecuencia sinusoidal (como 280,1 Hz) estaría bien.

En ese caso, con un intervalo de control de 30 µs, si tienes una tabla de 119 muestras que repites una y otra vez, obtendrás una onda sinusoidal de 280,112 Hz. Dado que tiene un DAC de 12 bits, solo necesita 119 * 2 = 238 bytes para almacenar esto si lo envía directamente al DAC. Si lo usa como entrada para cálculos adicionales como los que menciona en los comentarios, puede almacenarlo como float o double según lo desee. En una MCU con RAM estática integrada, solo se necesitan unos pocos ciclos como máximo para cargar desde la memoria.

Si tiene algunos kilobytes de memoria disponibles, puede eliminar este problema por completo con una tabla de búsqueda.

Con un período de muestreo de 30 µs, 2500 muestras tendrán una duración total de 75 ms. Esto es exactamente igual a la duración de 21 ciclos a 280 Hz.

No he probado ni compilado el siguiente código, pero al menos debería demostrar el enfoque:

 double sin2500() { static double *table = NULL; static int n = 2499; if (!table) { table = malloc(2500 * sizeof(double)); for (int i=0; i<2500; i++) table[i] = sin(2 * M_PI * 280 * i * 30e-06); } n = (n+1) % 2500; return table[n]; }

Cómo generar un hermoso seno.

DAC es de 12 bits, por lo que solo tiene 4096 niveles. No tiene sentido enviar más de 4096 muestras por período. En la vida real, necesitará muchas menos muestras para generar una forma de onda de buena calidad.

1. Cree un archivo C con la tabla de búsqueda (usando su PC). Redirija la salida al archivo ( https://helpdeskgeek.com/how-to/redirect-output-from-command-line-to-text-file/ ).

 #define STEP ((2*M_PI) / 4096.0) int main(void) { double alpha = 0; printf("#include <stdint.h>\nconst uint16_t sine[4096] = {\n"); for(int x = 0; x < 4096 / 16; x++) { for(int y = 0; y < 16; y++) { printf("%d, ", (int)(4095 * (sin(alpha) + 1.0) / 2.0)); alpha += STEP; } printf("\n"); } printf("};\n"); }

https://godbolt.org/z/e899d98oW

2. Configure el temporizador para activar el desbordamiento 4096*280=1146880 veces por segundo. Configure el temporizador para generar el evento de activación DAC. Para el reloj temporizador de 180 MHz, no será preciso y la frecuencia será de 279,906449045 Hz. Si necesita una mayor precisión, cambie el número de muestras para que coincida con la frecuencia de su temporizador o cambie la frecuencia del reloj del temporizador (los temporizadores H7 pueden funcionar hasta 480 MHz)

3. Configure DAC para usar DMA y transfiera el valor de la tabla de búsqueda creada en el paso 1 a DAC en el evento desencadenante.

4. Disfrute de una hermosa onda sinusoidal usando su osciloscopio. Tenga en cuenta que el núcleo de su microcontrolador no se cargará en absoluto. Lo tendrás para otras tareas. Si desea cambiar el período, simplemente reconfigure el temporizador. Puedes hacerlo tantas veces por segundo como desees. Para reconfigurar el temporizador, utilice el modo de ráfaga DMA del temporizador, que volverá a cargar los registros PSC y ARR en el evento de actualización automáticamente sin perturbar la forma de onda generada.

Sé que es programación STM32 avanzada y requerirá programación a nivel de registro. Lo uso para generar formas de onda complejas en nuestros dispositivos.

Es la forma correcta de hacerlo. Sin bucles de control, sin cálculos, sin carga central.

Creo que sería posible usar un módulo porque sin() es periódico.

Entonces no tienes que preocuparte por los problemas.

 double time = 0; long unsigned int timesteps = 0; double sine; void controll_loop() { timesteps++; time += 30e-6; if( time > 1 ) { time -= 1; } sine = sin( 2 * M_PI * 280 * time ); ... }

¿Qué tal una variante del concepto basado en módulos de otros?

 int t = 0; int divisor = 1000000; void control_loop() { t += 30 * 280; if (t > divisor) t -= divisor; double sine = sin(2 * M_PI * t / (double)divisor)); ... }

Calcula el módulo en un número entero y luego no provoca errores de redondeo.

Estoy bastante sorprendido por las respuestas existentes. El primer problema que detecta se resuelve fácilmente y el siguiente problema desaparece mágicamente cuando resuelve el primero.

Necesita una comprensión básica de las matemáticas para ver cómo funciona. Recuerda, sin(x+2pi) es simplemente sin(x) , matemáticamente. El gran aumento en el tiempo que ve ocurre cuando su implementación sin(float) cambia a otro algoritmo, y realmente quiere evitar eso.

Recuerda que float tiene solo 6 dígitos significativos. 100000.0f*M_PI+x usa esos 6 dígitos para 100000.0f*M_PI , por lo que no queda nada para x .

Entonces, la solución más fácil es realizar un seguimiento de x usted mismo. En t=0 , inicializa x a 0.0f . Cada 30 us, incrementas x+= M_PI * 280 * 30e-06; . ¡El tiempo no aparece en esta fórmula! Finalmente, si x>2*M_PI , decrementas x-=2*M_PI; (Ya que sin(x)==sin(x-2*pi)

Ahora tiene una x que se mantiene muy bien en el rango de 0 a 6.2834 , donde sin es rápido y los 6 dígitos de precisión son todos útiles.

Existe un enfoque alternativo para calcular una serie de valores de seno (y coseno) para ángulos que aumentan en una cantidad muy pequeña. Esencialmente se reduce a calcular las coordenadas X e Y de un círculo, y luego dividir el valor de Y por alguna constante para producir el seno, y dividir el valor de X por la misma constante para producir el coseno.

Si está satisfecho con generar una "elipse muy redonda", puede usar un siguiente truco, que se atribuye a Marvin Minsky en la década de 1960. Es mucho más rápido que calcular senos y cosenos, aunque introduce un error muy pequeño en la serie. Aquí hay un extracto del Documento Hakmem, Artículo 149. Se describe el algoritmo del círculo de Minsky.

PUNTO 149 (Minsky): ALGORITMO DEL CÍRCULO He aquí una forma elegante de dibujar casi círculos en una pantalla de trazado de puntos:

 NEW X = OLD X - epsilon * OLD Y NEW Y = OLD Y + epsilon * NEW(!) X

Esto hace una elipse muy redonda centrada en el origen con su tamaño determinado por el punto inicial. epsilon determina la velocidad angular del punto de circulación y afecta ligeramente la excentricidad. Si épsilon es una potencia de 2, entonces ni siquiera necesitamos la multiplicación, ¡y mucho menos las raíces cuadradas, los senos y los cosenos! El "círculo" será perfectamente estable porque los puntos pronto se volverán periódicos.

¡El algoritmo circular se inventó por error cuando intenté guardar un registro en un truco de visualización! Ben Gurley tuvo un increíble truco de visualización utilizando solo unas seis o siete instrucciones, y fue una gran maravilla. Pero estaba básicamente orientado a la línea. Se me ocurrió que sería emocionante tener curvas, y estaba tratando de obtener un truco de visualización de curvas con instrucciones mínimas.

Aquí hay un enlace al hakmem: http://inwap.com/pdp10/hbaker/hakmem/hacks.html

Me gustaría abordar los problemas de programación integrados en su código directamente: la respuesta de @ 0___________ es la forma correcta de hacer esto en un microcontrolador y no volveré a pisar el mismo terreno.

• Las variables que representan el tiempo nunca deben ser de coma flotante. Si su incremento no es una potencia de dos, los errores siempre se acumularán. Incluso si lo es, eventualmente su incremento será más pequeño que el incremento más pequeño y el temporizador se detendrá. Siempre use números enteros para el tiempo. Puede elegir un tamaño de entero lo suficientemente grande como para ignorar el traspaso: un entero de 32 bits sin signo que represente milisegundos tardará 50 días en traspasarse, mientras que un entero de 64 bits sin signo tardará más de 500 millones de años.
• Generar cualquier señal periódica en la que no le importe la fase de la señal no requiere una variable de tiempo. En su lugar, puede mantener un contador interno que se restablece a 0 al final de un período. (Cuando usa DMA con una tabla de búsqueda, eso es exactamente lo que está haciendo: el contador es el puntero de próxima lectura del controlador DMA).
• Siempre que use una función trascendental como seno en un microcontrolador, su primer pensamiento debería ser "¿puedo usar una tabla de búsqueda para esto?" No tiene acceso al lujo de un sistema operativo moderno que distribuye su carga de manera óptima en un procesador multinúcleo de más de 4 GHz. A menudo se trata de un solo hilo que se detendrá esperando que su microcontrolador de 200 MHz saque la FPU del modo de espera y realice el algoritmo de aproximación. Hay un costo significativo para las funciones trascendentales. Las LUT también tienen un costo, pero si está presionando la función constantemente, es muy probable que le gusten mucho más las compensaciones de la LUT.

Hilo fascinante. El algoritmo de Minsky mencionado en la respuesta de Walter Mitty me recordó un método para dibujar círculos que se publicó en Electronics & Wireless World y que conservé. (Crédito: https://www.electronicsworld.co.uk/magazines/ ). Lo adjunto aquí por interés.

Sin embargo, para mis propios proyectos similares (para la síntesis de audio) utilizo una tabla de búsqueda, con suficientes puntos para que la interpolación lineal sea lo suficientemente precisa (¡haga los cálculos!)