Matematica Interpolacion sen(x) / x o similares (osciloscopio)

No hay un tema matematicas en el foro asi que lo puse aca pero si consideran que va en otro lado por favor moverlo.

Estoy terminando el soft de mi osciloscopio digital y me encuentro con un problema donde no tengo la suficiente matematica en mi cabeza para resolverlo, por lo tanto acudo aqui por si alguien puede darme una mano con esto.

Resulta que tengo un tiempo minimo entre muestras de 25ns y quiero, solo matematicamente, llevar el osciloscopio a valores inferiores de escala de tiempo. Aclaro esto para evitar caer en soluciones electronicas como el sampleo de señales periodicas con retraso de trigger (oversampling creo que es el termino correcto).

Entonces, si tengo los puntos Y de una señal senoidal para X = 1, X = 2, X = n... necesito fabricar con una funcion matematica los puntos X = 1.5, X = 2.5, etc.etc (haciendo el osciloscopio el doble de rapido matematicamente, luego el triple y asi sucesivamente hasta que sea razonable y no tenga una distorsion inaceptable de la señal real)

Estuve viendo los siguientes articulos de wikipedia:
http://en.wikipedia.org/wiki/Trigonometric_interpolation
http://en.wikipedia.org/wiki/Trigonometric_polynomial
http://en.wikipedia.org/wiki/Discrete_Fourier_transform

y otros, pero lamentablemente mi nivel matematico es demasiado bajo para entender estos articulos.

Encontre de pura suerte esto: http://www.arachnoid.com/polysolve/ que es una pagina con un codigo javascript para generar el polinomio de unos puntos de sample dados y que permite calcular cualquier punto necesario Y en el tiempo X. Aparentemente funciona bien segun unas pruebas que hice con unas muestras pero creo que esta resolucion es mas bien lineal, no senoidal.

Tambien vi que osciloscopios digitales comerciales tienen una opcion (por defecto activada) de interpolacion sen(x) / x . Busque acerca de esto pero no pude o no supe encontrar la informacion de como se aplica (no veo variable Y en esa ecuacion)

¿Alguien sabria orientarme como puedo hacer esta interpolacion? ¿Como podria aplicarla? Soy de profesion programador asi que eso no es problema, solo me falta la matematica para hacer la funcion!!

Desde ya, muchas gracias por leer esto
 
No entiendo lo que querés hacer :confused: :confused:
Si vos muestreás cada 25ns.. pues eso es todo lo que tenés, y no podés muestrear señales más rápìdas que 20 Mhz por que Shanon se vá al diablo. Nunca vas a saber que hay entre esos 25ns.. y si lo inventás, es solo eso... un invento. Las interpolaciones son para reconstrucción de la señal muestreada, pero no para generar muestras "intermedias" que desconocés cuanto valen.... ni hablar de hacer "más rapido" el osciloscopio, por que SI o SI debés tener la señal acotada en frecuencia a la mitad de la máxima fcia de muestreo. Si querés "reconstruir" los puntos intermedios, la única reconstrucción válida es la sen(x)/x, pero eso es solo la conversión digital --> analógico aunque no has aumentado ninguna velocidad ni nada por el estilo. Claro que luego podés moverte entre los puntos reconstruidos "como si fueran" los de la señal original pero no podés muestrear señales de mayor frecuencia.
 
No hay un tema matematicas en el foro asi que lo puse aca pero si consideran que va en otro lado por favor moverlo.

Estoy terminando el soft de mi osciloscopio digital y me encuentro con un problema donde no tengo la suficiente matematica en mi cabeza para resolverlo, por lo tanto acudo aqui por si alguien puede darme una mano con esto.

Resulta que tengo un tiempo minimo entre muestras de 25ns y quiero, solo matematicamente, llevar el osciloscopio a valores inferiores de escala de tiempo. Aclaro esto para evitar caer en soluciones electronicas como el sampleo de señales periodicas con retraso de trigger (oversampling creo que es el termino correcto).

Entonces, si tengo los puntos Y de una señal senoidal para X = 1, X = 2, X = n... necesito fabricar con una funcion matematica los puntos X = 1.5, X = 2.5, etc.etc (haciendo el osciloscopio el doble de rapido matematicamente, luego el triple y asi sucesivamente hasta que sea razonable y no tenga una distorsion inaceptable de la señal real)

Estuve viendo los siguientes articulos de wikipedia:
http://en.wikipedia.org/wiki/Trigonometric_interpolation
http://en.wikipedia.org/wiki/Trigonometric_polynomial
http://en.wikipedia.org/wiki/Discrete_Fourier_transform

y otros, pero lamentablemente mi nivel matematico es demasiado bajo para entender estos articulos.

Encontre de pura suerte esto: http://www.arachnoid.com/polysolve/ que es una pagina con un codigo javascript para generar el polinomio de unos puntos de sample dados y que permite calcular cualquier punto necesario Y en el tiempo X. Aparentemente funciona bien segun unas pruebas que hice con unas muestras pero creo que esta resolucion es mas bien lineal, no senoidal.

Tambien vi que osciloscopios digitales comerciales tienen una opcion (por defecto activada) de interpolacion sen(x) / x . Busque acerca de esto pero no pude o no supe encontrar la informacion de como se aplica (no veo variable Y en esa ecuacion)

¿Alguien sabria orientarme como puedo hacer esta interpolacion? ¿Como podria aplicarla? Soy de profesion programador asi que eso no es problema, solo me falta la matematica para hacer la funcion!!

Desde ya, muchas gracias por leer esto

sen(x)/x es la función sinc, que es la que según el teorema del muestreo permite reconstruir una senoidal perfecta a partir de sólo 2 muestras. En la práctica es suficiente upsamplear (introducir muestras 0 alternativamente entre cada sample, y después pasar un filtro FIR conformador), de esta forma simulas la función SINC. Por cierto, la FFT de SINC es un filtro pasabaja.

