Interpolación Spline Cúbica - Resolver X desde Y

Question

Soy programador, no matemático, pero tengo un problema del mundo real que estoy tratando de resolver que está fuera de mi alcance, y mis habilidades en Google hasta ahora me han fallado.

Tengo una forma de onda analógica que ha sido muestreada, parte de la cual contiene una onda sinusoidal, como se muestra a continuación: La frecuencia de muestreo de los datos y la frecuencia de la onda sinusoidal son constantes, pero no están relacionadas. Necesito tomar esta forma de onda muestreada, y "upsamplearla", es decir, generar más muestras entre los puntos de muestreo existentes. Para hacer esto, he estado utilizando la interpolación de Spline cúbico. He tomado una biblioteca de C ++ existente de la siguiente ubicación: http://kluge.in-chemnitz.de/opensource/spline/ . Las fórmulas se dan en la página. Paso los datos muestreados existentes como puntos de control, con "X" ajustado al número de muestra, e "Y" ajustado al valor de la muestra. Dado un valor de entrada de X, esta biblioteca me genera el valor "Y" correspondiente, por lo que pasar una X de 2,5 generaría un valor en la spline generada a medio camino entre las muestras 2 y 3. Todo esto está funcionando muy bien hasta ahora, y produce una forma de onda como la que se muestra arriba, pero ahora necesito hacer una cosa extra en la que estoy un poco atascado.

Lo que necesito hacer ahora es "sincronizar" con la onda sinusoidal en los datos de la muestra. Para ello, quiero detectar el punto en el que la onda sinusoidal cruza un umbral determinado, como se muestra a continuación: En otras palabras, dado un valor inicial de X, quiero calcular el siguiente valor de X en el que Y es igual a un valor determinado. En la actualidad, "fuerzo" una solución aproximada, incrementando en un tamaño de paso dado entre dos valores de X hasta que cruza el punto de umbral. Esto plantea dos problemas:

1. Precisión. La solución es sólo aproximada, y necesito un resultado muy preciso. Esto conduce a un tamaño de paso muy pequeño, sin dejar de tener un margen de error.
2. Rendimiento. Tengo literalmente miles de millones de estas operaciones para realizar, y quiero algo que pueda evaluar en un período de tiempo razonable. Ahora mismo tengo que sacrificar la precisión para mejorar el rendimiento.

En lugar de utilizar una búsqueda de fuerza bruta imperfecta como hago actualmente, me preguntaba si hay una solución matemática para este problema. Básicamente quiero una función que, dada una posición inicial X y un valor de muestra Y objetivo, devuelva la siguiente X en la que Y sea igual al valor Y objetivo. Creo que esto debería ser posible de lograr, pero me faltan los conocimientos matemáticos para derivar la fórmula yo mismo. ¿Puede alguien decirme si tal función es teóricamente posible de escribir, y si es así, cuál sería una posible implementación?

Aquí hay algunas restricciones potencialmente útiles en mis datos de entrada, que cualquier solución aceptable puede asumir:

• Las coordenadas X de los puntos de control son valores linealmente crecientes (es decir, todos están espaciados uniformemente y nunca se "duplican")
• Conozco los límites superior e inferior de las coordenadas X e Y de los puntos de control al generar la spline

Por favor, dígame si puede ver una solución a este problema.

Answer 1

Básicamente, tienes un problema de búsqueda de raíces de polinomios.

Usted tiene un $y$ -valor, digamos $y = \bar{y}$ y quieres encontrar $x$ tal que $y(x) = \bar{y}$ . En otras palabras, se quiere encontrar una raíz de la función $f(x) = y(x) - \bar{y}$ . La función $f$ es un polinomio cúbico.

Algunas ideas:

(1) Busca una función de búsqueda de raíces de uso general. Busca en el libro "Numerical Recipes", por ejemplo. El método de Brent parece ser una opción popular. Puedes introducir tu función $f$ directamente, por lo que hay que escribir menos código que en las sugerencias siguientes.

(2) Busca algún código de búsqueda de raíces de polinomios. De nuevo, hay algunos en el libro "Recetas Numéricas". Esto puede requerir que encuentres los coeficientes polinómicos de los segmentos de tu spline cúbico, sin embargo. No es difícil, pero es trabajo.

(3) Existen fórmulas para encontrar las raíces de los polinomios cúbicos (ver Wikipedia), pero codificarlas correctamente es sorprendentemente difícil. Hay un buen código aquí y aquí . Esto probablemente le dará el mejor rendimiento con los cúbicos. De nuevo, requerirá que encuentres los coeficientes polinómicos de los segmentos de tu spline cúbica.

(4) Pasar a utilizar splines cuadráticos en lugar de cúbicos. Entonces sólo tendrás que encontrar raíces de polinomios cuadráticos, no cúbicos, por lo que el problema de búsqueda de raíces se vuelve trivial. Este será un cambio bastante grande, pero valdrá la pena si necesitas un alto rendimiento.

Answer 2

No tengo ni idea de cuál es la mejor manera de solucionar este problema, así que espero que alguien venga y le dé algo de magia, tras lo cual borraré esto con mucho gusto.

