senoidal - Generar señal sinusoidal en C sin utilizar la función estándar
graficar espectro de frecuencia en matlab (11)
Quiero generar una señal sinusoidal en C sin usar la función estándar sin () para desencadenar cambios sinusoidales en el brillo de un LED. Mi idea básica era utilizar una tabla de búsqueda con 40 puntos e interpolación.
Aquí está mi primer acercamiento:
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;
}
}
Desafortunadamente, esta función a veces devuelve valores mucho mayores que 1. Además, la interpolación no parece ser buena (utilicé esto para crear cambios de brillo en forma de seno de un LED, pero estos son muy poco sofisticados).
¿Alguien tiene una mejor idea para implementar un generador de seno en C?
¿Ha considerado modelar la porción de la curva sinusoidal de [0..PI] como una parábola? Si el brillo del LED solo debe ser observado por un ojo humano, las formas de las curvas deberían ser lo suficientemente similares para que se detecte una pequeña diferencia.
Solo necesitarías averiguar la ecuación apropiada para describirla.
Mmm, ...
Vértice en (PI / 2, 1)
Intersecciones del eje X en (0, 0) y (PI, 0)
f(x) = 1 - K * (x - PI/2) * (x - PI/2)
Donde K estaría ...
K = 4 / (PI * PI)
A menos que su aplicación requiera una precisión real, no se suicide ideando un algoritmo para una onda sinusoidal o coseno de 40 puntos. Además, los valores en su tabla deben coincidir con el rango de entrada pwm de su LED.
Dicho esto, le eché un vistazo a tu código, se imprimió y pensé que no estabas interpolando entre puntos. Con una pequeña modificación, lo arreglé, y el error entre la función de signo de excel y el suyo está desactivado en un máximo de aproximadamente 0.0032 o aproximadamente. El cambio es bastante fácil de implementar y se ha probado utilizando tcc, mi opción personal para las pruebas del algoritmo C
En primer lugar, agregué un punto más a tu matriz sinusoidal. El último punto se establece en el mismo valor que el primer elemento de la matriz senoidal. Esto corrige las matemáticas en su función seno, en particular cuando establece x1 en (int) fase% 40, y x2 en x1 + 1. No es necesario agregar el punto extra, ya que podría establecer x2 en (x1 + 1)% 40, pero elegí el primer enfoque. Solo estoy señalando diferentes maneras en que podrías lograr esto. También agregué el cálculo de un resto (Básicamente fase - (fase) fase). Estoy usando el resto para la interpolación. También agregué un valor temporal de seno y una variable delta.
const int sine_table[41] =
{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, 0};
int i = 0;
int x1 = 0;
int x2 = 0;
float y = 0;
float sin1(float phase)
{
int tsv,delta;
float rem;
rem = phase - (int)phase;
x1 = (int) phase % 40;
x2 = (x1 + 1);
tsv=sine_table[x1];
delta=sine_table[x2]-tsv;
y = tsv + (int)(rem*delta);
return y;
}
int main()
{
int i;
for(i=0;i<420;i++)
{
printf("%.2f, %f/n",0.1*i,sin1(0.1*i)/32768);
}
return 0;
}
Los resultados se ven bastante bien. La comparación de la aproximación lineal con la función senoidal de punto flotante del sistema me dio la gráfica de error que se muestra a continuación.
El principal problema de OP es generar el índice para la búsqueda de tablas.
El código de OP intenta acceder a la matriz sine_table[40] conduce a un comportamiento indefinido . Arregla eso al menos.
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
Cambio sugerido
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])*...
Existen enfoques mucho mejores, pero para centrarse en el método de OP, ver más abajo.
#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;
}
}
Salida
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
Un código mejor no pasaría los valores i, x1, x2, y como variables globales, sino como parámetros de función o variables de función. Quizás ese es un artefacto de la depuración de OP.
¿Alguien tiene una mejor idea para implementar un generador de seno en C?
Esto es bastante amplio. ¿Mejor que en velocidad, precisión, espacio de código, portabilidad o mantenimiento? sine() son fáciles de hacer. Los de alta calidad requieren más esfuerzo.
Aunque borroso, el uso de una pequeña tabla de consulta por parte de OP es un buen comienzo, aunque veo que se puede hacer sin ninguna matemática de punto flotante. Recomiendo a OP que construya una solución probada y de trabajo y la publique en Code Review para obtener ideas de mejora.
El truco clásico para dibujar un círculo (y por lo tanto generar una onda sinusoidal también) es Hakmem # 149 de Marvin Minsky . P.ej,:
#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);
}
}
Será un poco excéntrico, no un círculo perfecto, y puede obtener algunos valores ligeramente superiores a 1, pero puede ajustarlos dividiendo por el máximo o redondeando. Además, se puede usar la aritmética de enteros, y puedes eliminar la multiplicación / división usando una potencia negativa de dos para e , así que en su lugar se puede usar el desplazamiento.
La generación de una función sinusoidal precisa requiere una cantidad de recursos (ciclos de CPU y memoria) que no están justificados en esta aplicación. Su objetivo de generar una curva sinusoidal "suave" no está considerando los requisitos de la aplicación.
Mientras que cuando dibuja la curva puede observar imperfecciones, cuando aplica esa curva a una unidad LED PWM, el ojo humano no percibirá esas imperfecciones en absoluto.
Tampoco es probable que el ojo humano perciba la diferencia de brillo entre valores adyacentes incluso en una curva de 40 pasos, por lo que no es necesaria la interpolación.
En general, será más eficiente si genera una función seno que genere los valores de la unidad PWM apropiados directamente sin punto flotante. De hecho, en lugar de una función sinusoidal, un coseno elevado escalado sería más apropiado, de modo que una entrada de cero da como resultado una salida de cero, y una entrada de la mitad del número de valores en el ciclo da como resultado el valor máximo para su unidad PWM .
La siguiente función genera una curva de coseno aumentada para un PWM FSD de 8 bits a partir de una búsqueda de 16 valores (y 16 bytes) que genera un ciclo de 59 pasos. Por lo tanto, es eficiente en memoria y rendimiento en comparación con su implementación de punto flotante de 40 pasos.
#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 ;
}
Para una entrada de 0 <= theta < PWM_THETA_MAX la curva se ve así:
Que es lo que sugiero bastante suave suficiente para la iluminación.
En la práctica puedes usarlo así:
for(;;)
{
for( unsigned i = 0; i < PWM_THETA_MAX; i++ )
{
LedPwmDrive( RaisedCosine8bit( i ) ) ;
Delay( LED_UPDATE_DLEAY ) ;
}
}
Si su rango de PWM no es de 0 a 255, simplemente escale la salida de la función; La resolución de 8 bits es más que suficiente para la tarea.
Le ayudaría si explica por qué no quiere la función incorporada, pero como han dicho otros, la serie de Taylor es una forma de estimar el valor. Sin embargo, las otras respuestas parecen estar utilizando realmente la serie Maclaurin, no Taylor. Debe tener una tabla de búsqueda de seno y coseno. Luego encuentre x 0 , el valor x más cercano en su tabla de búsqueda al x que desea, y encuentre d = xx 0 . Entonces
sin (x) = sin (x 0 ) + cos (x 0 ) * d-sin (x 0 ) * d 2 /2-cos (x 0 ) * d 3/6 + ...
Si su tabla de búsqueda es tal que d <.01, entonces obtendrá más de dos dígitos de precisión por término.
Otro método es usar el hecho de que si x = x 0 + d, entonces
sin (x) = sin (x 0 ) * cos (d) + cos (x 0 ) * sin (d)
Puede usar una tabla de búsqueda para obtener sin (x 0 ) y cos (x 0 ), y luego usar la serie Maclaurin para obtener cos (d) y sin (d).
Me gustaría ir con la aproximación de Bhaskara I de una función sinusoidal. Usando grados, de 0 a 180, puede aproximar el valor de esta manera
float Sine0to180(float phase)
{
return (4.0f * phase) * (180.0f - phase) / (40500.0f - phase * (180.0f - phase));
}
Si quieres tener en cuenta cualquier ángulo, agregarías
float sine(float phase)
{
float FactorFor180to360 = -1 * (((int) phase / 180) % 2 );
float AbsoluteSineValue = Sine0to180(phase - (float)(180 * (int)(phase/180)));
return AbsoluteSineValue * FactorFor180to360;
}
Si quieres hacerlo en radianes, añadirías
float SineRads(float phase)
{
return Sine(phase * 180.0f / 3.1416);
}
Aquí hay una gráfica que muestra los puntos calculados con esta aproximación y también los puntos calculados con la función seno. Apenas se pueden ver los puntos de aproximación que sobresalen de los puntos sinusoidales reales.
Para un LED, probablemente podría hacerlo con 16 o más pasos sin siquiera interpolar. Dicho esto, puedo ver al menos dos cosas extrañas en su función sin1() :
1) Tiene 40 puntos de datos en sine_table , pero está tomando el índice x1 módulo 41 de la entrada. Esa no parece la forma correcta de manejar la periodicidad, y permite que x1 apunte más allá del último índice de la matriz.
2) Luego estás agregando +1, por lo que x2 puede estar incluso más allá de los límites de la matriz.
3) Estás usando i en la función, pero solo está configurado en el programa principal. No puedo decir qué se supone que debe hacer, pero usar un global como ese en una simple función de cálculo parece estar sucio como mínimo. Tal vez se supone que debe proporcionar la parte fraccionaria para la interpolación, pero no deberías usar la phase para eso.
Aquí hay un simple interpolador, que parece funcionar. Ajustar al gusto.
#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]);
}
Algunas salidas de muestra arbitrarias:
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
(Dejaré al lector interesado que reemplace todas las ocurrencias de 3 con una constante simbólica correctamente definida, para generalizar aún más la función y para implementar el cálculo de la fase negativa).
Podrías usar los primeros términos de la expansión de la serie de Taylor del sin . Puede usar la menor cantidad de términos que sea necesario para alcanzar el nivel de precisión deseado. Unos pocos términos más que el siguiente ejemplo deberían comenzar a toparse con los límites de una flotación de 32 bits.
Ejemplo:
#include <stdio.h>
// Please use the built-in floor function if you can.
float my_floor(float f) {
return (float) (int) f;
}
// Please use the built-in fmod function if you can.
float my_fmod(float f, float n) {
return f - n * my_floor(f / n);
}
// t should be in given in radians.
float sin_t(float t) {
const float PI = 3.14159265359f;
// First we clamp t to the interval [0, 2*pi)
// because this approximation loses precision for
// values of t not close to 0. We do this by
// taking fmod(t, 2*pi) because sin is a periodic
// function with period 2*pi.
t = my_fmod(t, 2.0f * PI);
// Next we clamp to [-pi, pi] to get our t as
// close to 0 as possible. We "reflect" any values
// greater than pi by subtracting them from pi. This
// works because sin is an odd function and so
// sin(-t) = -sin(t), and the particular shape of sin
// combined with the choice of pi as the endpoint
// takes care of the negative.
if (t >= PI) {
t = PI - t;
}
// You can precompute these if you want, but
// the compiler will probably optimize them out.
// These are the reciprocals of odd factorials.
// (1/n! for odd n)
const float c0 = 1.0f;
const float c1 = c0 / (2.0f * 3.0f);
const float c2 = c1 / (4.0f * 5.0f);
const float c3 = c2 / (6.0f * 7.0f);
const float c4 = c3 / (8.0f * 9.0f);
const float c5 = c4 / (10.0f * 11.0f);
const float c6 = c5 / (12.0f * 13.0f);
const float c7 = c6 / (14.0f * 15.0f);
const float c8 = c7 / (16.0f * 17.0f);
// Increasing odd powers of t.
const float t3 = t * t * t;
const float t5 = t3 * t * t;
const float t7 = t5 * t * t;
const float t9 = t7 * t * t;
const float t11 = t9 * t * t;
const float t13 = t9 * t * t;
const float t15 = t9 * t * t;
const float t17 = t9 * t * t;
return c0 * t - c1 * t3 + c2 * t5 - c3 * t7 + c4 * t9 - c5 * t11 + c6 * t13 - c7 * t15 + c8 * t17;
}
// Test the output
int main() {
const float PI = 3.14159265359f;
float t;
for (t = 0.0f; t < 12.0f * PI; t += (PI * 0.25f)) {
printf("sin(%f) = %f/n", t, sin_t(t));
}
return 0;
}
Ejemplo de salida:
sin(0.000000) = 0.000000
sin(0.785398) = 0.707107
sin(1.570796) = 1.000000
sin(2.356194) = 0.707098
sin(3.141593) = 0.000000
sin(3.926991) = -0.707107
sin(4.712389) = -1.000000
sin(5.497787) = -0.707098
sin(6.283185) = 0.000398
...
sin(31.415936) = 0.000008
sin(32.201332) = 0.707111
sin(32.986729) = 1.000000
sin(33.772125) = 0.707096
sin(34.557522) = -0.000001
sin(35.342918) = -0.707106
sin(36.128315) = -1.000000
sin(36.913712) = -0.707100
sin(37.699108) = 0.000393
Como puede ver, todavía hay espacio para mejorar la precisión. No soy un genio con aritmética de punto flotante, así que probablemente parte de esto tenga que ver con las implementaciones floor / fmod o el orden específico en el que se realizan las operaciones matemáticas.
Ya que intenta generar una señal, creo que usar una ecuación diferencial no debería ser una mala idea. da algo asi
#include <stdlib.h>
#include <stdio.h>
#define DT (0.01f) //1/s
#define W0 (3) //rad/s
int main(void) {
float a = 0.0f;
float b = DT * W0;
float tmp;
for (int i = 0; i < 400; i++) {
tmp = (1 / (1 + (DT * DT * W0 * W0))) * (2 * a - b);
b = a;
a = tmp;
printf("%f/n", tmp);
}
}
Aún configurando la amplitud y frecuencia de la señal es un dolor en el cuello: /
... una mejor idea para implementar un generador de seno en C?
Edición: sugiera la primera lectura de este artículo para obtener una apreciación de lo que está pidiendo OP.
El algoritmo CORDIC ( Ordenador digital de rotación coordinada ) es muy adecuado para su uso en implementaciones pequeñas uP, y FPGA que tienen capacidades de cálculo matemático limitadas, ya que calcula el seno y el coseno de un valor utilizando solo aritmética básica (suma, resta y cambios) . Aquí se proporciona más información sobre CORDIC y cómo usarla para producir seno / coseno de un ángulo.
También hay varios sitios que proporcionan ejemplos de implementación de algoritmos. CORDIC simple es uno que incluye explicaciones detalladas sobre cómo generar una tabla que luego se puede precompilar para usar en su dispositivo de destino, así como un código para probar la salida de la siguiente función (que usa matemática de punto fijo):
(Ver documentación de seguimiento, y otras funciones en enlace)
#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;
}
Editar:
Encontré la documentación para usar los ejemplos en el sitio de CORDIC simple muy fácil de seguir. Sin embargo, una pequeña cosa con la que me topé fue al compilar el archivo cordic-test.c el error: se usó el identificador no declarado ''M_PI'' . Parece que al ejecutar el archivo gentable.c compilado (que genera el archivo cordic-test.c ) la línea:
#define M_PI 3.1415926535897932384626
aunque se incluyó en sus propias declaraciones, no se incluyó en las declaraciones de printf utilizadas para producir el archivo cordic-test.c . Una vez que esto fue remediado, todo funcionó como se anuncia.
Como se documentó, el rango de datos producidos genera 1/4 de un ciclo sinusoidal completo (-π / 2 - π / 2). La siguiente ilustración contiene una representación de los datos reales producidos entre los puntos azules claros. El resto de la señal sinusoidal se fabrica a través de la duplicación y la transposición de la sección de datos originales.