Ahora bien, igual de sencillo es hacer una interpolación por SPLINES. Tienes un vector de P puntos, pongamos que en el punto N, puedes aproximar la curva, con un polinomio de tercer grado, cuyos coeficientes se resuelven por un sencillo sistema de ecuaciones. Nos interesa que la curva sea suave, por lo que aproximas la derivada de ese polinomio en el punto N promediando las pendientes en los puntos N-1 y N+1 por la pendiente ((V[N+1]-V[N-1])/2). Ahora ya sabes que Ax^3+Bx^2+CX+D=0, necesitas 4 valores para resolver las incónigas A B C y D, es decir, necesitas 4 variables. Sabes que la onda pasará por el los tres puntos, ya tienes 3 condiciones, y la cuarta es la derivada de la gráfica en el punto medio. Haces un sistema de ecuaciones y sacas los cuatro coeficientes del polinomio.
 
Desde ya, les agradezco su tiempo en responder. Voy por partes:

Zoidberg: Si, no puedo tener los datos reales de una señal mas rapida que la mitad de la velocidad de mi ADC por lo descripto por vos. Una de las tecnicas para logralo es hacer capturas de pagina sucesivas con retardo de disparo de trigger entre ellas. Solo para señales periodicas logicamente.

Sin embargo, tal vez me exprese mal, no es que con matematica voy a hacer mas rapido el bicho electricamente, sino que puedo construir en la grafica en pantalla mas puntos de los que realmente dispongo, ejemplo, el punto que existe entre dos muestras sucesivas. Por supuesto que esto seria un "invento" ya que no es un sample real pero la matematica tiene que dejarlo muy cerca.

Pongo por ejemplo el osciloscopio DSO-5200A cuyo minimo valor de base de tiempo es de 2ns por division, o sea 2ns / 10. Un tiempo entre samples de "solo" 2ns significaria un ADC de 500mhz y eso sin contar que tenes 10 samples en 1 division. Por lo tanto, aca la matematica esta implicita (eso, o muestras con retardo de trigger pero no me imagino la precision del trigger para lograr tales velocidades), Decanto que tiene que ser matematica porque un osciloscopio relativamente barato como este no puede tener tremendo hardware.

Asi como puedo generar una onda senoidal para ser reproducida en un DAC a partir de samples, reconstruyendola con sen(x)/x como me indicas... ¿no sera posible reconstruirla de la misma forma en pantalla es en definitiva lo que busco?

Scooter: Lo que decis tiene mucha verdad, una cuadrada reconstruida a partir de formulas de senoidal no va ser buena, sin embargo podria intentar implementar un analisis previo de la señal para detectar al menos la forma basica y aplicarle las formulas que mejor ajuste al tipo de señal en la reconstruccion. (No pretendo que pueda ver una señal mas rapida que el ADC y que encima quede perfecta! sino aproximada)

Palurdo:

Hay cosas que mencionas que se escapan de mi conocimiento matematico, como derivadas, pero, me has dado datos que me sirven para googlear (que desconocia nombre o existencia) y que me pueden ayudar a entender un poco mas las ecuaciones para lograr algo.

En uno de los links que pase (http://www.arachnoid.com/polysolve/) el creador hizo funciones de programacion que generan los terminos calculados de un polinomio de grado P y a partir de alli hace el calculo de un punto inexistente X en la tabla de datos de entrada.

Bajo este criterio, yo pensaba implementar eso, dandole como tabla de entrada una cierta cantidad de muestras reales obtenidas para fabricar el polinomio que mejor ajuste a la señal que se esta muestreando en ese momento y a parti de ahi, calcular los puntos "fantasma" intermedios. De toda la pagina de samples que forman una pantalla, uso por ejemplo la mitad de ellos (x0, x1, x2, xn/2) y calculo los intermedios. Entonces, si pudiera ver una señal de 4mhz claramente con el ADC en modo "normal"

Muchas gracias por el material! a estudiar :) Si tienen mas info les agradezco mucho.
 
Creo que mas o menos ya te contestaron todos.

Lo ideal seria usar la funcion sinc, sen(x)/x, pero computacionalmente me imagino que es pesadita de calcular (aunque podrias usar tablas precalculadas). Otro problema que tiene sinc es que muy probablemente tengas muy feas deformaciones al comienzo y al final del conjunto de muestras. Me parece que para los efectos practicos una interpolacion polinomica que pase por cuatro puntos vecinos a la muestra deberia estar bien.
 
Bajo este criterio, yo pensaba implementar eso, dandole como tabla de entrada una cierta cantidad de muestras reales obtenidas para fabricar el polinomio que mejor ajuste a la señal que se esta muestreando en ese momento y a parti de ahi, calcular los puntos "fantasma" intermedios. De toda la pagina de samples que forman una pantalla, uso por ejemplo la mitad de ellos (x0, x1, x2, xn/2) y calculo los intermedios. Entonces, si pudiera ver una señal de 4mhz claramente con el ADC en modo "normal"

Cuidado con eso, porque el problema de los polinomios de interpolación, es que a mayor grado del polinomio, mayor error de interpolación cuando el punto a interpolar se aleja del punto del medio. Es muchísimo mejor interpolar muchas veces con pequeños polinomios cogiendo grupos de 3 (con derivada del punto central mejor) o 4 puntos (con o sin derivadas en puntos centrales), que interpolar una vez todos los puntos intermedios con un polinomio muy grande.

De todas formas implementar SINC(x) es relativamente fácil en el dominio de la frecuencia. Pongamos que tienes 512 puntos. Haces la FFT y tienes un espectro de 512 frecuencias complejas, añades otros 512 puntos de valor 0 en la parte alta de la FFT y entonces tienes una FFT de 1024 frecuencias, la mitad superior son 0. Entonces haces la FFT inversa y ya tienes tu señal duplicada mediante la función SINC.
 
Me suena que usé un osciloscopio con algo así; habían filtros básicos para aplicar; senoidal, cuadrada, triangular... con "cuatro samples" salía una onda perfecta. Lo malo es que era igual de perfecta si te equivocabas de interpolación.
 