He pensado un poco en las ventajas que tiene la búsqueda por bisección, pero me he convencido, quizás erróneamente, de que no va a superar a la fórmula cúbica debido a algunos casos complicados en los que el umbral interseca dos veces el mismo trozo de spline. Como este es un caso bastante plausible, creo que la fórmula cúbica superará a la búsqueda por bisección. Si esto no es cierto, una búsqueda de bisección va a ser más fácil de escribir y probablemente más barato en promedio que la fórmula cúbica. Newton-Raphson tiene la mejor convergencia, pero no creo que sea muy buena en este caso, puesto que ya te falta una buena estimación. Sin embargo, una búsqueda inicial de bisección seguido de Newton-Raphson sería dinamita para la convergencia, pero creo que en el caso de un cúbico que va a realizar muy cerca del mismo número de operaciones, así que de nuevo estoy adivinando la fórmula cúbica ganará. La fórmula cúbica requiere un montón de multiplicaciones y algo de aritmética compleja, así que realmente no tengo claro cuál es la mejor, y desde un punto de vista numérico definitivamente no tengo ni idea. Estoy publicando esto principalmente porque nadie más ha respondido y tal vez una pequeña respuesta es mejor que nada.

Pero, no tengo 100% claro si hay que encontrar todos los puntos de intersección, o sólo un punto de intersección. Si sólo tiene que encontrar un punto de intersección para fines de sincronización, por favor, utilice la bisección y el método de Newton. Va a ganar.

Si quieres dar una oportunidad a la fórmula cúbica Si no tienes una biblioteca de números complejos, tendrás que conseguir los coeficientes de la spline modificando la biblioteca para que los vectores de los coeficientes sean públicos. Estarás tomando la raíz cúbica de un número negativo.

Si su spline tiene coeficientes $ax^3 + bx^2 + cx + D$ y su umbral está en $y = Y$ entonces quieres encontrar los ceros de $ax^3 + bx^2 + cx + d$ donde $d = D-Y$ .

Calcular el discriminante $\Delta = 18abcd - 4b^3d + b^2c^2 - 4ac^3 -27a^2d^2$ . Esto lo quieres para saber si tu cúbica tiene una o tres raíces reales. Supongo que es teóricamente posible que puedas conseguir algo de magia de punto flotante para tener una raíz doble ( $\Delta = 0$ ).

A continuación, utilice la forma $x_k = -\frac{1}{3a}(b + \zeta^kC + \frac{\Delta_0}{\zeta^k C})$ para $k = 0,~ 1,~ 2$ . Tenga en cuenta el comentario de wiki sobre la simplificación del cálculo de $C$ desde que calculaste el discriminante.

Así que tu fórmula te da tres raíces, y tu discriminante dice que tienes

1. Una raíz real. Busca la raíz en tu conjunto de tres que se ajuste a la ventana en la que estás mirando actualmente. Teniendo en cuenta las parcelas que tienes en tu post no creo que vaya a haber ningún caso límite aquí. Simplemente encuentra la parte real correcta y descarta la parte imaginaria ya que es un error de punto flotante.
2. Tres raíces reales. Es posible que tu umbral cruce la spline en dos puntos (como una tapa o una depresión). La buena noticia es que puedes descartar todas las partes imaginarias como error de punto flotante y sólo buscar las partes reales en tu rango.

3. Una raíz doble.

   a. Las dos raíces que están muy cerca en magnitud están en su rango. Esto significa que tu umbral está justo en la punta de un pico o una depresión. Creo que no hay nada malo en ignorar este caso.

   b. Las dos raíces que están muy cerca en magnitud, sólo una está en su rango. Probablemente también puedes ignorar este caso.

   c. La tercera raíz no cercana en magnitud está en su rango. Bien, lo tienes, descarta la parte imaginaria como error de punto flotante.

Por si acaso esto no fuera suficiente para ti, tengo que señalar que la biblioteca almacena sus coeficientes en la forma $a(x-x_i)^3 + b(x-x_i)^2 + c(x-x_i) + y_i$ . Así que para saber si la raíz está "al alcance" hay que traducirlo por la ventana en la que se encuentra actualmente.

Answer 3

Como alternativa a las fórmulas generales de Cardano, que pueden ser eficientes cuando se dispone de un coprocesador de punto flotante, puedes intentar lo siguiente:

• comprobar primero que la cúbica no tiene ningún extremo en el intervalo considerado, dejemos $[p,q]$ para ello, localiza las raíces de $$3ax^2+2bx+c=0;$$

  • para ello, comprueba si la raíz de $6ax+2b=0$ se encuentra en el intervalo, comparando los signos de $3ap+b$ y $3aq+b$ .

  • si es verdadero, dividir en los subintervalos $[p,-\dfrac b{3a}]$ y $[-\dfrac b{3a},q]$ ;

  • evaluar los signos de $3ax^2+2bx+c$ en los puntos finales de los intervalos correspondientes (uno o dos). Cuando encuentres una raíz de la cuadrática, evalúala y divide el intervalo correspondiente.

• después de esta discusión preliminar, conocerás de cero a tres subintervalos con exactamente una raíz.

• finalmente, resolver mediante unas cuantas iteraciones del método de la cuerda o de Newton.

---

Otro enfoque interesante (que puede implementarse utilizando sólo aritmética de números enteros) es expresar el polinomio cúbico como un Bezier cúbico (hay fórmulas sencillas para traducir de las formas canónicas a las de Bezier). Entonces se puede utilizar

• la propiedad hull para comprobar si la curva cruza una ordenada dada (los puntos de control deben estar en ambos lados),

• el algoritmo de subdivisión (conocido como De Casteljau) para encontrar los puntos de control de las mitades izquierda y derecha de las curvas,

• recursión a las mitades en las que se detectan cruces.

Se puede organizar el cálculo para realizar sólo sumas y comparaciones, y el proceso de búsqueda de raíces es de tipo dicotómico. Se puede terminar con una interpolación lineal.

---

Como dijo @bubba, la precisión extrema es exagerada ya que su interpolante es en sí misma una aproximación.