Pregunto.. Mas allá de como hacer la interpolacion senoidal o lineal... Vos lo que intentas hacer, no tendrá algo que ver con muestreo en tiempo equivalente? (googlea "equivalent time sample rate")
 
Necesitas hacerlo en 2 pasos... primeramente tienes que calcular los coeficientes para generar un polinomio que cubra los puntos que quieres "amplificar" y despues evaluar ese polinomio en los puntos intermedios que quieres obtener

http://www.sc.ehu.es/sqwpolim/MATEMATICASII/cal.ppt
http://es.wikipedia.org/wiki/Interpolación_polinómica
http://es.wikipedia.org/wiki/Interpolación_de_Taylor

Te recomiendo que uses interpolaciones de Taylor ya que son mas faciles de calcular en un microprocesador, ya que permiten aproximar cualquier funcion usando e las 4 ecuaciones basicas (suma, resta, multiplicacion y division)

si el procesador es mas potente tambien puedes usar minimos cuadrados..

http://es.wikipedia.org/wiki/Ajuste_de_curvas
http://www2.dis.ulpgc.es/~lalvarez/...nsparenciasTema7_InterpolacionFunciones_2.pdf
 
Última edición:
Ahhh... la interpolación era para dibujar la onda en la pantalla!!!!
Te sugiero que atiendas la recomendación de palurdo y de Chico en cuanto a intepolación por tramos, por que no podés buscar calzar una parva puntos con un solo polinomio, por que el error va a aumentar y te vas a demorar en calcular el sistema de ecuaciones (va a ser un poco "grande"). Tirale la interpolación por 3 o 4 puntos y dibujá con eso.... debería andar muy bien en la medida de que el error del ajuste polinómico sea menor que la cuantización de los pixels en la pantalla (vamos... pixels "y medio" no hay :D).... y a fin de cuentas, siempre va a ser un "punto inventado".
 
.......
Tambien vi que osciloscopios digitales comerciales tienen una opcion (por defecto activada) de interpolacion sen(x) / x . Busque acerca de esto pero no pude o no supe encontrar la informacion de como se aplica (no veo variable Y en esa ecuacion)

¿Alguien sabria orientarme como puedo hacer esta interpolacion? ¿Como podria aplicarla? Soy de profesion programador asi que eso no es problema, solo me falta la matematica para hacer la funcion!!

Esa interpolación se basa en que una señal que no posea componentes de frecuencia superior a la de Nyquist puede expresarse en el dominio temporal de forma exacta por:

[LATEX]y(t)=\displaystyle\sum_{k=-\infty}^{\infty}{\dfrac{\sin(\omega(t-k\delta))}{\omega( t-k\delta)}\,y(k\delta)}[/LATEX]
siendo [LATEX]\delta[/LATEX] el intervalo entre muestreos.

Obviamente, la sumatoria no se hace entre +/- infinito, tomando 10 puntos para cada lado ya da buen resultado (más de 20 para cada lado es malgastar CPU porque las mejoras son insignificantes).

Como truncar la sumatoria ("filtro rectangular") acarrea el problema de oscilaciones en la gráfica --> efecto no sólo desagradable sino que genera confusión pues esas oscilaciones no existen en la señal, los coeficientes de la sumatoria se ponderan con funciones "de suavizado", que nos atenúan un poco las frecuencias altas, pero eso se compensa con más términos.
Hay diferentes funciones: Ventana rectangular, triangular, Hamming, Von Hann, Black... cada una con sus ventajas y desventajas. A mi gusto, la más general es la de Hamming.

El hecho que la expresión final de cada coeficiente resulte un bolonki no es ningún inconveniente, ya que son valores que se calculan una sola vez o pueden estar tabulados. Lo único que termina haciendo la rutina de graficación es multiplicar N pares de valores y sumarlos.
 
Bueno, pues te aconsejo que implementes varias formas de interpolado y permitas desde un menú elegir el tipo de interpolación, porque cada una tiene sus ventajas e inconvenientes. A veces la solución mejor y más simple es unir mediante rectas los puntos a dibujar en pantalla.

Sin embargo, prácticamente todas las interpolaciones presentan problemas en las discontinuidades de las formas de onda, lo que se conoce como efecto "gibbs", ya que las interpolaciones esperan construir una función de onda continua (suave y sin saltos).

Prueba la forma convolutiva que te dice eduardo, aunque es más eficiente hacer la FFT (enventanando la serie para minimizar el efecto gibbs, aunque distorsionas las frecuencias) añadir 0 en las frecuencias altas y hacer la FFT inversa. Más que nada porque la convolución tiene orden N^2 y la FFT tiene orden N*Log(N) que es más rápida.

De todas formas para micros sencillos es mucho más rápido una interpolación por polinomios pequeños, tipo splines.

Lo de la derivada no es problemático de entender. Imagina que por un punto de una gráfica pasa una recta. Pues la derivada de la recta en ese punto es la inclinación o pendiente de esa recta con respecto al eje X. Si el valor es muy alto, la recta estará muy hacia la vertical, y si es muy bajo, el valor será casi horizontal. Eso es la derivada de una función en un punto, la inclinación de la curva alrededor del punto donde se calcula la derivada.

Como para que una sucesión de curvas interpoladas sean suaves, la inclinación de la curva de la derecha del punto tiene que ser la misma que la inclinación de la curva a su izquierda, de lo contrario tendrías un pico o escalón en la unión de ambas curvas.

Por eso aproximamos la inclinación de un punto, por la pendiente de la recta que cruzaría ese punto y su anterior (por 2 puntos sólo puede pasar una recta), que al tratarse de pasos de 1 en 1, para el punto Y[1], su aproximación de derivada (inclinación que pasa por la recta de los puntos X[1],Y[1] y X[0],Y[0] sería (Y[1]-Y[0])/(X[1]-X[0])=(Y[1]-Y[0])/1=Y[1]-Y[0] (las X muestras están igualmente espaciadas en el tiempo, por lo que podemos suponer que los pasos valen 1 y nos simplifican el cálculo).

Podemos elegir 2 puntos consecutivos cualquiera de la secuencia de muestreo, aproximar sus derivadas como he comentado (por lo que hará falta un tercer punto anterior a Y[0], que será Y[-1]), y al tener los valores de los dos puntos y las dos derivadas ya tenemos las condiciones para resolver los cuatro coeficientes de la ecuación A*X^3+B*X^2+C*X+D=Y, para el intervalo de X=0...1, por ejempo podemos asignar a X, 10 valores que fueran 0 0,1 0,2 0,3 ..., y evaluar el polinomio en esos valores para calcular Y[0],Y[0,1],Y[0,2]... que serían 10 valores interpolados dentro del segmento de dos puntos.

Los coeficientes son, usando la variable auxiliar T:

T = Y[1]+Y[-1]-2*Y[0]
A = -T
B = 2*T
C = Y[0]-Y[-1]
D = Y[0]

En la imagen que adjunto te pongo la demostración del cálculo de los coeficientes:

Calculo Polimomio Cubico.jpg

Y una vez sabiendo como calcular valores intermedios entre dos puntos de una serie usando polinomios cúbicos, se puede hacer una función que te interpole un vector de M segmentos, cada segmento en N puntos, por lo que tendrías a la salida un vector de M*N elementos.

He implementado la función en Matlab y te darás cuenta de lo sencilla que es para lo que hace:

Código:
function salida = splines(samples, npuntos)

  Y_1=samples(length(samples));
  X=((1:npuntos)-1)/npuntos;

  for i=1:(length(samples)-1)

	Y0=samples(i);
	Y1=samples(i+1);

	T=Y1+Y_1-2*Y0;
	A=-T;
	B=2*T;
	C=Y0-Y_1;
	D=Y0;
	
	salida(i,:)=A*X.^3+B*X.^2+C*X+D;

	Y_1=Y0;
	end
 
 salida=salida';
 salida=salida(:)'; 

 return
end

A la función se le pasa un vector de samples, y un número de interpolaciones entre 2 samples, de manera que si por ejemplo el vector samples tiene 25 muestras y npuntos es 10, el vector resultado será un vector interpolado de 250 puntos.

Te pongo unos ejemplos usando la función de interpolación de arriba, de cómo funciona de bien (o de mal) este sistema. Los circulitos verdes son los puntos del vector de entrada , y la interpolación es sobre 10 puntos calculados por cada segmento de 2 puntos del vector de entrada, y unidos por segmentos rectos (en lugar de 10 podían ser 100 puntos, y no cambiaría nada ya que el polinomio de interpolación no depende del número de puntos en el que se calcule dentro del intervalo de interpolación).

Vector de cuatro puntos aparentemente sin relación:

4 Puntos interpolados.jpg

Puntos que dibujan una onda senoidal:

Seno interpolado.jpg

Un diente de Sierra:

Diente de Sierra.jpg

Onda Triangular:

Triangular.jpg

Onda Cuadrada:

Cuadrada.jpg

Como ves, para 4 puntos, la curva que los une es muy suave (no es la única curva posible, pero es una buena aproximación). Aunque no son perfectos, el Seno y la onda Triangular son bastante aceptables, porque no tienen saltos ni discontinuidades serias. La onda más problemática es el diente de sierra, donde en la discontinuidad aparece un bucle inferior hasta que se recupera el ascenso de la onda, bucle que también aparece aunque más atenuado en la onda cuadrada.
 
Excelente, les agradezco enormemente lo que estan brindandome, voy a detallar un poco mas aunque creo que todos captaron exactamente lo que quiero lograr y me estan ayudando aunque para comprenderlos al 100% me falta estudiar.

El osciloscopio que fabrique (se viene la publicada aca cuando me organice el despelote que es) tiene 2 canales, con 2 ADC corriendo en distintas bases de tiempo seleccionables: La mas rapida samplea a 40mhz, es decir, 25ns por sample. La memoria de almacenamiento es tipo fifo de 910 bytes, de los cuales uso 896 bytes debido a que los envio en paquetes de 64 bytes hacia la PC luego de terminado el proceso de sampleo de la señal. Tengo 10 divisiones en pantalla de 100 pixeles cada una, vale decir: 1000 samples y la diferencia entre lo que dispongo de memoria y la pantalla se completa con el valor 128, o sea, 0 volts (la onda "termina" antes en pantalla pero es totalmente util y tolerable ese detalle)

Dicho esto, el osciloscopio tiene entonces base de tiempos en: 2.5us/DIV, 5us/DIV y asi sucesivamente hasta los 360us/DIV a partir de los cuales el retardo entre samples ya no se hace por hardware sino por software (el PIC)

Entonces, para dominios de tiempos menores es donde quiero implementar esto. Tambien se me ocurrio lo siguiente:

Si tengo 1 sample cada 25ns y cada sample representa 1 pixel en pantalla, entonces si yo en vez de usar los 896 samples de un canal utilizo los primeros 448 pero en el eje X voy saltando de a 2 pixels (haciendo que una division use 50 samples en vez de 100, pero los muestre como 100) lo que estaria logrando es tener una base de tiempos "virtual" de 1.25uS /DIV

Hasta ahi perfecto, pero a medida que me voy yendo mas abajo en la base de tiempos, esta separacion entre puntos de sample se hace mas y mas visible y es por eso que quiero inventar matematicamente los puntos intermedios.

Cabe aclarar que este bicho esta funcionando de pelos! me sorprende lo que logre estirando los componentes mas alla de sus especificaciones de datasheet (por ej la memoria daria hasta 27mhz y se la banca de 10 a 40mhz!! no probe mas alla porque el ADC es de 40mhz)

Adjunto solo para compartir un par de imagenes del programa mostrando una senoidal de 20khz en canal A y una de 10khz en canal B. Una de ellas a 20us/DIV y la otra a 2.5us/DIV (maximo por hardware)

Señal mas rapida no puedo tirar porque la placa de sonido de la compu no puede ir mas alla asi que cuando termine el software me voy al taller de una amigo a darle caña con el generador de señales.

Y aclaro tambien que el procesamiento lo hace la compu asi que no tengo las limitaciones del micro del osciloscopio

Volviendo al thread y las matematicas: La tecnica de muestreo equivalente la probe generando un retraso de trigger de 83.33ns para el tiempo entre samples de 25ns. Funciona bien pero solo para el canal que tenga el trigger asociado, en el otro canal veo cualquier cosa porque la señal esta desenganchada del trigger. Ademas, solo logro 1 retraso posible en la configuracion de hardware actual entonces no podria pasar de 1.25uS / DIV.

Chico3001, gracias por la data y los links! El micro es el de la compu asi que tengo toda la potencia de calculo que haga falta :)

Como dice palurdo "A veces la solución mejor y más simple es unir mediante rectas los puntos a dibujar en pantalla." Si, la verdad que lo pense hacer asi, totalmente simple. A medida que voy mas abajo en la escala de tiempos tendria cada vez mas aliasing. La verdad tendria que hacer una prueba rapida en el programa a ver que tan "feo" se ve, quiza es tolerable. Me meti en esta solucion matematica para aprender mas que nada. Estoy viendo ecuaciones que asustan debido exclusivamente a mi falta de formacion en esa area, pero las ire masticando poco a poco hasta entenderlas.

Palurdo, lo que me brindas es algo muy apreciado y les prometo que antes de seguir con el thread me voy a sentar con un buen cafe a estudiar cada uno de los links y las respuestas que me han dado todos para comprender y aplicar esto al programa. Aprecio enormemente el tiempo que han invertido en ayudarme y algo bueno voy a sacar de todo esto!

Digo "antes de seguir con el thread" porque como aprecio lo que estan haciendo, no quiero quitarles el tiempo con preguntas antes de procesar esta informacion primero :)

Esto promete! con un poco de tiempo voy a armar el thread de este osci y explicando los 7 intentos fallidos previos para que todo el que quiera armar uno con cosas que andan por ahi sueltas pueda hacerlo. (La mayoria de los componentes son recicle de placas relacionadas a computadoras, excepto los ADC y los OP)
 

Adjuntos

  • 20-10khz-a-2.5usdiv.png
    20-10khz-a-2.5usdiv.png
    74.1 KB · Visitas: 23
  • 20-10khz-a-20usdiv.png
    20-10khz-a-20usdiv.png
    76.9 KB · Visitas: 23
Palurdo,

Nunca use mathlab, trato de bajar el trial pero no me lo dan :( no importa, me ayudo con un pequeño tuto que encontre de mathlab e intento implementar tu funcion mathlab en vb.NET para poder entenderla viendo que sucede con ella dados distintos parametros de entrada. Aproximadamente entiendo lo que hace pero, abusando de tu ayuda voy a hacerte un par de preguntas:

Las mismas las pongo en contexto de la funcion mathlab asi se entiende que es lo que pregunto, en estilo comentario C // :)

Código:
function salida = splines(samples, npuntos)

// Aca entiendo que asignas a Y_1 el valor del ultimo de los samples del vector inicial
  Y_1=samples(length(samples));

// En la siguiente linea, no comprendo la sintaxis 1:npuntos, parece ser la inicializacion de un vector de
// rango de 1 a 9 en caso de que npuntos = 10. Esto dividido por la cantidad de puntos, ahi me perdi
  X=((1:npuntos)-1)/npuntos;

// Aqui recorro el vector de samples desde el primero al anteultimo
  for i=1:(length(samples)-1)

	Y0=samples(i);  // Obtengo el valor del punto actual del ciclo for
	Y1=samples(i+1); // Y el siguiente punto

	T=Y1+Y_1-2*Y0; // Esto parece ser simple matematica con los valores, no con vector
	A=-T; // Parece ser que A toma el valor T *-1
	B=2*T; // B=T*2
	C=Y0-Y_1; // C = valor sample actual - ultimo sample?
	D=Y0; // D= sample actual
	
        // Aca aplicas el polinomio pero segun vi en un tutorial, el "X." se refiere a una matriz o 
       // vector, no a un dato simple?
	salida(i,:)=A*X.^3+B*X.^2+C*X+D;
        // en salida(i;:) estas especificando que el resultado del calculo se almacena en la posicion I del
        // vector de salida? o quiza aqui es donde el resultado del calculo NO es un dato simple sino un
        // vector, y lo asignas A PARTIR de la posicion i y que ocupe la cantidad de elementos del vector 
        // resultado. Siendo asi, me perdi porque el bucle FOR solo va de 1:totalsamples-1. O tal vez, en 
        // salida no es un vector sino una matriz multidimensional donde i es el puntero y luego tengo, en
        // este ejemplo, 10 datos en esa posicion.

	Y_1=Y0; // Asignando sample actual a Y_1 para la proxima vuelta.
	end
 
// Aqui no comprendo la asignacion de salida prima, pero supongo que es parte de la sintaxis de mathlab
 salida=salida';
 salida=salida(:)'; 

 return
end

Te pido disculpas si consideras que estas preguntas estan fuera de lugar. Estoy de todos modos tratando de conseguir un mathlab para comprenderlo. Aparentemente la parte de la funcion donde aplicas el calculo trabaja con matrices, multiplicando una matriz o vector y elevando a potencia de 3, etc.etc. Si asi fuera entonces como primera medida tengo que ver como traducir lo que es tan simple en mathlab a lenguaje "normal" que no tiene multiplicacion de matrices como operandos directos.

Si pudieras contestar esas dudas estoy muy agradecido, si no podes, no hay problema porque de todos modos me has dado (todos ustedes) una base muy grande desde donde trabajar. Hasta esta hecho! :)
 
Pues en realidad has entendido bastante bien el funcionamiento del programa, lo único es que, lo entiendo, parece el tema un poco lioso por las operaciones con matrices, pero en realidad no se están haciendo operaciones matriciales, sino operaciones elemento a elemento. Es una forma de evitar en matlab hacer un bucle for añadido, además de que el programa queda mucho más sencillo.

Te doy comentado el código (con los comentarios al modo de matlab):

Código:
function salida = splines(samples, npuntos)

% Y_1 ya sabes que significa Y[-1] pero no se puede poner el signo - en una variable en matlab
% El valor inicial de este elemento lo hago el último sample porque este valor sería el que iba antes
% del primero. Yo asumo que la señal es periódica y cuando acaban los datos vuelven a empezar 
% pero podía haber puesto Y_1=0 y con una ligera variación de la inclinación inicial, hubiera servido igual.
  Y_1=samples(length(samples));

%Efectivamente esto crea un vector de N elementos que van de 0 a 1. Digamos que son las coordenadas
%en las que se parte el intervalo de 0 a 1 en N trozos. Si N=10 X=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9] 
%Si N=4 X=[0 0.25 0.5 0.75] creo que ahora lo ves claro.

  X=((1:npuntos)-1)/npuntos;

% Aqui recorro el vector de samples desde el primero al anteultimo
  for i=1:(length(samples)-1)

	Y0=samples(i);  % Obtengo el valor del sample actual del ciclo for Y[0]
	Y1=samples(i+1); % Y el siguiente sample Y[1], ya teníamos de antes el sample anterior en Y[-1]

	T=Y1+Y_1-2*Y0; % Con la formula para sacar los coeficientes A,B,C,D
	A=-T; % Salen directos a partir de los samples Y[-1],Y[0],Y[1].
	B=2*T; 
	C=Y0-Y_1; % C = valor sample actual - valor sample anterior
	D=Y0; % D= sample actual
	
        %Se aplica el polinomio sobre todos los valores del vector X para guardar el resultado en
        %la coluna i, conteniendo esta columna el mismo número de elementos que X. A,B,C,D son números
        %y cuando un número multiplica un vector, multiplica todos sus elementos a la vez. Y en matlab
        %si pones un . ante un operador, como es el operador potencia ^, entonces matlab interpreta que
        %no quieres hacer una potencia matricial, sino que quieres una potencia numérica sobre todos los
        %valores de la matriz (en este caso es un vector columna).
	salida(i,:)=A*X.^3+B*X.^2+C*X+D;
        %Tienes razón, salida se convierte en una matriz de length(samples) filas por npuntos columnas.
	Y_1=Y0; %el sample actual Y[0] será el sample anterior Y[-1] en la siguiente vuelta.
	end
 
%La prima sirve para transponer una matriz, es decir, lo que eran filas pasan a ser columnas y viceversa.
salida=salida';

%al poner salida(:) transformo la matriz en un vector columna (uniendo las columnas en orden) y al
%poner la prima, paso la matriz columna a una matriz fila, es decir, a un vector normal. 
salida=salida(:)'; 

 return
end

Si no entiendes algo, no te cortes en preguntarmelo. Es muy fácil traducir este algoritmo a cualquier otro lenguaje, pero matlab permite jugar con esta función a través de su consola de comandos. Busca cómo dibujar gráficas, con las funciones plot, axis, hold, etc. En la consola de matlab, cuando quieres ayuda sobre un comando, sólo tienes que poner help [comando] por ejemplo para dibujar gráfica "help plot"



Por cierto, la versión GNU de Matlab y supuestamente compatible con él, es Octave, lo que no sé si habrá versión para windows.
 
Última edición:
Bueno, te confirmo que con GNU Octave bajo Windows funciona la función igual exactamente que en Matlab.

Para instalar la función en Octave, la guardo en un archivo llamado splines.m y ese archivo lo copio dentro de una carpeta que he creado y he llamado proyectos, dentro de la ruta de octave, en mi caso 'C:\Software\Octave-3.6.4\proyectos'

Te voy a poner 2 ejemplos hechos en la ventana de comandos de Octave:

Código:
GNU Octave, version 3.6.4
Copyright (C) 2013 John W. Eaton and others.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  For details, type `warranty'.

Octave was configured for "i686-pc-mingw32".

Additional information about Octave is available at http://www.octave.org.

Please contribute if you find this software useful.
For more information, visit http://www.octave.org/get-involved.html

Read http://www.octave.org/bugs.html to learn how to submit bug reports.

For information about changes from previous versions, type `news'.

 - Use `pkg list' to see a list of installed packages.
 - MSYS shell available (C:\Software\Octave-3.6.4\msys).
 - Graphics backend: gnuplot.

octave-3.6.4.exe:1> % antes que nada compruebo que tengo el
octave-3.6.4.exe:1> % archivo splines.m con mi función, copiada
octave-3.6.4.exe:1> % dentro de la carpeta proyectos donde voy
octave-3.6.4.exe:1> % a trabajar.
octave-3.6.4.exe:1> cd proyectos
octave-3.6.4.exe:2> dir  
.          ..         splines.m
octave-3.6.4.exe:3> % Defino samples a mi antojo
octave-3.6.4.exe:3> Y=[1 1 2 2 3 3 4 4 1 1 2 2 3 3 4 4 5 5 6 6 5 4 3 2 1]
Y =

 Columns 1 through 19:

   1   1   2   2   3   3   4   4   1   1   2   2   3   3   4   4   5   5   6

 Columns 20 through 25:

   6   5   4   3   2   1

octave-3.6.4.exe:4> % Preparo las coordenadas X para plot.
octave-3.6.4.exe:4> % interpolaremos 25 puntos por segmento.
octave-3.6.4.exe:4>
octave-3.6.4.exe:4> X=25*(1:length(Y))-24
X =

 Columns 1 through 13:

     1    26    51    76   101   126   151   176   201   226   251   276   301

 Columns 14 through 25:

   326   351   376   401   426   451   476   501   526   551   576   601

octave-3.6.4.exe:5> % el primer elemento comienza en la posicion 1
octave-3.6.4.exe:5>
octave-3.6.4.exe:5> % ploteamos los length(Y) puntos
octave-3.6.4.exe:5> length(Y)
ans =  25
octave-3.6.4.exe:6> % los plotearemos en color rojo como circulos
octave-3.6.4.exe:6> plot(X,Y,'ro');
octave-3.6.4.exe:7>
octave-3.6.4.exe:7> % mantenemos el contenido de la grafica para nuevos plots
octave-3.6.4.exe:7> hold on
octave-3.6.4.exe:8> % plotearemos las 25x25=625 interpolaciones en azul
octave-3.6.4.exe:8> plot(splines(Y,25),'b-'); % unidas por lineas
octave-3.6.4.exe:9>
octave-3.6.4.exe:9> % borraremos el contenido de la grafica para nuevos plots
octave-3.6.4.exe:9> hold off
octave-3.6.4.exe:10>
octave-3.6.4.exe:10> %creemos otro vector
octave-3.6.4.exe:10>
octave-3.6.4.exe:10>
<10 +10 -9 +9 -8 8 -7 7 -6 6 -5 5 -4 4 -3 3 -2 2 -1 1 0 0]
U =

 Columns 1 through 15:

   -1    1   -2    2   -3    3   -4    4   -5    5   -6    6   -7    7   -8

 Columns 16 through 30:

    8   -9    9  -10   10   -9    9   -8    8   -7    7   -6    6   -5    5

 Columns 31 through 40:

   -4    4   -3    3   -2    2   -1    1    0    0

octave-3.6.4.exe:12> % ahora en lugar de 25, usaremos 100 puntos interpolados
octave-3.6.4.exe:12> V=100*(1:length(U))-99;
octave-3.6.4.exe:13> 
octave-3.6.4.exe:13> % Imprimiremos los puntos en la grafica como antes:
octave-3.6.4.exe:13> plot(V,U,'ro');
octave-3.6.4.exe:14> 
octave-3.6.4.exe:14> %congelamos la grafica otra vez
octave-3.6.4.exe:14> hold on
octave-3.6.4.exe:15>
octave-3.6.4.exe:15> %imprimimos los puntos interpolados como antes:
octave-3.6.4.exe:15> plot(splines(U,100));
octave-3.6.4.exe:16>

Y aquí sus gráficas:

Ejemplo 1:
Interpolacion1.png

Ejemplo 2:
Interpolacion2.png

Instalate Octave que es gratis, desde sourceforge (versión compilada con virtual studio) y juega con la función splines para ver cómo funciona...
 
Bueno, pude implementar la funcion mathlab en el programa vb.net del osciloscopio con muy buenos resultados!

Ante todo, el nucleo de la cosa, que era:

Código:
salida(i,:)=A*X.^3+B*X.^2+C*X+D;

Se transformo en esto (a lo chancho nomas sin optimizaciones):

Código:
For ix = 0 To NPuntos - 1
  dRet = A * X(ix) ^ 3 + B * X(ix) ^ 2 + C * X(ix) + D
  If dRet < 0 Then dRet = 0
  ret(p) = CByte(dRet)
  p = p + 1
Next

Este bucle for esta dentro del for principal de la funcion, Y donde la variable ret() es un vector del tamaño final (evito hacer traspasos de matrices, etc) y pasa derecho al graficado.

La señal original en canal A (verde) de 20khz y en canal B (amarillo) de 10khz, la medi en una base de tiempos alta de 80uS / DIV para poder tener una senoidal "apretada" que pudiera luego expandir como si fuese bajando la base de tiempos (todo esto por no tener un generador mas alla de 20khz).

Esa señal la almacene en el archivo que adjunto "testdata.txt" donde los primeros 896 samples son canal A y los restantes canal B.

896 samples por canal a 80us / DIV escalados a 100 samples / DIV dan lo siguiente:

20khz80usdiv_real.png

Y ahora es donde la funcion de palurdo empieza a trabajar. Lo siguiente es una muestra de los primeros 896 / 2 samples de cada canal, con una interpolacion de npuntos = 1 (o sea, intercalar un punto calculado entre 2 reales, como si bajase la base de tiempos a 40us/DIV)

interpolacion1.png

Adjunto 2 interpolaciones mas: una con npuntos = 2 y otra con npuntos = 5. Recien con npuntos = 127 aproximadamente se empieza a deformar de forma notable la senoidal. Esto representa una multiplicacion considerable de la señal que puede graficar llegando la base de tiempos a cerca de 20ns / DIV

Por lo tanto, podria decir que esto es mucho mas de lo que esperaba sacar de este bicho :) Entiendo perfectamente que no estaria viendo transitorios o cosas por el estilo pero para ser algo casero va muy bien!

Sin embargo, atento a la idea de, en vez de usar matematicas para calcular puntos intermedios, tomar el sample nro 1, dibujarlo en el pixel X = 0, luego tomar el sample nro 2 y dibujarlo en el pixel X = 4 (seria como npuntos = 5) uniendolos por una linea recta, tambien se obtiene un muy buen resultado! La imagen es la siguiente:

salteando5.png

Esto ultimo resulta extremadamente mas liviano porque no involucra calculos extra de ningun tipo. Por supuesto eso no es un problema en una PC.

Asi que ahora estoy en la disyuntiva, creo que voy a aplicar la funcion de interpolacion como una opcion de menu para los casos donde simplemente unir con lineas rectas los samples no alcance.

Todo esto me abre a nuevas posibilidades, la implementacion de esta funcion no la veo complicada de hacer directamente en un PIC, con lo cual se podria hacer un osciloscopio portatil usando el ADC lento del micro y que "invente" los samples intermedios para bases de tiempo rapidas. (Misma mecanica: Almacenar a todo trapo en ram, luego calcular y luego graficar)

Obviamente que teniendo los ADC de 40mhz y todo ya hecho no voy a volver atras, pero es algo interesante para experimentar mas adelante.

Mas alla de todo esto y que el problema principal esta resuelto, toda la documentacion que me han brindado me abrira mas perspectivas y otras posibles formas de hacerlo.

Por ultimo, les adjunto una foto del engendro. Fue hecho en placas simple faz con SMD del lado cobre y DIP del otro lado, en forma de "islas" en cada integrado unidos con puentes y rodeado con plano de masa. Todo cable de interconexion excepto los que sean digitales y de señal no cambiante (los de control) es mallado con un extremo de la malla conectado a masa y el otro libre, pasando los conductores por adentro.

En la foto se ve la fuentecita que transforma los 5v del usb en 5v -5v y 2.5v, se ve el pic al lado del cable usb, se ve el oscilador de cristal de 80mhz y junto a el el generador de base de tiempos. En la parte analogica solo se ven aqui los 3 mux, una llave 4066 y los 4094 los cuales 1 es de control (me quede sin pines en el pic) y el otro es el generador de nivel de voltaje para el trigger. Los operacionales y un flip-flop del trigger estan del otro lado en smd.

Muchas gracias a todos por los valiosos datos que me han brindado!
 

Adjuntos

  • interpolacion2.jpg
    interpolacion2.jpg
    43.7 KB · Visitas: 6
  • interpolacion5.png
    interpolacion5.png
    57.7 KB · Visitas: 6
  • testdata.txt
    8 KB · Visitas: 5
  • hardware.jpg
    hardware.jpg
    90.7 KB · Visitas: 25
pue espera que todavia se puede mejorar mucho (vamos a limpiar esas senoidales jeje)

En cuanto tenga un rato te sigo explicando. Por lo demas enhorabuena. Saludos
 
Bueno, pues ya te puedo explicar. La función hace una interpolación hacia atrás, es decir, dada la muestra que conocemos en el instante 1, usa las muestras anteriores 0 y -1 para obtener el polinomio interpolador. Esto es útil en tiempo real, cuando se va interpolando a la vez que se captura cada sample. Esto es lo que se conoce como una función "causal" puesto que considera tiempos pasados y presentes. Aquella que considera tiempos futuros, como la muestra 2 que todavía no se conocería en tiempo real, es "anticausal" y en general no causal. La idea para mejorar la interpolación es obtener dos interpolaciones simétricas y promediarlas de forma que los errores de una cancelen los de la otra. Usando la misma función se puede hacer procesando 2 veces, la primera la usual para generar el vector causal, la segunda, invertir el orden de los datos de entrada en la función, y luego después generar el vector interpolado volver a invertir el orden de los datos de salida, para obtener el vector anticausal. Una vez obtenidos ambos vectores, se promedian (Vc+Vac)/2 y se obtiene la interpolación tanto hacia adelante como hacia atrás.

Lo que pasa es que hacer esto es un jaleo. Es mejor en lugar de promediar los vectores, calcular el polinomio resultante de promediar los 2 polinomios, el causal y el anticausal. Ahora necesitamos también tener en cuenta el valor de Y2 para calcular la pendiente del punto 1.

Nos ayudamos con una nueva variable auxilar S además de la T, y la fórmula queda ligéramente modificada:

S=Y[2]+Y[0]-2*Y[1];
T=Y[1]+Y[-1]-2*Y[0];

A=(S-T)/2;
B=(2*T-S)/2;
C=(Y[1]-Y[-1])/2;
D=Y[0];

Tan simple como esto, y la interpolación mejora sustancialmente.

Pongamos una muestra de 6 puntos de una onda senoidal. Si la representamos con rectas queda así:

Senoidal 6 puntos.png

Con la interpolación causal (la primera que se ha implementado antes) queda así:

Senoidal 6 puntos interp.png

Con la interpolación causal-anticausal, la senoidal queda así:

Senoidal 2.png

¿Te has dado cuenta qué bonita queda, pues con más puntos queda todavía más bonita?

Aquí te pongo el código de la función nueva, que he llamado splines2.m:

Código:
function salida = splines2(samples, npuntos)

  Y_1=samples(length(samples));
  X=((1:npuntos)-1)/npuntos;

  for i=1:(length(samples)-1)

	Y0=samples(i);
	Y1=samples(i+1);
	if (i==length(samples)-1) 
 	  Y2=samples(1);
	else 
          Y2=samples(i+2);
	end

	S=Y2+Y0-2*Y1;
	T=Y1+Y_1-2*Y0;

	A=(S-T)/2;
	B=(2*T-S)/2;
	C=(Y1-Y_1)/2;
	D=Y0;
	
	salida(i,:)=A*X.^3+B*X.^2+C*X+D;

	Y_1=Y0;
	end
 
 salida=salida';
 salida=salida(:)'; 

 return
end

Compárala si quieres con la primera y verás que sólo está la diferencia en el cálculo de los coeficientes del polinomio interpolador, lo demás es lo mismo.

Te adjunto las imágenes de los cálculos de interpolación de las gráficas de los post anteriores, pero recalculadas con la función de interpolación nueva:

Cuatro puntos interpolados:
4 puntos interpolados.png

Onda senoidal (15 puntos):
Onda Senoidal.png

Diente de Sierra:
Diente de Sierra.png

Triangular:
Triangular.png

Cuadrada:
Cuadrada.png

Ejemplo 1:
Ejemplo 1.png

Ejemplo 2:
Ejemplo 2.png

Como ves, la senoidal se queda prácticamente perfecta, la curva de los 4 puntos se ajusta más a los puntos que la otra, es más suave, el diente de sierra y la onda cuadrada, aunque se nota que todavía deforman, pero han mejorado mucho y ya parecen un diente de sierra y una cuadrada. La triangular aunque tiene todavía algo de deformidad en las puntas, la deformación ha mejorado también bastante y además ahora es simétrica, dando la impresión de una triangular de verdad. Los ejemplos que te he pasado antes también ajustan mejor a los puntos.

Bueno, pues ya me contarás si introduces estos pequeños cambios en tu función para ver qué tal quedan las ondas en el osciloscopio.

Un saludo.
 
Atrás
Arriba