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



## seaarg (Jul 23, 2013)

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


----------



## Dr. Zoidberg (Jul 23, 2013)

No entiendo lo que querés hacer  
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.


----------



## Scooter (Jul 23, 2013)

Sí sabes que la señal que estás midiendo es de un tipo podrías interpolar, lo que pasa es que no lo sabes...


----------



## palurdo (Jul 23, 2013)

seaarg dijo:


> 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.
> 
> ...



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.


----------



## seaarg (Jul 23, 2013)

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.


----------



## chclau (Jul 23, 2013)

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.


----------



## palurdo (Jul 23, 2013)

seaarg dijo:


> 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.


----------



## Scooter (Jul 23, 2013)

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.


----------



## hazard_1998 (Jul 23, 2013)

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")


----------



## Chico3001 (Jul 23, 2013)

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


----------



## Dr. Zoidberg (Jul 24, 2013)

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 ).... y a fin de cuentas, siempre va a ser un "punto inventado".


----------



## Eduardo (Jul 24, 2013)

seaarg dijo:


> .......
> 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.


----------



## palurdo (Jul 24, 2013)

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:



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:


```
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:



Puntos que dibujan una onda senoidal:



Un diente de Sierra:



Onda Triangular:



Onda Cuadrada:



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.


----------



## seaarg (Jul 24, 2013)

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)


----------



## seaarg (Jul 24, 2013)

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 // 


```
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!


----------



## palurdo (Jul 24, 2013)

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):


```
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.


----------



## palurdo (Jul 24, 2013)

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:


```
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:


Ejemplo 2:


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...


----------



## seaarg (Jul 25, 2013)

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:


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

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


```
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:



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)



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:



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!


----------



## palurdo (Jul 25, 2013)

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


----------



## palurdo (Jul 25, 2013)

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í:



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



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



¿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:


```
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:


Onda senoidal (15 puntos):


Diente de Sierra:


Triangular:


Cuadrada:


Ejemplo 1:


Ejemplo 2:


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.


----------



## palurdo (Jul 25, 2013)

Bueno, pues para hacer una prueba con datos reales, he cogido los primeros 60 puntos de tu fichero testdata, y los he decimado por un factor de 5 (es como si hubieras generado una onda de 20kHz*5=100kHz.

En la gráfica represento los samples en rojo, sus uniones mediante rectas en verde:

Para la primera función de splines, corresponde a la curva azul (interpolada a 25 puntos por segmento aunque el aspecto sería prácticamente el mismo si fueran 100 puntos):



Para la función de splines nueva, corresponde a la curva negra, igual que antes, interpolada a 25 puntos por segmento:



Como ves, la mejora con respecto a la función antigua y a la unión por rectas es sustancial, sobre todo al ampliar mucho la forma de onda en la pantalla.


----------



## seaarg (Jul 25, 2013)

Excelente! me estas brindando mucho mas de lo que esperaba y lo agradezco de verdad. Actualmente le deje el codigo de la funcion splines implementado al programa. Esta noche cuando me ponga con mis ediciones nuevamente voy a aplicar lo que aqui me describis.

Como comentario: Por alguna razon, estimo que electrica, el primer dato de las series de samples para cada canal tiene un valor con mas "salto" que el resto de los samples, se me hace que es alguna inestabilidad en el momento que el trigger se dispara.

Dicho esto, en la muestra de samples que adjunte antes ese dato existe (porque grabo raw data como viene por usb) y por eso en tu curva ves la primera parte como mas "recta" pero en el programa del osciloscopio, descarto unas 3-4 primeras muestras (no recuerdo el numero pero si que descarto)


----------



## seaarg (Jul 26, 2013)

Bueno, me he tomado un tiempo de leer el resto del material que me han dado y por esos enlaces encontre un PDF que me aclaro algunos otros conceptos. http://www2.dis.ulpgc.es/~lalvarez/...nsparenciasTema7_InterpolacionFunciones_2.pdf

Intente aplicar lo que me especifico Eduardo sobre la funcion sinc







Sin embargo se me escapan algunos conceptos. Veo claramente que es
(sen(valor) / valor) * valordeunsample. Todo esto en sumatoria de los 10 samples (o calculos) anteriores y 10 posteriores

Lo que se me escapa creo es que no sabria que es ω. Entiendo que multiplica al valor de "tiempo" - la constante de entrada K (supongo los -10+10) * 1, siendo 1 si realizara una interpolacion de 1 elemento (intervalo entre muestreos? o te referis a tiempo?)

En fin, entonces aplique la 2da version de la funcion de palurdo con los siguientes resultados.

Señal de 20khz generando 16 puntos cada 2 samples, uniendolos solo con segmentos (linea de sample 1 a 2, o sea, este caso no genera puntos sino que saltea 16 pixels y dibuja linea entre ellos)



Lo siguiente es la nueva funcion de palurdo con Y2



Como se puede ver, hay un error al inicio del grafico. Sobre este error, estoy casi 100% seguro que se debe a que por alguna razon que aun no encontre, los primeros datos del vector de puntos estan incorrectos (tal vez este tomando en el canal B amarillo al inicio, los ultimos datos del canal A)

Entonces, implemente un descarte de los primeros 10 datos haciendo que el vector en X sea X+10 y completando el vector (por cuestiones de programacion y limite de array) con ceros antes de pasarlo a la funcion de interpolacion.

El resultado es este:



Como ven, va perfecto, solamente que al final, como tengo un cambio brusco a 0 (que no es real) la funcion tira cualquier cosa.

Me toca averiguar porque razon tengo ese inicio de grafica raro. Tal vez sea que no estoy parando la captura en el lugar que yo creo (tiempos de pic) o que estoy contando mal algun dato en los vectores del programa. De hecho, en la primera imagen donde no hay matematica sino simple union de puntos, se nota claramente que hay un dato del canal amarillo que no esta bien al inicio.

Tambien me resulta curioso que no haya demasiada diferencia visual entre unir con segmentos y calcular realmente los puntos. Sin embargo en las graficas que pone aqui palurdo se nota claramente la diferencia. Es posible que la diferencia se empiece a notar en señales de mayor amplitud, donde los saltos de Y entre un sample y otro son mayores.


----------



## palurdo (Jul 27, 2013)

seaarg dijo:


> Tambien me resulta curioso que no haya demasiada diferencia visual entre unir con segmentos y calcular realmente los puntos. Sin embargo en las graficas que pone aqui palurdo se nota claramente la diferencia. Es posible que la diferencia se empiece a notar en señales de mayor amplitud, donde los saltos de Y entre un sample y otro son mayores.



En realidad se empieza a notar en señales de mayor frecuencia. Cuando la distancia de los puntos es corta (hay muchos puntos para dibujar una onda), las inclinaciones de las rectas de los puntos contiguos es bastante parecida, por lo que dos rectas que se unen con una ligera inclinación, parecen una curva. Es por eso por lo que te dije que para señales de baja frecuencia con muchos samples simplemente uniendo con rectas es suficiente (Yo te diría que para una senoidal de 60 puntos no hay mucha diferencia entre unir con rectas y unir interpolando), ya que la inclinación de la curva en el punto 1 (por decir un punto cualquiera de la muestra) es muy similar a la del punto 2, por lo que la recta es una buena aproximación. 

El problema es cuando tienes pocos puntos por periodo de señal, entonces la inclinación entre un punto y otro de la gráfica es muy alta y es cuando el margen de error de unir con rectas se puede apreciar bastante. Uno de los casos extremos es el que te presenté cuando la onda tiene sólo 6 puntos (en el caso de frecuencias cercanas a nyquist, un periodo de onda sólo tiene un par de puntos). En ese caso unir con rectas no da una idea de la onda senoidal, sin embargo interpolando puedes ampliar todo lo que quieras la onda, ya que seguirá pareciendo senoidal.

Es por eso por lo que decimé tus datos por 5, para simular una onda senoidal de 100kHz, o lo que era lo mismo, de 12 samples por periodo. Si hubiera dejado la onda original (los 60 samples) habría obtenido el mismo resultado que tu, una onda muy similar interpolando y uniendo con rectas, ya que con 60 puntos  por periodo senoidal las dos aproximaciones son buenas.

Haz la prueba, en lugar de dibujar todos los samples, de dibujar sólo 1 de cada 10 samples en pantalla y unirlos con rectas (a 2X dibuja el sample 0 en posición 0, el 10 en posición 20, el 20 en 40...), y después interpla por 10 (o a 2X por 20), y verás el cambio entre una y otra gráfica. En este caso estarías haciendo como si hubieras muestreado una senoidal de 200kHz.



Te voy mostrar un ejemplo:

Dos señales senoidales de duración 28 muestras. La primera de ellas tiene un periodo de 27 muestras, lo que con 28 muestras se dibuja la onda en la pantalla completamente. En la siguiente gráfica dibujo, en la parte de arriba, la onda unida con rectas y en la parte de abajo, la interpolada. Aunque la interpolada es más suave y en la de arriba se notan trozos rectos de segmento, la diferencia más o menos es como en tu gráfica (onda de 60 puntos que al ampliar por 2 es como si tuvieras una onda de 30 puntos).



Ahora la segunda senoidal. La longitud de la muestra es lo mismo, 28 muestras. La amplitud es la misma, entre -1 y 1. Pero el periodo de la senoidal es un periodo de 5,4 muestras, por lo que cada 27 muestras, cuando en la primera senoidal teníamos una onda en la segunda tenemos 5 ondas.



La imagen habla por sí misma. Cabe destacar que las variaciones entre los puntos consecutivos son mayores cuanto mayor es la frecuencia de la señal (es decir, menor es el periodo), aunque la amplitud total de la onda sea la misma a una frecuencia o a otra.


----------



## palurdo (Jul 27, 2013)

Bueno, como supongo mi función de interpolación puede ser útil para ser implementada en un procesador de baja potencia, voy a tratar de optimizar la ecuación de interpolación para utilizar aritmética de números enteros (mucho más rápida que la aritmética de punto fijo y punto flotante), así mismo tratar de evitar todas las operaciones de multiplicación y división que sean posibles, a fin de que el algoritmo de evaluación del polinomio interpolador sea eficiente en una CPU o MCU que no tenga implementado por hardware la multiplicación y la división, como podría ser un PIC, por ejemplo.

Para ello, vamos a transformar la expresión matemática, tanto de los coeficientes del polinomio interpolador como la del polinomio en sí, por una expresión equivalente.

En aritmética de enteros, es importante aplazar las operaciones de división a las últimas etapas del cálculo, ya que si se ha de multiplicar dos enteros, la multiplicación siempre será entera, pero si se han de dividir 2 enteros, la división no siempre será entera, y entonces habría que descartar la parte decimal, es decir, introducir un error. Por lo tanto cuanto más temprano se haga una división más error saldrá en el resultado final de la fórmula. 

Por otro lado, nos interesa que las divisiones y multiplicaciones sean del tipo 2^u, porque multiplicar un número N por 2^u es desplazar N hacia la izquierda u bits, y dividir N por 2^u, desplazar u bits hacia la derecha (con lo cual no se están haciendo multiplicaciones y divisiones verdaderas sino sólo bit-shift).

Así vamos a definir unos coeficientes alternativos a los A,B,C,D, pero que sean números enteros, a los que llamaré A',B',C' y C':



Ahora iremos a la ecuación original Y=AX^3+BX^2+CX+D. Resulta que tenemos un problema, ya que X tiene que ir desde el rango 0..1, es decir que X son todo decimales, y nosotros buscamos hacer cálculos enteros así que tenemos que transformar nuestras X a algo que podamos usar con aritmética de enteros.

Pongamos que queremos partir un segmento en cuatro segmentos. El primer punto empieza en X=0, el segundo en X=0,25, el tercero X= 0,5 y el cuarto en x=0,75. Pero es que podemos ver que es lo mismo que X[0]=0, X[1]=1/4, X[2]=2/4, X[3]=3/4. En general, si partimos un segmento en un número potencia de 2, como sería 2^n, tendremos que X[0..(2^n)-1]=(0..(2^n)-1)/2^n. Si partimos en 4 trozos, estamos partiendo en 2^2 trozos. Si partimos en 16 trozos, en 2^4 trozos, 64 trozos->2^6, y 256 trozos->2^8.

Hasta aquí queda claro que podemos partir en X[0..1] en 2^n trozos mediante (0..(2^n)-1)/2^n. Entonces la nueva X iría desde 0 hasta 2^n-1, sustituyendo en la ecuación:



Como se ve, en la última ecuación hemos utilizado la sustitución de los coeficientes por sus primas, de esa manera la única división que aparece en todo el cálculo es el denominador de la fracción, que además es del tipo 2^(3n+1), lo que significa que para dividir el numerador, sólo hay que desplazar hacia la derecha 3n+1 bits. Los bits decimales que resultaran de ahí se descartan ya que no nos interesan ya que si representan pixeles a dibujar en pantalla, no hay pixeles fraccionales.

Otra cosa es que en pantalla se esté haciendo un zoom vertical(pongamos un zoom de 2^p), por lo que en lugar de desplazar hacia la derecha 3n+1 bits, se podrían desplazar simplemente 3n+1-p bits y de nuevo descartar decimales.

Ahora bien, el polinomio de arriba, con las X variando de 1 en 1 desde 0..2^n-1, tiene varias exponenciaciones y multiplicaciones. Los multiplicadores del tipo 2 elevado a algo es trivial simplemente desplazando hacia la izquierda tantos bits como sea ese algo. Ahora, las exponenciaciones hay un método muy bueno para evitarlas y transformarlas en multiplicaciones sencillas:



Simplemente agrupando y sacando factores comunes, acabamos por eliminar las potencias de X. La operación es sencilla. Se va guardando en la variable de resultado las operaciones que van desde los parentesis interiores hacia los exteriores hasta que se completa el polinomio. 

De de manera que ahora sólo nos hace falta multiplicar por X, pero ¿como lo hacemos en un procesador que no dispone de multiplicación?, no queda más remedio que hacer una multiplicación por software ,pero OJO, lo importante es el rango de las X que nos va a definir cómo de compleja será esa multiplicación, con respecto a las sumas y desplazamientos necesarios.

Vamos a hacer números para saber cuantos bits necesitamos de memoria para cada variable y para el resultado, teniendo en cuenta varias consideraciones:

-El rango de los datos es de 0 a 255 (8 bits). 
-Vamos a interpolar como máximo 2^6 puntos entre 2 muestras reales (es decir, sacar 64 puntos entre 2 muestras. Luego en la fórmula el índice de la X recorrerá desde 0 hasta 63 para cada pareja de puntos a interpolar.
-Luego calcularemos lo que necesitamos si decidimos un rango de 2^8 puntos.

Bueno, pues sabiendo esto, y haciendo unos cálculos, tenemos que 

max(A')=1020; max(B')=1530;max(C')=255;max(D')=510
min(A')=-1020; min(B')=-1530;min(C')=-255;min(D')=0

Por lo que necestamos 11 bits para A', 12 bits para B', 9 bits para C' y 10 bits para D' (en teoría para D' son 9 bits también, pero consideramos más util poder operar como si min(D') fuera -510.

Por consideraciones prácticas, para cada coeficiente asignamos 2 bytes de memoria (16 bits).

Ahora en la expresión fórmula, para resolver el polinomio numerador, necesitamos saber el tamaño del resultado del numerador. Resulta que el término que domina es A*X^3, ya que X^3 se hace pronto mucho más grande que X^2 y X^1, aunque B sea mayor que A. Sin embargo el hecho de estar multiplicados por 2^(u*n) donde u=0,1,2,3, iguala todos los coeficientes en rango, por lo que habría que hacer la suma de A'+B'+C'+D' máximos y multiplicar por 2^(3n) para saber cuantos bits nos hacen falta. Así necesitamos 13+3n bits. Para el caso de 64 particiones, necesitamos 31 bits en total, que podemos empaquetar en 4 bytes (considerando enteros con signo en complemento a 2). Es decir, tenemos que implementar sumas de 4 bytes (no es muy complicado el código ya que son 4 sumas con arrastre en cascada byte a byte). Para el caso de 256 particiones necesitamos 37 bits, por lo que las sumas serán en paquetes de 5 bytes (40 bits).

Por otro lado, X va a ir desde 0 a (2^n)-1, es decir, que para n=6 cuando se multiplica por n, habría que implementar la multiplicación con 6 sumas y 6 desplazamientos como máximo (si X es 0, son 0 sumas y 0 desplazamientos). Si n es 8 (particiones de 256 interpolaciones), las multiplicaciones serán 8 sumas y 8 desplazamientos.

Es decir, que ahora yendo a la última ecuación, para n=6, se efectuan 3 multiplicaciones por X, 3 desplazamientos de bits, y 3 sumas. Además un desplazamiento hacia la derecha para efectuar la división de 2^(3n+1) que será un desplazamiento de 19 bits, lo que nos deja un resultado final de 31-19=12 bits posibles para el valor de interpolación (11 bits mas el signo, es decir, que la onda podría ir en el peor de los casos hasta el valor 2048 y en el peor a -2048, aunque por lo general se quedará dentro del intervalo 0-255). Pero cada multiplicación como máximo son 6 sumas y 6 desplazamientos, por lo que tenemos en total que para calcular un punto interpolado se necesitan :

6*3+3=21 sumas de 4 bytes y 6*3+3+1=22 desplazamientos. Supongamos que todo el proceso nos ocupa 43*4(bytes por número)=172us en un pic a 4MHz, podríamos calcular aproximadamente unos 5800 puntos interpolados cada segundo (sin contar los cálculos de los coeficientes de los polinomios).

Bueno, espero que estos cálculos les sea útiles a alguien, porque puede venir muy bien para implementar un osciloscopio LCD portátil con un procesador que no sea muy potente.

Un saludo.


----------



## cosmefulanito04 (Jul 27, 2013)

Interesante el tema.

¿Qué micro usás que el ADC puede trabajar a 40MHz?


----------



## seaarg (Jul 27, 2013)

Palurdo:

En tu primer mensaje se ve claramente la diferencia y es muy considerable. Lo que voy a hacer es ponerme a corregir el programa del micro para poder hacer un descarte efectivo de las primeras muestras. Tengo 910 bytes de memoria por canal de los cuales estoy usando 896 asi que bien puedo descartar, digamos una potencia de 2 al inicio (porque la cantidad de samples "corruptos" parece depender de la escala de tiempo) y usar samples del final, de esos 910 bytes, que no estoy usando. Esto haria que se me corra el trigger algunos uS en el tiempo pero es tolerable ya que de por si la mecanica que uso produce este efecto. El trigger dispara la captura automatica y luego el micro se "da cuenta" que disparo el trigger y temporiza el final de la captura. Hasta que el micro se entero de que el trigger disparo pasan un minimo de 160nS aprox y eso puede estar corrompiendo los datos iniciales.

Sobre el segundo desarrollo, esta muy bueno para implementarlo como decis en osciloscopio portatil. El micro mio corre a 48mhz o sea, 83.33ns por instruccion asi que bien se podria bancar ese algoritmo sin hacerse muy lento, pero disponiendo de la potencia de calculo de la PC, lo dejo para una version portatil con LCD y ADC mas lentos.

cosme: Bienvenido al thread! nos hemos cruzado en varios de estos de osciloscopio asi que se que el tema te interesa. Para responderte, el micro es el PIC18F2550 corriendo a 48mhz, pero el ADC no es el del micro sino un TLC5540 de texas instruments que puede correr hasta 40mhz. Entre el micro y el adc hay una memoria FIFO D42101G. A pesar de que el datasheet de la misma dice que sus ciclos son de 34ns, se banca perfectamente ciclos de 25ns para escritura.

El osciloscopio es de 2 canales duplicando adc y memoria. Bien se podria haber implementado lo mismo con solo 1 conjunto adc-memoria y haciendo un distribuidor de señal para los 2 canales pero como disponia de ellos, use 2 para simplificar. Es automata, es decir, el micro le da la orden al circuito de ponerse en modo captura y cuando la señal de trigger aparece el circuito funciona solo (ahi el micro esta a la espera de fin de captura nada mas). Cuando termina la captura, el micro detiene la secuencia y se pone a leer las 2 memorias por el puerto B, para transferirlas a la PC. En el circuito de captura, hay un mux digital que selecciona la base de tiempos de la siguiente forma:

Oscilador 40mhz -> contador 4040 para hacer division de tiempo -> mux digital para elegir la base de tiempos.

Cuando esto no es suficiente (desde los 320us/DIV para arriba) el micro toma el control y temporiza el mismo la base de tiempos (ahi hubo que contar instrucciones generadas de C a ASM)

Entonces, tengo divisiones en escala de tiempo de: 2.5us, 5us, 10us, 20us, 40us, 80us, 160us, 320us, 500us, 1ms, 2ms, 10ms, etc.etc.

El algoritmo de interpolacion se activa a partir de valores de escala menores a 2.5us, con 1 punto de interpolacion tengo 1.25us/DIV y asi sucesivamente.

En cuanto a la parte analogica, dispongo de selectores AC/DC, divisores 1:1, 1:2, 1:5 y multiplicador 1x5 para los dos canales + 1 trigger externo.

El trigger lo implemente con 2 flip flops integrados en el 74F74. Uno esta asociado al selector de nivel de trigger (que es un DAC hecho con un 4094 y una escalera r2r) y el otro tiene una base de tiempos RC para que dispare automaticamente el trigger si no se produce un evento de trigger real (trigger mode auto, queda desenganchado de la señal)

El flip flop actua sobre la compuerta NAND que habilita el circuito de reloj de captura.

Como reflexion, la parte digital la verdad que se ve compleja pero es sencillisima. Lo que fue un desafio para mi fue la parte analogica. Segun mediciones contrastadas con un osciloscopio comercial DSO-5200A funciono muy bien!

La desventaja del sistema de trigger (precisamente, con el nivel de trigger) es que tengo de 0 a 5volts por lo tanto en escalas de division se me reduce el rango del mismo.
Tambien le voy a implementar un divisor 1:2 en la entrada que se active automaticamente cuando selecciono escalas de mayor vpp que -5+5 (en este momento, si pongo mas de 10vpp en la entrada tengo recorte de señal, logico.)

Otro errorcito de diseño es que hice que la referencia del ADC sea de 5v, lo que me da 256/5v de resolucion. Tendria que haber implementado una referencia de 2v y hacer que la entrada adapte el nivel a fin de tener mejor resolucion vertical.

Por ultimo, los operacionales de entrada son LMH6644 lo que me da un ancho de banda analogico teorico de 130mhz pero ese valor es para señal de salida de 200mvpp. Un dato interesante es que se consiguen muy facilmente y baratos los operacionales LM318 que son de unos 15mhz lo cual es interesante.


----------



## palurdo (Jul 27, 2013)

La verdad es que tu proyecto es muy interesante seaarg. Voy a disfrutar del diseño cuando lo tengas acabado y lo publiques. Una pregunta. Cómo haces para muestrear? me refiero al sistema de Sample&Hold que utilizas. Creo haber leido por ahí que usas un 4066 por algún lado. Qué fin tiene en tu circuito?. Un saludo y espero tus novedades.


----------



## seaarg (Jul 27, 2013)

Gracias!

Me voy a volver loco para pasar los esquematicos jeje mas que nada por la cantidad de lineas de bus de datos entre las 2 memorias hacia el micro. El programa livewire no se si tiene la opcion para hacer un bus.

Adjunto aca lo que seria el nucleo del sistema de captura:



Representa conexiones desde el micro y trigger hacia el control de memoria y ADC.

RST = Reset de memoria, tanto escritura como lectura. Recordemos que son memorias FIFO que permiten grabar por un lado y leer por otro a la vez o a distintos tiempos y en vez de tener lineas de direccion, solo tienen un clock de escritura y otro de lectura que hacen que cambie la direccion secuencialmente. Como ven, conecta directamente 1 pin del micro a los reset de lectura y de escritura de las 2 memorias.

SCK = Señal de control con doble funcion: Por un lado, es el clock de lectura cuando transfiero los datos de la memoria hacia el PIC, y por otro lado, es el clock de grabacion para bases de tiempo mayores a 320 us/div cuando el pic es el que controla el tiempo entre muestras. Se observa que esta señal controla el clock del ADC (ADC x CK) y el clock de grabacion de las memorias WCK X (negado con respecto al adc)

Cuando estoy en modo grabacion automatica (base de tiempo de 320us/div hacia abajo) este pin controla la habilitacion o no del sampleo. Es decir, deja de ser un clock para ser un pin de control.

25mhz = Dice 25 porque originalmente usaba un oscilador de 50mhz fijo, pero en realidad esto no es una señal de esa frecuencia, sino que es de 40mhz, 20mhz, 10mhz, 5mhz, 2.5mhz y asi sucesivamente, dependiendo de la base de tiempo seleccionada. Como conte mas arriba, esta señal de reloj esta formada por un oscilador fijo de 80mhz que entra al clock de un contador 74HC4040 y sus salidas desde Q0-Q7 son las distintas frecuencias mencionadas. Estas pasan por un mux digital para seleccionar la que nos interese en la base de tiempos.

FCK EN = Habilitacion o no de la toma de samples automatica (fast clock). Esta señal no viene del micro (originalmente si) sino que viene de un flip-flop del sistema de trigger. Es decir, para hacer una pagina de samples, el micro pone SCK = 1 y habilita el trigger para que dispare la señal FCK EN, dejando pasar la frecuencia de la base de tiempos al ADC y memoria. Transcurrido un tiempo determinado por la base de tiempos, el micro pone SCK = 0 y el control de trigger = 0 para detener la captura.

~WE = No mucho que decir: Habilita la escritura de memoria. Este pin esta constantemente en cero.

RCK = Read clock: Clock de lectura de memoria, se usa solo en transferencia memoria-pic.

~RE = Como veran, uso la ultima nand que me queda en el 74F00 como inversor. Esto es simplemente el selector canal A, canal B para la lectura de la memoria hacia el pic. (En el momento de escritura las memorias se escriben a la vez)

De a poquito voy a ir armando el esquematico completo de todo esto para compartir con todos ustedes. Es sencillo solamente que un poco largo y mis habilidades con los programas de esquematicos son limitadas  El osciloscopio lo hice mayormente directo en PCB wizard, usando un poco el proteus para simular algunas condiciones. (entrada analogica, trigger,etc pero la parte de captura digital la hice "en la cabeza")

Respondiendo concretamente a tu pregunta, el sample & hold lo hace solo el ADC que mencione mas arriba. Solamente hay que darle un clock y el hace su trabajo.

La llave 4066 no interviene en la captura, sus funciones son:
- Seleccionar o no trigger canal B / Externo (el canal A se selecciona desde un mux que mezcla la fuente de trigger)
- Seleccionar o no el auto-trigger (que es un flip flop que dispara el trigger automaticamente si no viene el disparo normal de trigger desde un operacional usado como comparador)

Les adjunto tambien un zip con los 3 pcb que forman este osciloscopio. No creo que les sirvan de nada sin los esquematicos pero si alguien quiere pegarle una mirada, ahi estan.

Por ultimo: Buenas noticias! Esta tarde implemente un descarte de las primeras 6 muestras y un retraso de algunos nanosegundos (dependientes de la escala de tiempos) antes de detener la captura y ahora las ondas estan perfectas, ya no tengo el salto que habia al inicio o al final de la pagina de samples, por lo tanto, ahora puedo hacer una mejor implementacion de los algoritmos brindados por palurdo sin que esos valores corruptos afecten el polinomio. De hecho, en un par de pruebas pude ver la diferencia entre los segmentos y el algoritmo y si, definitivamente la curva es suave.

Ahora tengo que hacer algunos numeros para resolver lo siguiente:

- Dado el tiempo minimo entre muestras de 25ns, y usando 100 pixels por division, tengo 2.5us / div como base de tiempos minima. Usando el algoritmo de interpolacion, empiezo a bajar esa base de tiempos. Si intercalo solo 1 sample calculado entre 2 reales, tengo 1.25us / div, pero cuando intercalo 3 samples calculados entre 2 reales, tengo 625ns/div y esos numeros no me gustan mucho. Quiero hacer cuentas para ver cuantos samples tengo que intercalar para hacer, digamos, en vez de 625ns/div: 500ns/div, 200ns/div, etc. (numeros mas redondos)

Una pregunta aparte: Como uso 896 samples por canal, a 100 pixels por division en una pantalla cubro 8.96 divisiones. Si yo quisiera hacer que el programa simplemente calcule los pixeles restantes para llegar a 1000 samples (o sea, las 10 divisiones normales  de cualquier osci) esta funcion polinomica no me sirve verdad? Supongo que no debido a que necesita 2 samples para calcular los intermedios y cuando X > 895 ya no dispongo de mas samples. Se que podria hacerlo de alguna otra forma pero no lo tengo muy claro. Alguna sugerencia?

Otra: Una vez perfeccionado el soft en la parte de interpolacion, planeo que el mismo me de datos utiles, como ser, una etiqueta que me diga de que frecuencia es la señal. Se me ocurrio que, mientras voy recorriendo los samples, calcular el valor minimo y maximo: Luego, cuando tengo esos datos, contar cuantos samples pasan desde un minimo pasando por un maximo y hasta un minimo de nuevo. Ahi podria calcular la frecuencia. Tendria que poner algo con un error de un cierto % debido a que esos minimos y maximos bien pueden ser picos en la señal y no la señal en si. No me gusta... es muy "fuerza bruta". ¿Hay algo mejor para hacer esto?

Disculpen! me entusiasme escribiendo!


----------



## seaarg (Jul 27, 2013)

Bueno, mas resultados!

Logre encontrar la causa de los datos iniciales corruptos. Me habia olvidado que el ADC tiene un pipeline de 4 muestras (a pesar de ser flash) por lo que cualquier lectura antes de los 4 clocks debe tomarse como invalida. Entonces, modifique el firmware para hacer algunas lecturas previas que no se envian a la PC.

Solucionado esto, le meti una señal de 20khz en ambos canales pero en una escala lenta de 500us/div.

Esto es lo que se ve:



Entonces, para tener una referencia de como deberia verse la interpolacion, capture una pantalla de la misma señal pero a 40us/div (aclaro: no se va a ver igual porque utilizare interpolacion 16 y no da exactamente 40us/div sino mas abajo, en relacion a 500us/div)

Esta es la señal:



Vuelvo la escala de base de tiempos a 500us/div pero esta vez, aplicando interpolacion 16 simplemente uniendo con 1 segmento entre samples: (salteo 15 pixels)



Y ahora, la interpolacion 16 pero con la funcion polinomica: (fabrico 15 puntos en Y)



Si comparo con la señal original en 40us/div, como pueden ver esta muy parecida! y las curvas son suavecitas exceptuando al final donde no tengo mas samples entonces la funcion utiliza un cero para Y2. Esta perfecto.

Toca empezar a probar con señales cuadradas a ver que pasa, me fabricare un oscilador con un inversor y un cristal de 4mhz como fuente de señal. Ademas, como dije antes, me queda aun mucha matematica que implementar para que el programa me de los datos relevantes de la señal: Frecuencia, vmin, vmax, vpp, etc. asi como tambien hacer escalas de tiempo un poco mejores, mas redondas. Y luego de tener un soft bastante funcional, toca hacer los esquematicos para compartir por aqui.

Muchisimas gracias!


----------



## palurdo (Jul 30, 2013)

seaarg dijo:


> - Dado el tiempo minimo entre muestras de 25ns, y usando 100 pixels por division, tengo 2.5us / div como base de tiempos minima. Usando el algoritmo de interpolacion, empiezo a bajar esa base de tiempos. Si intercalo solo 1 sample calculado entre 2 reales, tengo 1.25us / div, pero cuando intercalo 3 samples calculados entre 2 reales, tengo 625ns/div y esos numeros no me gustan mucho. Quiero hacer cuentas para ver cuantos samples tengo que intercalar para hacer, digamos, en vez de 625ns/div: 500ns/div, 200ns/div, etc. (numeros mas redondos)



Bueno, pues este tema es muy interesante que lo propongas, porque gracias a que tienes polinomios interpoladores, puedes poner la base de tiempo arbitraria que necesites. Al tener la onda formada por una combinación de funciones continuas, puedes evaluar cada tramo de la función en el valor que necesites, no necesariamente una partición entera del intervalo original. Pongamos un ejemplo:

Si intecalas 4 samples entre 2 reales, tendrás una base de 500ns/div, ya que creas 5 segmentos por lo que partes un segmento de 2500ns/div entre 5. Esto equivale a ir recorriendo las X en pasos de [0 0,2 0,4 0,6 0,8 1] (los puntos 0 y 1 son reales, el resto son interpolados). Si te has dado cuenta, lo que hacemos es dividir la escala nueva por la escala vieja para obtener el paso necesario para pasar de un punto interpolado al siguiente, es decir, tenemos 500ns/2500ns=0,2. Ahora, si en lugar de 500ns queremos una escala de 1000, pues hacemos lo mismo, 1000ns/2500ns=0,4. Ahora los pasos para dibujar serán [0 0,4 0,8 1,2 1,6 2], donde 0 y 2 son puntos reales y los demás puntos interpolados.

¿Pero qué ha pasado con el punto 1?, pues que al quedarse entre medias no se imprime en pantalla, pero ya se imprimen los valores cercanos los puntos que están alrededor, es decir 0,8 y 1,2, por lo que la onda sigue fielmente representada en pantalla. De todas formas la función de interpolación sólo funciona en el intervalo 0 a 1, por lo que tenemos que partir el cálculo de los puntos interpolados [0,4 0,8 1,2 1,6] en 2 partes: La primera Y0 e Y1 serán el punto 0, y el 1, y los puntos a calcular serán 0,4 y 0,8. Y la segunda Y0 será ahora el punto 1 y Y1 el punto 2, y los puntos a calcular serán 0,2 (1,2-1) y 0,6 (1,6-1). Parece un poco lío, pero de forma algorítmica se resuelve simple y elegántemente con un bucle que vaya dibujando los pixeles en pantalla (o los guarde en un vector de gráfica para luego imprimirlos). 

De hecho el incremento de cada paso no tiene por qué ser un valor tan limpio, puede ser cualquier valor real, de manera que podemos hacer una opción para personalizar la base de tiempos a nuestro antojo, por ejemplo si el usuario requiere que la base de tiempo sea de 720ns, pues 720/2500 = 0,288, así la secuencia interpolada será [0 0,288 0,576 0,864 1,152 1,440 1,728 2,016 2,304...] (los cálculos serán de 0,288, 0,576, 0,864 sobre el intervalo 0...1,  0,152, 0,440 0,728, sobre el intervalo 1...2, 0,016, 0,304... sobre el intervalo 2...3 y así hasta completar los pixeles en pantalla.

Así podría definir el algoritmo en pseudo C de la siguiente manera:


```
double incremento=escalanueva/escalavieja,j=0,X=0;

for(int i=0;(i<MAX_PIXEL)&&(j<MAX_DATOS);i++) {

  X=trunc(j); //X será la parte entera de J.

  if (X==j) Dibujapunto(i,DATOS[X]); //si j es entero, dibuja Y0 en el punto DATOS[X];

  else {

    Y0=DATOS[X];     //Preparamos los coeficientes para pasarlos a la función
    Y1=DATOS[X+1];   
    Y2=DATOS[X+2];
    Y_1=DATOS[X-1];

    X=j-X; //ahora X es la parte decimal, ya que X debe de ir entre 0 y 1.

    Y_interpolado=splines(Y_1,Y0,Y1,Y2,X); //Evaluamos X en el polinomio interpolado.

    Dibujapunto(i,Y_interpolado);

  }

  j+=incremento;
}
```



seaarg dijo:


> Una pregunta aparte: Como uso 896 samples por canal, a 100 pixels por division en una pantalla cubro 8.96 divisiones. Si yo quisiera hacer que el programa simplemente calcule los pixeles restantes para llegar a 1000 samples (o sea, las 10 divisiones normales  de cualquier osci) esta funcion polinomica no me sirve verdad? Supongo que no debido a que necesita 2 samples para calcular los intermedios y cuando X > 895 ya no dispongo de mas samples. Se que podria hacerlo de alguna otra forma pero no lo tengo muy claro. Alguna sugerencia?



Lamentablemente, mientras interpolar da una buena idea de una función si esta no es una función muy irregular (por ejemplo una función fractal no es interpolable ya que es infinitamente discontinua), ya que el objetivo de interpolar es representar una curva suave, el calcular puntos fuera de un intervalo, lo que ya no es interpolar sino extrapolar, no es que sea dificil si no que la señal que obtienes extrapolando no es que no sea real (que lo es, al igual que los puntos interpolados dentro de un intervalo), es que ni siquiera lo parece ya que cuanto más te alejas de los puntos utilizados para calcular el valor irreal, más se desvía dicho valor de algo parecido a la función original. Por ejemplo puedes interpolar sin problemas una senoidal de 100 puntos, pero cuando vas a calcular los 50 puntos siguientes de fuera de la senoidal, te encuentras cualquier cosa menos la senoidal. 

Pues sobre ideas, una que no me gusta mucho, y otra que me gusta más. La que no me gusta mucho es suponer que al ser la señal periódica, reproducir una parte del principio a continuación del final. No me gusta esto porque cuando la señal no es periódica estás engañando a quien ve la captura. Por ejemplo, imagina que estamos viendo un tren de bits modulados. Cuando terminara esa zona de pantalla con los datos reales, lo siguiente sería o bien extrapolación o duplicado del principio, lo que confundiría a la hora de analizar los datos, ya que por ejemplo si los datos son 1 0 0 1 0 1 0 0 pero en pantalla solo se han capturado la representación modulada de 1 0 0 1, veríamos esto 1 0 0 1 1 0... lo cual no es la señal original y puede dar muchos dolores de cabeza a quien esté analizándola.

La idea que se me pasa por la cabeza es hacer la ventana más pequeña de manera que en lugar de que te falten datos te sobren, por ejemplo: haces una división cada 50 pixeles y cubres 500 muestras en 10 divisiones. El resto de samples los usas como buffer de pantalla, de manera que el usuario pueda recorrer la ventana mediante flechas o botones desde 0 hasta 895 samples (por ejemplo podría situar la ventana entre el sample 165 y el 665). ¿Que la ventana se queda muy pequeña?, no hay problema, haces zoom x2, división cada 100 pixeles e interpolas 1 punto entre muestras para tener en total 1792, y haces la ventana de 1000 puntos (a recorrer entre 0 y 1791).

Esta solución es la más elegante y la que creo que implementan muchísimos osciloscopios.



seaarg dijo:


> Otra: Una vez perfeccionado el soft en la parte de interpolacion, planeo que el mismo me de datos utiles, como ser, una etiqueta que me diga de que frecuencia es la señal. Se me ocurrio que, mientras voy recorriendo los samples, calcular el valor minimo y maximo: Luego, cuando tengo esos datos, contar cuantos samples pasan desde un minimo pasando por un maximo y hasta un minimo de nuevo. Ahi podria calcular la frecuencia. Tendria que poner algo con un error de un cierto % debido a que esos minimos y maximos bien pueden ser picos en la señal y no la señal en si. No me gusta... es muy "fuerza bruta". ¿Hay algo mejor para hacer esto?



Vamos a ver. Parámetros interesantes para calcular cuando tenemos una onda:

-Valor Máximo
-Valor Mínimo
-Frequencia principal
-Nivel de continua
-Valor eficaz (RMS).

Con el tema de la frecuencia, tienes razón, el que comentas es un método poco fiable para poder medir la frecuencia, ya que el hecho de que hayan máximos y mínimos locales (crestas y valles) no significa que la frecuencia sea la inversa de la distancia entre cresta y cresta. Podrías estar calculando sin querer la frecuencia de un armónico (la señal de normal será la suma de varios armónicos que hacen que la onda vaya subiendo y bajando muchas veces más que la frecuencia original).

Para tener una indicación fiable no hay más remedio que usar una transformación de la serie temporal de los datos a una serie frecuencial, es decir calcular la FFT. Aunque tienes 896 puntos, puedes calcular la FFT óptima de 1024 puntos ya que añadir los puntos de valor 0 que faltan no varía las componentes frecuenciales de la FFT. 

La FFT te va a calcular puntos en el plano complejo, cada punto tiene una parte real y otra imaginaria, pero a ti lo que te interesa es conocer el módulo de cada frecuencia (amplitud). Dado un número complejo, a+bi, su módulo es sqrt(a^2+b^2). Ya de entrada te digo que una vez calculada la FFT, la componente primera, es la de frecuencia 0, entonces esa primera componente es el valor medio, o la componente en continua, de la señal muestreada. Buscando un valor que nos interesa nos hemos encontrado con otro valor que también nos interesa.

Ahora, si tenemos 1024 puntos, y hacemos una FFT, la primera frecuencia es la componente que se repite cada 0 samples(la continua). La segunda es la que se repite cada FFT/2, la tercera se repite cada FFT/4, la cuarta cada 3/4 FFT, y así, hasta la componente 512 que es la frecuencia que se repite cada 2 samples (la frecuencia de nyquist). Del 513 al 1024 es la misma FFT que de 1 a 512, pero reflejada en orden descendente, por lo que no hace falta que te fijes en ese tramo... Así sabiendo la frecuencia de la FFT que te interesa, sabes cual es la frecuencia de la señal que necesitas. 

La frecuencia principal de una onda, quitando su valor de continua FFT[0], es la frecuencia de la componente de la FFT que tenga mayor amplitud. Así que para saber la frecuencia de tu señal, recorres las amplitudes de las 512 frecuencias de la FFT, y te quedas con la frecuencia de mayor amplitud. Así sabrás cada cuantos samples la onda se repite. 

El problema es que la FFT te va a dar una aproximación digital de la frecuencia, nunca una frecuencia que no sea racional (porque usando series discretas sólo podemos tener frecuencias racionales, es decir, frecuencias que se puedan expresar en una fracción de la frecuencia de muestreo). Cuando tienes muchos puntos, esta frecuencia digital será bastante próxima a la frecuencia real de la señal, (por ejemplo, pongamos que la frecuencia real es 3.141593 MHz, y al calcular la FFT de la señal muestreada, nos da una frecuencia de 3.142857 MHz..., que es el resultado de que la razón de la señal con respecto a la frecuencia de muestreo sea de 22/7. Entonces pongamos que cada 7 puntos la señal se repite en el análisis digital, pero en realidad se debería de ir repitiendo cada 7.0028174 puntos. Tienes ahí una pequeña fracción de 0,0028174 que te va a introducir un error en el cálculo de la frecuencia (esto es sólo una estimación, el error real entre el muestreo y la frecuencia real va creciendo conforme la onda va acercándose a la frecuencia de muestreo). 

Entonces podemos ayudarnos con las funciones polinómicas de interpolación para aproximar el punto X exacto donde la onda se repite. Ponemos que una onda vale 5 en el punto 0, y la FFT nos dice que su frecuencia principal es de 16 samples. Pues en el sample 16 deberíamos esperar encontrar un 5, pero nos encontramos con un 6. En el sample anterior, el 15, tenemos un 4,5. Entonces tenemos un intervalo, desde X=15,Y=4,5 hasta X=16,Y=6, y entre un valor de ,llamemosla Xu, de entre 15 a 16, debe de haber un valor de Xu donde Yu sea 5. Pongamos que podemos construir dentro de ese intervalo un polinomio interpolador que llamaremos P(X). Podemos construir una función F(X)=P(X)-Yu, de forma que F(X) sea 0 para el valor de Xu esperado. Entonces debemos aproximar X a Xu para conseguir que F(Xu)=0. Esto lo podemos hacer bien resolviendo directamente las soluciones del polinomio interpolador de tercer grado (un peaso de fórmula increiblemente grande, pero se puede hacer), o bien numéricamente en varios pasos con el método de la bisección, el de la secante, o el de newton (el de la bisección da buen resultado ya que tenemos un intervalo definido). Entonces encontramos una Xu=15,37 por ejemplo, la cual F(15,37)=0, por lo tanto P(15,37)=5= valor del sample 0. Entonces 15,37 samples será una mejor aproximación de la frecuencia que el valor de 16 samples que nos había dado la FFT.


Sobre los valores máximos y mínimos: Pues puedes recorrer todos los valores de la muestra e ir investigando si el valor actual es mayor o menor que el máximo y mínimo provisionales, hasta finalizar el vector de muestras. El problema de nuevo lo tenemos en que estamos midiendo unos valores máximos y mínimos sobre muestras discretas, y sin embargo disponemos de funciones continuas donde podemos calcular el máximo y el mínimo de cualquier tramo de la función, por lo que encontrar la cresta y el valle en esas funciones nos sería muy útil ya que casi nunca va a coincidir el valor de la muestra con el valor de cresta o valle de la onda. 

Bueno, pues resulta que para encontrar la cresta y el valle de un polinomio de tercer grado (el punto donde la pendiente de subida y bajada se hace 0), lo que hay que hacer es derivarlo e igualar la derivada a 0, de manera que obtenemos una ecuación sencilla de segundo grado que podemos resolver teniendo dos soluciones Xa y Xb, mediante la conocidisima fórmula de ecuaciones de segundo grado, pero adaptada a los coeficientes de la derivada de nuestro polinomio interpolador:






Esas soluciones existirán siempre (excepto cuando lo que multiplica a X^2 es cero, porque entonces tenemos solución única que será X=-c/b. Pero tanto Xa como Xb tienen que encontrarse dentro del intervalo 0 y 1 para que nos sirvan, de forma que si Xa y Xb vale por ejemplo -1,34 y 2,37, pues en el intervalo 0 y 1 no tiene un punto de valle ni de cresta, entonces para saber dentro de ese intervalo cual es el máximo o el mínimo, sólo hay que ver el valor en el punto 0 y el 1 y ver cual es mayor que el otro. En definitiva, para saber el máximo y el mínimo de un intervalo 0...1, tenemos considerar 4 puntos y ver cual es el mayor y cual es el menor dentro de ese intervalo, es decir X=0, X=Xa, X=Xb, X=1. Y luego ir recorriendo el vector hasta encontrar el valor máximo absoluto y el mínimo.

El programa lo he implementado y probado en Octave:


```
function [xmax ymax xmin ymin] = maxminsplines(samples)

  Y_1=samples(length(samples));

xmax=0;
ymax=0;
xmin=0;
ymin=0;

  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;
                  
        if  A==0
		X1=-C/(2*B);
		X2=X1;
	else
		X1=(-2*B+sqrt(4*B^2-12*A*C))/(6*A);
		X2=(-2*B-sqrt(4*B^2-12*A*C))/(6*A);
	end

	Y1=A*X1^3+B*X1^2+C*X1+D;
	Y2=A*X2^3+B*X2^2+C*X2+D;

	if X1>0 && X1<1 && imag(X1)==0
		if Y1 >ymax
			ymax=Y1;
			xmax=i+X1;
		end
		if Y1 <ymin
			ymin=Y1;
			xmin=i+X1;
		end
	end

	if X2>0 && X2<1 && imag(X2)==0
		if Y2 >ymax
			ymax=Y2;
			xmax=i+X2;
		end
		if Y2 <ymin
			ymin=Y2;
			xmin=i+X2;
		end
	end

	if D>ymax
		ymax=D;
		xmax=i;
	end 
	if (A+B+C+D)>ymax
		ymax=A+B+C+D;
		xmax=i+1;
	end
	if D<ymin
		ymin=D;
		xmin=i;
	end 
	if (A+B+C+D)<ymin
		ymin=A+B+C+D;
		xmin=i+1;
	end 
	Y_1=Y0;
	end
 return
end
```

Como ves, ahora ya no hace falta introducir en el programa el número de puntos a interpolar, ya que lo que se calcula es el valor máximo y mínimo absoluto en un X que no tiene que ser una partición fracional del intervalo 0..1 sino que se encontrará en cualquier parte de ese intervalo. Ahora tampoco se devuelve el vector de interpolación sino los valores X e Y máximos y mínimos.

Ejemplo, tenemos estos puntos representados en la siguiente gráfica:



Y para comparar con la función interpolada N=25 que sería esta:



Ahora usamos la función maxminsplines para encontrar los X e Y máximos y mínimos de los polinomios interpoladores de los samples y tenemos este resultado:


```
octave-3.6.4.exe:49> [xmax ymax xmin ymin]=maxminsplines(Y)
xmax =  12.574
ymax =  0.99667
xmin =  17.603
ymin = -0.99680

octave-3.6.4.exe:50> [y x]=max(Y)
y =  0.96350
x =  13
octave-3.6.4.exe:51> [y x]=min(Y)
y = -0.96781
x =  18
```

Es decir, sin interpolar, el máximo lo da en el sample 13, con valor 0,96350, y el mínimo en el sample 18 con el valor -0,96781. Usando los máximos y mínimos encontrados por interpolación, encontramos un máximo en el sample 12,574 con valor 0,99667 y un mínimo en el sample 17,603 con valor -0.99680. Como ves, uno se acerca al sample 13 y el otro al 18. Los máximos y mínimos se acercan más al máximo y mínimo teóricos de una senoidal que es 1 y -1.

Bueno, pues sobre el valor eficaz de la señal ya hablaré más adelante y te iré explicando qué es y cómo se calcula (es que se me ha hecho ya tarde). Por lo pronto te enseño una función que he hecho otra vez en matlab para calcular el valor eficaz de las funciones interpoladas de una secuencia de samples:


```
function salida = RMSsplines(samples)

  Y_1=samples(length(samples));

  intparcial=0;
  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;

	intparcial=intparcial+A^2/7+(A*B+2*B*D+C^2)/3+(2*A*C+B^2)/5+(A*D+B*C)/2+C*D+D^2;
	Y_1=Y0;
	end

salida=sqrt(intparcial/length(samples));

 return
end
```

Esta función es una función que devuelve un valor que es el valor RMS REAL, no aproximado, de la secuencia de polinomios interpoladores. Cuanto más parecidos sean los polinomios interpoladores a la onda real, más se parecerá el valor RMS real de los polinomios al valor RMS de la onda que aproximan.

Saludos.


----------



## asherar (Jul 30, 2013)

seaarg dijo:


> Intente aplicar lo que me especifico Eduardo sobre la funcion sinc
> 
> 
> 
> ...



Los detalles de la fórmula anterior se explican en el  Teorema del Muestreo.


----------



## seaarg (Jul 30, 2013)

palurdo, como decimos por aca, sos un grande! 

Una explicacion digna de analizarla con paciencia e ir aplicandola. Por lo pronto, les cuento las ultimas modificaciones:

Pude resolver lo que pregunte sobre las bases de tiempo haciendome un dibujito en un papel, y llegando a la misma conclusion que vos me decias aqui: Que si tengo un tiempo entre samples de 25ns, para poder tener 1us/div debo intercalar 4 samples, con un pequeño truco:

- Selecciono en el mux una frecuencia de sampleo para 5us/div y utilizo 5 samples, 2 reales, 3 calculados. (Esto claro, x100) Entonces, entre los 2 samples reales tengo 50ns / 5 = 10ns * 100 = 1000ns/div

Es decir, la escala minima "real" de mi osciloscopio es de 2.5us/div pero no utilizo esa para 1us/div sino la siguiente 5us/div y uso interpolacion.

A partir de alli para abajo, si uso la 2.5us/div como frecuencia de muestreo y hago:
- 2500ns / 5 (2 reales, 3 calculados) = 500ns/div
- 2500ns / 10 = 250ns/div
- 2500ns / 25 = 100ns/div
- 2500ns / 50 = 50ns/div

y ahi paro nomas, porque en esa ultima ya tengo tan solo 2 samples reales por division. Podria seguir pero ya es suficiente por el momento.

No me puse a investigar el tema de los calculos de datos utiles (muchas gracias! me diste mucha informacion para aplicar) porque me encontre con un problema:

Me fabrique un oscilador de onda cuadrada con un cristal de 16mhz, inversor schmitt trigger 74HC14 y un 74HC4040 para hacer divisiones. Lo alimente con una bateria que tenia por ahi de 6v.

Resulta que, midiendo con osciloscopio comercial, veo que la señal tiene entre 0v-5.86v y yo, como un tonto, tengo una referencia del adc de 5v entonces veo la parte de abajo de la señal pero arriba "hace tope" con vcc.

Entonces, aproveche el divisor 2:1 que tenia hecho con los mux y lo deje por defecto activado, reemplazando la cuenta 5 / 256 en el programa por 10 / 256 para saber a cuantos volts corresponde cada incremento en el resultado del ADC.

Entonces, incremente el voltaje de los operacionales desde +-5v a +-7v (son rail to rail) y puse un regulador 78L05 con un preset de ajuste a la referencia del ADC.

Con esto logre que, si tengo una señal de 10vpp en la entrada del osciloscopio (5 positivos, 5 negativos), en la entrada del ADC tengo 5vpp centrados en 2.5v (gracias al ex divisor 2:1, ahora llamado 1:1)

Lo que antes era 1:1 ahora pasa a llamarse 1:2 (multiplicador ahora) quedandome:
1:1
1x2
2.5:1 (ex divisor 5:1)
1x10 (ex multiplicador 1x5)

De esta forma, logre que la señal en pantalla, en las mismas escalas, coincida con la de un osciloscopio digital comercial (prestado jeje) que estoy usando para calibrar esta cosa. Su señal de calibracion de 1khz 2vpp se ve exactamente igual en el mio que en el comercial 

Ahora la mala noticia: Lo unico que pude ver bien fue una señal cuadrada de 500khz. Una de 1mhz se ve pero ya no tan cuadrada (en el comercial se ve perfecta) y mas alla van cada vez peor hasta que una de 4mhz ya no se ve (quedan solo picos abajo).

Creo que es porque sin querer arme un pasabajos en la entrada principal del osciloscopio: Antes del primer operacional tengo una resistencia en serie con el op-amp de 10K en paralelo con un capacitor de 22n y luego de esto, los diodos 1n4148 de proteccion hacia VCC y VEE y ahi si va al + del operacional.

Asi que, antes de proseguir con el soft voy a reparar esta parte del hard. Creo que los operacionales que uso, segun su datasheet y a +-7v deberian ir sobrados para una frecuencia de 4mhz, y por otro lado, el ADC corriendo a 40mhz deberia darme 10 puntos claros bien cuadrados (aclaro: su ancho de banda analogico es de 75mhz) Estoy casi seguro que cometi algun error de diseño y formo un pasabajos accidental por algun lado. Aclaro que en esta prueba no use interpolacion.

De todos modos, estoy conforme con el resultado pero tengo mucho hardware para tan poco asi que lo voy a exprimir al maximo 

Por otro lado, apenitas resuelva esto me pongo con las mejoras de software que me propones, son excelentes!

asherar: Muchas gracias por el link! Habia leido varias veces tu post sobre la placa de captura a 32mhz y obtuve datos utiles de el para el diseño de esto. Lamentablemente sitio donde tenias detalles hace rato que esta caido.

Tambien podria seguir aumentando la capacidad de esto haciendo muestreo equivalente. De hecho lo hice con un trigger inexacto (aplicando 1 NOP de 83.33~ previo al inicio de captura en la segunda pasada) y luego intercalando los samples de 1 y otra pasada, pero la verdad que el resultado que obtuve con interpolacion 1 matematica fue practicamente el mismo y el hardware se complica bastante si tenemos que generar retrasos variables exactos (de fraccion de la frecuencia de muestreo).


----------



## asherar (Jul 30, 2013)

Yo también estuve probando las funciones de matlab que ha ido subiendo *palurdo* y andan de 10. 
Es realmente notable lo que se logra con un algoritmo tan sencillo. 
Aunque los métodos de interpolación parecen muy enredados si se los ve en forma analítica (como el de la función sinc), pero pueden ser llevados a la práctica con algunos truquitos ingeniosos. 

En tu caso, *seaarg, *trabajando con los datos una vez que están en la PC, no se tiene problemas con la velocidad. Pero si se hicieran a bordo del pic, además de simple, el algoritmo tiene que ser rápido. 

Con respecto a mi sitio, ya veré de subir todo otra vez. 
Mi intención no era interpolar sino hacer muestreo de tiempo equivalente, y para eso necesitaba un retardo 
como el que mencionas al final: ajustable, estable y preciso, con una duración de fracción del retardo (1/32MHz ~ 30 ns), y eso no sé si existe. 
Varias veces intenté continuar el proyecto pero algunos temas de hardware me "ganaron" y ya no tengo el tiempo libre que tenía antes. Todo se puede resolver, pero lleva tiempo. 

PD: Acá hay algo interesante sobre calcular funciones trigonométricas con numeros enteros 
(pero todavía me falta entender cómo es el uso práctico). 
En forma de presentación.
En forma de texto plano.


----------



## asherar (Jul 31, 2013)

Estuve pensando en este tema y me surgieron algunas dudas. 

Por lo que se ve en las rutinas que subió palurdo a la cantidad de muestras originales N, la rutina le devuelve un vector de N*M. Usar M = N la verdad es que me parece un poco mucho. 
Tal vez con cuadruplicar (M=4) la cantidad de muestras originales alcanzaría para ver una curva suave. Si nó la curva mostrada puede ser demasiado irreal. 
Digo esto porque en el intervalo ente dos muestras uno puede intercalar un punto, o dos y ya con eso ganar suavidad en la curva. Pero por la misma cuestión que existe la frecuencia de Nyquist, no tiene sentido agregar más puntos. Si se agregan demasiados puntos se está asumiendo que lo que hay entre medio de dos muestras ES suave. 
Además se carga el sistema con cálculos que en realidad resultan innecesarios. 
Entiendo que usando un algoritmo sencillo no queda otra, pero tampoco es bueno exagerar con la información que uno agrega a dedo. 
Tal vez, la salida elegante sería intercalar más o menos puntos, de modo tal que la curva completa nunca pase de un cierto número, digamos N*M=1000 o algo así. 
Pero reitero: que se agregen siempre una cantidad fija (M=N) de puntos ente cada 2 de un conjunto de N no me parece algo conveniente. 

Otra cosa. 

Cuando te falla la última parte de la curva porque faltan las muestras, es porque te faltan las muestras del "futuro", necesarias para aplicar el algoritmo splines2. 
En ese caso (es una sugerencia), y sólo para ese tramo final, se podría probar con la aproximación que tiene en cuenta solamente muestras del "pasado" (la primera versión de splines). Es una idea. No sé qué pueda pasar. 

Bueno, no sigo (no molesto) más.

....

Bueno sí, pero espero que no sea molestia. 
Estuve haciendo una prueba rápida con la función sinc de Matlab y uno de los ejemplos que viene en el Help. 
Les dejo el texto de la ayuda, donde explica que la función sinc es la anti-transformada de Fourier del pulso rectangular que vale 1 en (-pi,pi) y 0 en (-inf,-pi) y (pi,inf).

Copia de pantalla de la ayuda:


Ahora un código modificado por mí a partir de un ejemplo.

```
%
% Reconstruccion de una señal a partir de su muestreo
% usando la funcion sinc de Matlab
%
clear all
clc

N=20;
t = (1:N)';             % Un vector columna de coordenadas de tiempo
s = - t .* t .* cos(t); % Un vector columna de muestras (coseno)
                        % El producto .* opera elemento a elemento
figure(1);
xlabel('tiempo');
ylabel('amplitud');
hold on                 % Para que no se borre 
plot(t,s,'o');          % La grafica de las muestras
texto=sprintf('N=%d',N);text(0,100,texto);

N1 = -2;                % Limite izquierdo de la grafica
N2 = N+2;               % Limite derecho de la grafica
M = 100;                % Numero de puntos a calcular
texto=sprintf('M=%d',M);text(0,-100,texto);

ts = linspace(N1,N2,M)'; % Coordenadas de tiempo interpoladas

% Operacion usando la funcion sinc de Matlab
y = sinc(ts(:,ones(size(t))) - t(:,ones(size(ts)))')*s;

figure(1)
plot(ts,y,'b-')         % La grafica de la funcion reconstruida
```
Esto da: 


Lo interesante es que la función sinc que viene implementada usa la misma cantidad de puntos que el vector de muestras dato, o mayor como en el ejemplo. 

M es la cantidad de muestras del vector de interpolación que incluye los puntos de la zona inicial y final que serían "puntos desperdiciados". 

Si se hace: 

ts = linspace(0,N,M)'; % Coordenadas de tiempo interpoladas

la función "linspace" no genera puntos desperdiciados. 



Acá, además de la curva en rojo se muestran los puntos calculados como x azules. 
Es para ver que aunque se agregaron solo 3 puntos por cada uno de muestra la curva se ve suave.

Para dejar este tema redondito quedaría hacer este mismo cálculo sin  usar la función sinc de Matlab, es decir recorrer el arreglo elemento  por elemento.


----------



## seaarg (Jul 31, 2013)

asherar:

Ante todo, no es molestia! Todo suma al conocimiento y es bienvenido al menos por mi.

En el planteo de la cantidad de muestras inventadas estoy de acuerdo: Estariamos asumiendo que la señal es limpia y la realidad casi siempre tiene espurias por lo tanto, es conveniente no pasarse en cantidad de datos calculados.

En mi caso particular, fabrico 1us/div, 500ns/div, 250ns/div, 100ns/div y 50ns/div. Sin embargo, no me olvido cuando hago pruebas de que hasta 500ns/div tengo una idea mas o menos acertada de que lo que veo "debe" ser como lo real. De ahi para abajo ya no confio tanto jeje, estan ahi, mas que para que sean reales, para poder "hacer zoom" en señales mas rapidas y que sean mas comodas de ver en pantalla. Por otro lado, en mi implementacion hay un checkbox en el soft que decide si aplica interpolacion lineal o splines por lo tanto ante cualquier sospecha se deja lineal.

En lo de la falla en la ultima parte puede que tengas razon, cuestion de probar que el algoritmo anterior se use en ese caso pero creo recordar (no mucho) que hacia lo mismo. Probare y les cuento cuando resuelva el problema que tengo del "pasabajos" en algun lado.

Al respecto de lo que pusiste sobre sinc(x), originalmente antes que palurdo me mostrara splines buscaba informacion precisamente sobre como aplicar esto, ya que el lenguaje de programacion que estoy usando que es vb.net no tiene dicha funcion (y ningun lenguaje de proposito general creo) entonces quise comprender la definicion de la funcion para poder hacerla yo en codigo (como vos decis, recorrer elemento x elemento) y ahi fue cuando el cerebro me hizo division por cero 

Adelanto que voy a proseguir con esto, no esta abandonado, solamente que tengo que encontrar porque en mi diseño frecuencias cuadradas de 4mhz se ven como triangulo y ya vi que asi estan llegando al ADC. Por lo pronto ya vi que solo la punta de osciloscopio que tengo (dice que es de 20mhz) de por si ya le quita un poco de amplitud a la señal original, lo cual me parece raro pero para nada imposible. Despues, luego del primer operacional de entrada vuelvo a tener la señal como la lee un osciloscopio comercial pero antes de el esta muy reducida en amplitud (algo raro pasa y tengo que enchufar el soldador para encontrar el culpable, me juego que es el conjunto R paralelo a C que esta en serie con la entrada del primer operacional)

Para peor, los cambios que hice para corregir la poca amplitud de entrada disponible funcionan de 10 en cuanto a su objetivo, pero me afectaron el funcionamiento del trigger (ya no dispara bien) y tengo que ver porque razon 

Hoy no pude dedicarme por el trabajo pero se vienen nuevos avances supongo. Por lo pronto tambien intente empezar a hacer esquematicos al menos en proteus para compartir con ustedes, solo que voy poquito a poquito.


----------



## Sebastian1989 (Ago 1, 2013)

Me alejo un poco del tema de este tópico pero ¿Has simulado la red RC de entrada con algún programa?.
Podrías subir el esquema de la entrada de tu osciloscopio para ver si surgen ideas para alguna mejora.


----------



## palurdo (Ago 1, 2013)

asherar dijo:


> Estuve pensando en este tema y me surgieron algunas dudas.
> 
> Por lo que se ve en las rutinas que subió palurdo a la cantidad de muestras originales N, la rutina le devuelve un vector de N*M. Usar M = N la verdad es que me parece un poco mucho.
> Tal vez con cuadruplicar (M=4) la cantidad de muestras originales alcanzaría para ver una curva suave. Si nó la curva mostrada puede ser demasiado irreal.
> ...



Gracias por tus apuntes. La verdad es que depende del fin que le des a la interpolación tendrá sentido o no que pongas muchas o pocas muestras. Como ves el vector interpolado es M veces N, donde M el el número de intervalos (puntos+1) a interpolar entre dos puntos reales. Por simplicidad en la función no incluí en el vector de salida el último punto, que sería la última muestra real del vector N, pero no afecta a la interpolación. El M que se utilice ya lo elijes dependiendo de los puntos que necesitas para representar en pantalla y las muestras reales representables. Es decir, que si por ejemplo tienes N muestras representables en pantalla en tu escala de tiempo y tienes disponibles 1000 posiciones en pantalla, entonces para esa escala utilizas 1000/N=M puntos a interpolar entre muestras de forma que M*N=1000. Dejar fijo M sería una pérdida de tiempo para calcular puntos que al final no vas a representar en pantalla porque no caben en los pixeles. Pongamos que tienes 100 puntos en 1us, pues para representar 1000 puntos en pantalla necesitas partir los tramos entre puntos reales en 10 segmentos interpolados (9 puntos intermedios) para tener los 1000 pixeles. Si ahora haces zoom x2, para 500ns, para representar en pantalla sólo tienes la mitad de puntos, por lo que para 1000 pixeles necesitas aumentar al doble el número de puntos interpolados. Hacer una matriz cuadrada M=N es una barbaridad, ya que si tienes por ejemplo 500 puntos, estarías interpolando 500 segmentos entre cada punto, por lo que tendrías 250000 muestras interpoladas, algo completamente innecesario, ya que con M=2, es decir, intercalando un punto interpolado entre 2 reales, ya consigues los 1000 pixeles.

Por lo mismo dejar un M demasiado bajo no es buena idea ya que si tienes pocas muestras, y pocos puntos interpolados, al unir con rectas sigues teniendo una curva irreal, pero con discontinuidades inventadas (esquinas en las uniones de las rectas). Lo ideal es tener un M tal que N*M=Total pixeles a dibujar.

Por otro lado, creo que dices que las curvas interpoladas con muchos puntos quedan demasiado sintéticas. A pesar de que estoy de acuerdo contigo, debo de comentar que discrepo en algunos matices pero propondría alguna solución a esto aunque sería empeorar la calidad de la interpolación. Discrepo porque el que use el osciloscopio debe de tener claro el ancho de banda de éste, de forma que frecuencias superiores al límite no pueden representarse. Por ejemplo, si tienes un ancho de banda de entrada 10MHz, y un muestreador a 200MHz, siempre tendrás curvas suaves de como mínimo de 10 muestras, aunque la señal original tenga transiciones repentinas entre esas 10 muestras, ya que esas transiciones tendrán frecuencia superior al ancho de banda de entrada del osciloscopio. No tiene sentido pensar en un osciloscopio con un ancho de banda de 10MHz y un muestreador a 200MHz, porque desperdicias capacidad de muestreo ya que el ancho de banda es muy pequeño en conmparación con la serie muestreada. Resulta que *interplolando, lo que haces al aumentar la cantidad de samples es aumentar la frecuencia de muestreo* ya que para el mismo periodo temporal tienes mayor cantidad de samples. *Pero esto no aumenta el ancho de banda *de la señal de entrada. Por lo tanto, en lugar de tener una entrada de 10MHz y un muestreador a 200MHz, se puede tener una entrada de 10MHz y un muestreador de 20MHz (o 25MHz para no complicarse la vida con el filtro antialiasing), y a partir de los samples a 20MHz interpolar para convertir la señal muestreada a 200MHz con ancho de banda 10MHz lo cual es el mismo resultado que haber muestreado directamente a 200MHz. En todo caso se puede recordar en pantalla que el límite superior a reprentar en la gráfica, ya sea interpolada o no, es el límite superior del ancho de banda de entrada al capturador.

El aspecto de la 'calidez' de las señales reales VS las interpoladas viene a raiz de que en pantalla aparecen componentes frecuenciales de banda ancha que no pertenecen a la señal. Dos motivos ocurren para esto:

-Ruido de captura.
-Ruido de cuantización.

El ruido de cuantización aparece en la señal real al limitar las muestras a 8 bits (LSB a -49dB). Con muestras de 24 bits (LSB a -144dB) esto no es un problema porque el valor mínimo es muy pequeño y la curva muestreada es prácticamente matemática. Cargad con Audacity cualquier archivo de audio directamente extraido a 16 bits de un CD de música, id haciendo zoom hasta ver la onda definida en pantalla, y vereis lo que os digo cuando los samples a 16bits parecen samples calculados matemáticos y sin embargo son samples reales de audio.

El ruido de fondo de captura (noise floor), es una porción de ruido blanco de banda ancha que aparece en las etapas amplificadoras antes de la captura, y ese suelo de ruido cuando es mayor que el ruido de cuantización es cuando aparece representado en la señal (por debajo del ruido de cuantización no tiene efecto ya que el que domina es el de cuantización). El ruido blanco aparece igualmente espaciado en todas las frecuencias, por lo que para saber su magnitud sólo se necesita conocer su valor eficaz RMS.

Reproducir el ruido de cuantización en los puntos interpolados no tiene mayor secreto que truncar (o redondear para eliminar un pequeño sesgo) la parte decimal de los valores interpolados. De hecho al dibujar en pantalla los puntos interpolados ya se está cuantizando porque no existen pixeles fracionales sino enteros.

Reproducir el ruido de fondo captura, a costa de degradar la señal, sería medir el valor RMS del ruido de captura cuando no hay señal, y aplicar un ruido blanco de ese valor a las muestras interpoladas. Este proceso, que empeora la calidad de la señal resultante (añade ruido) pero dependiendo de la aplicación es apropiado para enmascarar ciertos aspectos de una conversión de los datos para hacerlos más naturales, se le conoce como *DITHERING*.



asherar dijo:


> Otra cosa.
> 
> Cuando te falla la última parte de la curva porque faltan las muestras, es porque te faltan las muestras del "futuro", necesarias para aplicar el algoritmo splines2.
> En ese caso (es una sugerencia), y sólo para ese tramo final, se podría probar con la aproximación que tiene en cuenta solamente muestras del "pasado" (la primera versión de splines). Es una idea. No sé qué pueda pasar.



No habría problema, ya que cada conjunto de 3 puntos en splines y 4 puntos en splines2 es independiente del resto, así que es como si hasta 0...N-1 calcularas la gráfica en splines2 y de N-1 a N la gráfica para splines 1. Ese tramo seguiría siendo suave ya que la pendiente de la curva anterior a N-1 tiene la misma pendiente que la curva posterior. 

De todas formas yo para splines2 asumo que la señal es periódica, es decir que la primera muestra de la señal sigue después de la última. Creo que seaarg la muestra del futuro la hace 0, o igual es por el padding que hace en su señal para rellenar los puntos que le faltan en la  pantalla.

Intenté hacer una función de interpolación splines3, donde consideraba las derivadas segundas de los puntos Y[-1] e Y[2], para completar un polinomio de grado 5 con curvas más suaves (dos niveles de derivadas). Los cálculos de los coeficientes A,B,C,D,E,F tuve que hacerlos con Derive ya que resolver el sistema de ecuaciones a mano se hizo una pesadilla. El resultado, como era de esperar, fue una peor interpolación que usando splines de tercer grado. Si quieres comparar, te dejo la función que hice, pero no la uses para interpolaciones reales porque aproxima peor que splines2.


```
function salida = splines3(samples, npuntos)

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

  for i=1:(L-1)

	Y0=samples(rem(i-1+L,L)+1);
	Y1=samples(rem(i-1+L+1,L)+1);
	Y2=samples(rem(i-1+L+2,L)+1);
	Y_1=samples(rem(i-1+L-1,L)+1);
	Y_2=samples(rem(i-1+L-2,L)+1);
	Y3=samples(rem(i-1+L+3,L)+1);

	A1= -(Y_2 -9*Y_1 -7*Y1 +15*Y0)/2;
	B1= 3*(Y_2 -8*Y_1 -6*Y1 +13*Y0)/2;
	C1= -(3*Y_2 -19*Y_1 -13*Y1 +29*Y0)/2;
	D1=(Y_2+Y0-2*Y_1)/2;
	E1=Y0-Y_1;
	F1=Y0;

	A2= (Y3 - 9*Y2 + 15*Y1 - 7*Y0)/2;
	B2= -(2*Y3 - 21*Y2 + 36*Y1 - 17*Y0)/2;
	C2= (Y3 - 13*Y2 + 23*Y1 - 11*Y0)/2;
	D2= (Y2 - 2*Y1 + Y0)/2;
	E2= Y1 - Y0;
	F2= Y0;
	
	A=(A1+A2)/2;
	B=(B1+B2)/2;
	C=(C1+C2)/2;
	D=(D1+D2)/2;
	E=(E1+E2)/2;
	F=(F1+F2)/2;
	
	salida(i,:)=A*X.^5+B*X.^4+C*X.^3+D*X.^2+E*X+F;

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

 return
end
```




asherar dijo:


> Bueno sí, pero espero que no sea molestia.
> Estuve haciendo una prueba rápida con la función sinc de Matlab y uno de los ejemplos que viene en el Help.
> Les dejo el texto de la ayuda, donde explica que la función sinc es la anti-transformada de Fourier del pulso rectangular que vale 1 en (-pi,pi) y 0 en (-inf,-pi) y (pi,inf).
> 
> ...



El problema de sinc es que por cada muestra lo que hace es el sumatorio de toda la serie de puntos para sacar la funcion de onda donde después se remuestrea en la serie interpolada, ten en cuenta que en cálculo sinc(T)*s estás haciendo un producto matricial, lo que viene a ser una convolución de los valores del vector T con el vector S. Esto hace que el algoritmo sea muy lento de implementar. Consigues el mismo resultado implementando la interpolación por FFT que es mucho más eficiente. En matlab/octave existe la implementación de la función de interpolado FFT, llamada intepft. La onda de tu ejemplo usando interpft:



Por cierto, haz la prueba, pero t=linspace(N1,N2,M) es exáctamente lo mismo que hacer T=(N1N2-N1)/(M-1):N2), ya lo verás... (mas que nada para entender la implementación de la función linspace).

De todas formas, la interpolación por splines genera un resultado también bastante decente sin complicar tanto los calculos, aquí abajo la interpolación de dicha señal usando splines2:



Gráfica generada con las mismas condiciones que las interpolaciones anteriores, lo único es que no se generan los últimos 5 puntos, por lo que el vector interpolado tiene 95 en lugar de 100 puntos.

seaarg:

¿Has probado a atenuar la onda cuadrada de 4MHz al mínimo que te permita capturar aunque sea 1 bit?, más que nada para descartar un problema de Slew-Rate de los operacionales para señales rápidas de niveles grandes. Si pudieras poner unas capturas de lo que te sale con la onda "mala" te lo agradeceríamos. Un saludo y buen trabajo.


----------



## palurdo (Ago 1, 2013)

*asherar:*

Quería comentarte que cuando me refería que la cantidad de puntos a interpolar dependía de la aplicación, es porque yo inicialmente implementé este sistema de interpolado para otros menesteres. 

Hace como cosa de 15 años, me encontré con un problema. Quería pasar por un filtrado de barrido un clip de audio para recrear el efecto de "barrido de filtro resonante" muy popular en muchas muestras que se utiliza en música electrónica. Lo traté de hacer con el filtrado FFT del programa cooledit (ahora adobe audition). El filtro FFT del cooledit es muy bueno en filtrado estático, pero en filtrado de barrido es horrible ya que el "morphing" que hace del filtro a medida que pasa el tiempo no es suave sino que va a tramos, por lo que no se escucha el efecto de barrido sino un efecto como escalonado. 

Así que programe un programa que lo que hacía era tener un espacio gráfico donde había una linea horizontal que acababa en 2 vértices, y el usuario podía mover esos vértices entre los niveles de filtrado iniciales y finales, pero también añadir vértices interiores para construir una transición a tramos. Utilizaba la interpolación para generar la transición lo más suave posible. Por ejemplo, si el usuario usaba 5 puntos, se creaban 4 polinomios para crear la transición del filtro (se representaba en pantalla la curva interpolada por pixeles para que el usuario tuviera una idea de como quedaba la función de transición). La interpolación del coeficiente del filtro se hacía cada sample del clip de audio. Pongamos que había 1 segundo de sample a 44100 smp/s, pues con cada sample de entrada se recalculaban los coeficientes del filtro a partir del valor de tiempo que me daba la interpolación dentro del segmento de 2 puntos correspondientes de la gráfica. Al final, había calculado los coeficientes del filtro 44100 veces, de ellas 44095 eran posiciones temporales interpoladas, y 5 reales pertenecientes a la gráfica de transición. La verdad es que el programa funcionaba bastante rápido funcionando en un Pentium I@100MHz y el efecto era el pretendido, tener un filtro de barrido de dos polos movibles con una transición muy suave, al igual que lo hacía el REBIRTH-303 con su resonador. 

Cuando quise convertir el código en un plug-in VST, el Norton Antivirus me hizo un desastre en la partición del disco duro y perdí todos los proyectos (no tenía copia de seguridad, en aquella época era jóven y atrevido...)

Como ves, en esta aplicación fue bastante bien interpolar 44095 puntos desde 5 puntos reales, pero claro, poco tiene que ver aquella implementación con esta del osciloscopio, sólo que ambos usan splines.

Saludos.


----------



## seaarg (Ago 1, 2013)

Sebastian1989 dijo:


> Me alejo un poco del tema de este tópico pero ¿Has simulado la red RC de entrada con algún programa?.
> Podrías subir el esquema de la entrada de tu osciloscopio para ver si surgen ideas para alguna mejora.



No hay problema. Si la simule en proteus con un grafico de frecuencia pero no tengo el operacional exacto que estoy usando asi que no es confiable.

No tengo esquematico de la entrada (si una simulacion pero toda desprolija) sin embargo, la misma es casi una copia de esto: http://www.bitscope.com/design/hardware/pdf/Bs11-4.PDF

Solamente que el recuadro verde "vertical input buffers" + el primer mux son reemplazados por un operacional en donde tengo la sospecha que arme un pasabajos sin querer.

Tambien cambie con respecto a este esquematico el ultimo operacional que va al ADC, pero no es el problema (por ahora) porque a partir del 1 operacional de entrada la señal ya esta distorsionada a frecuencias altas.

Esta tarde agarro el soldador y elimino la R limitadora de entrada a ver que pasa



Palurdo, la verdad que no, no probe atenuar aun porque con el osciloscopio real medi que a la salida de mi primer op-amp ya veo distorsionado y eso me confundio pero bien puede ser como sugieren, que el slew rate sea insuficiente (130v/us!)

Esta tarde retomo las pruebas y les envio gustoso unas capturas de pantalla tanto del osciloscopio mio como de un comercial para la misma señal.


----------



## seaarg (Ago 1, 2013)

Bueno, les cuento que pude arreglar un poco la etapa de entrada analogica y si, diversas resistencias que habia puesto en simulacion andaban bien pero en la realidad me hacian de pasabajos. Lastima que ahora que arregle esto, se me jodio el trigger, anda bien en determinadas condiciones. Mañana tocara revisar esa parte y luego a seguir con las matematicas!

Aca una cuadrada de 1mhz 5.44vpp medida con un osciloscopio hantek, tomada desde la entrada de mi osciloscopio (para contar la punta del mismo que esta distorsionando un poco por alguna razon)



Y aca, mi osciloscopio midiendo la misma señal



Ahi la escala es 1us/div y estoy sampleandola a 20mhz (en vez de 40mhz) e interpolando 5 puntos en forma lineal. No adjunto la interpolacion con calculo splines debido a que empece a tener un error de desborde de vectores que tengo que corregir. De todos modos por lo que pude ver splines anda muy bien para cuadradas tambien.

Notese que en la parte de arriba esta bien plano, debido a que la señal se esta pasando de la escala vertical max y actuan los diodos de recorte. Pero noten la parte inferior como es muy parecida al osciloscopio comercial.

Tengo cada vez mas ganas de terminar este bicho de una vez para embarcarme a hacer uno portatil, mucho mas sencillo, que solo use el ADC del pic que es bien lento y un buffercito de entrada. Con la implementacion de splines de numeros enteros que desarrollo palurdo se podria obtener algo bastante interesante para ser chiquito, simple y portatil (con un GLCD como tantos osciloscopios que andan dando vueltas por ahi)

Ah! Adjunto tambien una simulacion en proteus de la entrada (1 solo canal + trigger), les pido disculpas por lo horriblemente chancha que esta su distribucion. Aclaro que no estoy usando en la realidad el mismo operacional por lo tanto, en el aparato, el preset de entrada esta hacia VEE en vez de VCC y el valor no es el mismo sino mas grande.


----------



## asherar (Ago 2, 2013)

seaarg dijo:


> Tengo cada vez mas ganas de terminar este bicho de una vez para embarcarme a hacer uno portatil, mucho mas sencillo, que solo use el ADC del pic que es bien lento y un buffercito de entrada. Con la implementacion de splines de numeros enteros que desarrollo palurdo se podria obtener algo bastante interesante para ser chiquito, simple y portatil (con un GLCD como tantos osciloscopios que andan dando vueltas por ahi)



Ojo que los GLCD son en general bastante lentos. 
He trabajado con algunos de 128x128 y 128x64 pixels y tienen cierto tiempo de refresco sólo 
tolerable por los que alguna vez sufrimos las XT 286, como yo. 

Hace un tiempo estoy poniendole mis horas sueltas a un aparatito construido con tu misma filosofia. 
No había pensado trabajarle el interpolado porque no hay mucho lugar dónde poner más puntos. 
Sacando algunas columnas para textos me quedan unos escasos 100 pixels de ancho para la gráfica. 
A ver si les gusta la pinta. 
La placa de adaptación es externa para poder probar cosas sin tener que desarmar todo. 
En este caso toma señal analógica de dos "electret".


----------



## seaarg (Ago 3, 2013)

Yo empece con una TI99 asi que estoy acostumbrado a la lentitud jaja. Otra idea un poco mas complicada pero aun factible es usar una pantalla de notebook vieja. Venden por ebay una plaquita que es una interfaz vga-pantalla lcd con el elevador de tension y todo. Si con un micro extra (tipo controlador de pantalla) se genera la señal vga se podria hacer un osciloscopio digital de mesa. Las posibilidades estan ahi pero primero voy a terminar este que ya esta casi para ir a un gabinete.


----------



## palurdo (Ago 3, 2013)

las pantallas de 7" de las tablets casi son mas baratas que los glcd. Otra cosa es con que controlarlas... Hay que pensarlo muy bien a ver si va a salir mas caro el collar que el perro jejejeje


----------



## seaarg (Ago 3, 2013)

palurdo dijo:


> ```
> function [xmax ymax xmin ymin] = maxminsplines(samples)
> 
> Y_1=samples(length(samples));
> ...



Palurdo, he implementado en mi programa la funcion para extraer datos de los samples crudos que me detallas mas arriba. Los resultados parecen ir teniendo forma.

Tuve que establecer el valor inicial de xmin = 9999 ymin = 256 (mi maximo es 255) para que esos datos tomen algun valor sino nunca entra a los IF. Por lo que veo en tu codigo los inicializas a cero, entonces supongo que tus "samples" incluyen numeros negativos (en mi caso, de 0-255)

Los valores que obtengo para YMax YMin se ven muy razonables y varian cuando cambio la amplitud de la señal, me falta traducirlos a volts segun la escala que este la entrada.

Sin embargo, para Xmax y Xmin obtengo datos un poco erraticos, esta fue mi implementacion para calcular frecuencia (pongo unidades aqui en la formula asi se entiende que es cada parte)

sampling_rate = [1=25ns,2=50ns,3=100ns]

dDiff = ((Abs(xmax - xmin) * 2 [dobleperiodo] ) * (25ns * sampling_rate)) / 1000000000 [cuantos ns es 1 segundo]

dDiff = 1 / dDiff

La frecuencia de la señal de prueba es otra vez 20khz y obtengo, cuando 1 solo periodo de la señal ocupa toda la pantalla (todos los samples) valores similares a 20.000hz pero que varian bastante rapido entre 19 y algo y 21 y algo...

Sin embargo, en escalas mayores donde tengo varios picos de la onda seno x pantalla, la frecuencia calculada (por ende, los valores xmax y xmin) pegan saltos muy grandes y muy rapidos.

Una aclaracion importante seria que en tu funcion existe esta linea:

if X1>0 && X1<1 && imag(X1)==0

Que yo la tuve que cambiar (para probar nomas) por esto:

if X1>0 && X1<1

Debido a que no se que hace la funcion imag(x) (sera numeros imaginarios o algo asi?)

Si pudieses aclararme que hace, ya que no esta definida en el codigo y supongo que es alguna funcion de octave/mathlab, quiza solo eso sea el problema que estoy teniendo, ya que sin esa ultima condicion, en casi toda vuelta del ciclo for entro al IF.

Muchas gracias!


----------



## asherar (Ago 3, 2013)

En Matlab imag(x) devuelve la parte imaginaria del número x. 
Lo necesita para asegurarse que la x (que sale de la raíz de una ecuación cuadrática) sea un número real.

Lo raro es que en varios lugares usa: xmin = i + X2, y similares. 
Los complejos no están ordenados, entonces no tiene sentido decir que un nro. complejo sea mínimo o máximo de un conjunto.


----------



## palurdo (Ago 4, 2013)

asherar dijo:


> En Matlab imag(x) devuelve la parte imaginaria del número x.
> Lo necesita para asegurarse que la x (que sale de la raíz de una ecuación cuadrática) sea un número real.
> 
> Lo raro es que en varios lugares usa: xmin = i + X2, y similares.
> Los complejos no están ordenados, entonces no tiene sentido decir que un nro. complejo sea mínimo o máximo de un conjunto.



Mea culpa. Ya sabes que en ingenieria usamos j en lugar de i para variable compleja asi que no meimporto usar i como varianble de indice. Igual habiendo puesto xmin=ind+X2 se habria visto mas claro. De todas formas me aseguro que todos los valores sean reales.

Seaarg: 

En tu caso, simplemente asegurate de que lo que hay dentro de la raiz cuadrada sea mayor que 0. En caso de queno sea asi no calcules X1, X2, Y1, Y2 porque esos valores no caen dentro de el eje real.



la funcion con inicializar xmin, xmax ymin e ymax con x e y del primer sample debe ser suficiente.


----------



## seaarg (Ago 4, 2013)

Bueno, en vb.net tenemos un constructor de clase numerica tipo "complex" al cual asigne (luego de lo que me explico asherar) el resultado de la operacion (defini variables X1 y X2 como complejos) y en todas los calculos donde intervienen esas variables, uso la propiedad real, que devuelve la parte real de un numero complejo. En el IF comparo la propiedad "imaginary" contra cero. Creo que logre hacer lo mismo que tenes hecho en tu codigo con estas herramientas.

Sin embargo el problema persiste y creo que encontre la causa aunque no sabria como solucionarlo:

Probe implementar tanto tu funcion, asi como tambien una funcion mia bien "bruta" donde simplemente promedio 3 valores continuos de Y en el tiempo, buscando maximo y minimo a lo largo de todos los samples. Y el resultado de ambas es similar en comportamiento (no asi en valores, estoy 99% seguro que los valores devueltos por tu funcion son mas acertados)

Resulta que, si tengo varios picos de la onda seno en pantalla (o sea, en el set de samples), la funcion va recorriendolos, buscando el maximo y minimo valor de Y (luego de los calculos claro) y sucede lo siguiente:

Para un set de 896 datos: 0-895 en eje x

- En una primera pagina, digamos que encuentra el minimo Y en el sample 80 y el maximo en el 160
- En una segunda pagina, encuentra el minimo en el mismo sample 80 pero el maximo por el 500 (debido a que el "maximo" del punto 160 no fue mayor a un pico que hubo en esta pagina, en el sample nro 500)
- En una tercera pagina, el minimo lo encuentra en 200 y el maximo en 896

Por supuesto los valores que pongo son de ejemplo extremo de lo que sucede. La cuestion es que claro, como la señal no es algo perfecto y una cresta puede tener 1 punto mas que otra tanto para arriba como para abajo, detecta el periodo entre 1 cresta y otra, o entre 2 crestas, o 3 y asi, haciendo que el calculo de distancia entre crestas que hago sea variante entre cada "pagina" de samples.

Intente hacer una especie de atajo para evitar esto, reemplazando:


```
if Y1 >ymax 
  ymax=Y1;
  xmax=i+X1;
end
```

por:


```
if Y1 >ymax + 10
  ymax=Y1;
  xmax=i+X1;
end
```

para que asi, solo cuando el valor es bastante mas grande cambie la variable ymax (seria como una deteccion de pico maximo y minimo digamos) con iguales resultados.

Tambien probe implementar un detector para que solo cuando la señal "sube" compare YActual con Ymax y solo cuando "baja" Yactual con Ymin con los mismos resultados (claro! me di cuenta despues de probar que es lo mismo)

Los valores de YMAX e YMIN se mantienen mas o menos constantes, con variaciones de 1 o 2 enteros nomas, por lo que el calculo de vmin y vmax de la señal ya estoy en condiciones de hacerlo pero el calculo de la frecuencia se me resiste por esa variabilidad.

Estoy pensando como puedo hacer para detectar 1 cresta y su valle siguiente, o viceversa, para capturar solamente 1 periodo (siempre duplicando la distancia en X claro porque la F es entre cresta y cresta)

Esto se pone cada vez mas lindo! 

Edito con una pregunta:

¿Que opinan de esto? Si en la funcion de palurdo pongo mas codigo que haga algo como:

- Detectar todas las crestas, es decir, para cada vez que y > ymax almaceno ese valor como ymax pero, una vez que el valor de Y sea menor a ymax (con cierto margen, o sea, cuando empieza a bajar) "reseteo" ymax para detectar la siguiente cresta... y asi sucesivamente hasta terminar los samples.

Luego, teniendo las varias distancias en X entre las distintas crestas, hacer un promedio entre todas ellas y que eso sea la diferencia xmax - xmin.

Igualmente, la deteccion de que termino una cresta es mas facil decirla que programarla.


----------



## palurdo (Ago 4, 2013)

Claro, es que la función detecta el máximo y mínimo global y lo que necesitas es detectar 2 máximos locales o dos mínimos locales. Se puede modificar la función para que además de los globales te detecte todos los locales de la señal. Como no sabemos cuantos máximo y mínimos locales puede tener la señal en la gráfica, tendrás que ir guardándolos en un vector de máximos y un vector de mínimos. 

Bueno, pues cada vez que se entra en uno de los if X1>0 y X1<1 entonces ese X1 (o X2) es o bien un máximo o un mínimo. Para saber si es un máximo o un mínimo haces el cálculo de 6*A*X1+2*B y si ese resultado es negativo tienes un máximo, y si es positivo tienes un minimo, y si es 0 no es ni un máximo ni un mínimo sino un punto de inflexión (que a ti no te es útil). Además X_ es un máximo local si X[i-1] y X[i+1] son menores, y es un mínimo local si son mayores. 

Por lo tanto por cada dos puntos consecutivos yo comprobaría si cada uno de los puntos X, X1 y X2 son máximos o mínimos y en caso de que algunos lo sean, añadirlos a un vector. 

Por definición, los puntos inicio y final no pueden ser máximos locales, por lo que no te fijes si X[0] y X[895] con mínimos y máximos locales.


Por otro lado este método sólo sirve para calcular la frecuencia de una senoidal pura o de una onda monótona (cuyas crestas y valles coinciden con su periodo), pero falla cuando la señal sea la suma de varias senoidales, por ejemplo, porque calcularás la frecuencia de un armónico de la señal. Lo que puedes hacer es calcular amplitud y frecuenci entre un máximo y un mínimo local, despues generar una senoidal con esa amplitud y frecuencia (y ajustando la fase para que coincida el máximo y el mínimo de la senoidal generada con los máximos y mínimos obtenidos) y le restas esa senoidal a tu onda original. repites el proceso, y si la amplitud de la señal resultante es menor que la de tu senoidal, pues ya tienes la componente principal. En caso contrario vuelves a calcular máximos y mínimos y vuelves a generar otra senoidal y la restas, y así hasta que la senoidal que generes sea mayor que la amplitud de la onda que te queda. En ese caso de la lista de senoidales generadas, la que tenga mayor amplitud es la componente principal de frecuencia._


----------



## asherar (Ago 4, 2013)

seaarg dijo:


> ... Otra idea un poco mas complicada pero aun factible es usar una pantalla de notebook vieja. Venden por ebay una plaquita que es una interfaz vga-pantalla lcd con el elevador de tension y todo. Si con un micro extra (tipo controlador de pantalla) se genera la señal vga se podria hacer un osciloscopio digital de mesa. .


Estuve gugleando y encontré formas simples de acoplar salidas VGA con entradas a LCD ( ENLACE al artículo 1) y también con RGB+Sync tipo salida de PlayStation (ENLACE al artículo 2).
Aparentemente todo pasa por tener el pinout del LCD. Si no ... 





​ 
Yo tengo una Sharp LM130SS1T611, que venía en una Compaq Presario, hoy fenecida de muerte natural. 
Si alguien tiene el pinout por favor le pido que lo suba por acá (y si no es demasiada molestia me avise por MP). 
Mil gracias.


----------



## seaarg (Ago 5, 2013)

Buscando programas que implementen funciones e interpolaciones, he caido en esto:
http://jonw0224.tripod.com/ppmscope.html

Se lo ve prometedor al chiquitin. El soft es bastante completo y tiene interpolaciones splines y sinc

Viene con soft para PC, el firmware del pic (en assembler!) PCB, lista de componentes, etc.etc.

Editado: Pongo este post en el thread de asherar sobre osciloscopios, ya que me parece mas adecuando que en este. Perdon!


----------



## seaarg (Ago 8, 2013)

Les cuento, al final pude obtener la frecuencia de la señal con un metodo medio feo que es un detector de cruce por cero. Cuento los samples entre 2 cruces por cero, eso lo multiplico por 2 y con el valor de tiempo x sample obtengo el periodo y de ahi la frecuencia.

No es un buen metodo, si es simple.

Adapte desde codigo C hacia vb.net unos 3 o 4 codigos que pude encontrar para calcular FFT, y al final me quede con el mas prolijo y completo, que es el que adjunto (el archivo se llama func.c) y que pertenece al soft del osciloscopio PPMScope

Creo haberlo adaptado bien, reemplazando los tipos de datos, las copias de vectores, los operadores especificos de C, etc.

Para el que lo quiera ver, la cosa esta en la funcion sigFreq() donde se aplica una ventana de Hann a los samples, luego se realiza la FFT y luego se busca la frecuencia fundamental pasandole las magnitudes calculadas en la FFT a traves de fundFreqBin, que por lo que entendi, calcula la frecuencia dividiendo la FFT en contenedores que representan cada frecuencia posible.

Dicho esto: NO HAY FORMA! algo esta mal en mi codigo o en mi interpretacion de los resultados, porque SIEMPRE, sea la que sea la frecuencia que estoy sampleando, obtengo que la mayor magnitud esta en la posicion 256 de la FFT, que es la mitad de los datos enviados a la misma. Si uso la mitad de los datos de la FFT (n/2) entonces la posicion siempre es 128 y asi sucesivamente.

Lo mismo me paso con otros algoritmos para hacer FFT que fui traduciendo desde C a net.

¿Cual puede ser mi concepto errado por el cual esto no anda? Es decir: Puede que este interpretando mal la salida de la FFT en mi cabeza o habra un error de codigo?

Mi interpretacion es: La salida de la FFT es el dominio de las frecuencias. El valor de vector de mayor amplitud indica la posicion de la frecuencia fundamental. Luego: Como se cual es dicha frecuencia? (relacion con sampling rate 100% seguro pero no se como aplicarla) Esto esta definido en fundFreqBin() pero no usa sampling rate.

Aclaro que no pido que me resuelvan el programa, sino que me ayuden a entender el resultado de la FFT.

Adjunto tambien un archivo FFTAnalisis.zip que es un excel donde puse, en la columna A los samples originales (señal 20khz sampleados a 5.000.000 hz) y en la columna B un dump de lo que devuelve la funcion faFourier, en el vector de magnitudes.


----------



## asherar (Ago 8, 2013)

Esta te devuelve el valor, escalado en frecuencia, del índice del máximo de la FFT. 


> /* Returns the fundamental frequency of the signal
> * arr[]: array in which to find the fundamental frequency
> * len: length of the array
> */
> ...


Y usa ésta otra para buscar el índice. 


> /* A Helper function. (Una función de ayuda)
> * Find the fundamental frequency bin, assumes no peak merging
> (Encuentra el índice de la frecuencia fundamental, supone que no hay picos mezclados)
> */
> ...


Esta te devuelve el índice de la muestra de mayor amplitud. 
Si el error es acá, debe ser que te queda siempre una muestra de valor muy alto al final, y por eso te queda ahí al salir. 

Otra cosa más: 
En la planilla hay un dato, en el casillero C1: (4,09E+75) que no me parece normal. 
Hay otros similares por ahí, y por eso la gráfica te da las abscisas entre 0 y 10^80. 
La forma de la FFT tiene muchos picos de ese orden de magnitud.
Si es la Transformada de la sinusoide de la primera gráfica, debería ser un solo pico en la f del seno. 

No probaste buscar algo ya hecho para PIC en la página de Microchip ?
Buscando aquí por ejemplo  De paso ya nos queda para el portátil !!! 

En el foro se ha tratado el tema: https://www.forosdeelectronica.com/f24/fft-pic-13856/

Implementation of Fast Fourier Transforms





Reading and Using Fast Fourier Transforms (*FFT*)


----------



## seaarg (Ago 9, 2013)

Bien, primero que nada voy a revisar si la traduccion de C a vb.net la hice bien (creo que si) de esta funcion


```
int fundFreqBin(double magn[], int len)
{
    double avg, a;
    int i;

    avg = sigMax(&magn[1], len-1)/2;

    i = 0;
    a = magn[i];
    while(magn[++i] < a)
        a = magn[i];
    a = magn[i];
    while(magn[++i] > a || a < avg)
        a = magn[i];
    i--;

    if(i >= len - 1)
         return 0;
    return i;
}
```

A esto en vb.net


```
Private Function fundFreqBin(ByVal magn() As Double, ByVal len As Integer) As Integer

        Dim avg, a As Double
        Dim i As Integer
        Dim tmagn(magn.Length - 2) As Double
        Array.ConstrainedCopy(magn, 1, tmagn, 0, len - 1)
        avg = sigMax(tmagn, len - 1) / 2

        i = 0
        a = magn(i)
        i = 1   '++i del while
        While magn(i) < a   ' La primera vez = 1
            a = magn(i)
            i = i + 1
        End While
        a = magn(i)
        i = i + 1
        While magn(i) > a Or a < avg
            a = magn(i)
            i = i + 1
        End While
        i = i - 1
        If i >= len - 1 Then
            Return 0
        Else
            Return i
        End If

    End Function
```

Luego, tengo un problemon y es el siguiente:

BIN = Fs/N = 40000000/1024 = 39062Hz.

Eso significa que, para 1024 muestras en la maxima frecuencia de sampleo del osciloscopio, la variacion entre 1 bin y otro es ENORME!

Se me ocurre que podria tomar, por ejemplo 1 muestra de cada 2 para asi "bajar" la F de sample a la mitad y completar el resto con ceros o algo por el estilo pero aun asi los saltos entre BIN son grandes.

Tambien sino, fabricar muchas mas muestras usando la interpolacion splines de palurdo para achicar de esta forma el bin size. El problema de esta forma seria que tendria que procesar una FFT para, digamos, 1024*10 muestras, con todo lo que eso significa en consumo de memoria y en velocidad de proceso (ya hasta la PC se pondria lenta creo)

Es decir: para hacer FFT me encuentro con no solo los problemas tecnicos de programacion, etc, sino de que la Fs del osciloscopio es demasiado grande en la escala de tiempos mas chica. ¿Es correcta esta interpretacion?


----------



## asherar (Ago 9, 2013)

A mí la traducción de las rutinas me parece correcta. 



seaarg dijo:


> Luego, tengo un problemon y es el siguiente:
> 
> BIN = Fs/N = 40000000/1024 = 39062Hz.
> 
> ...



No entiendo bien tu problema. 
Si no querés tener un BIN demasiado grande bastaría con expresar tus frecuencias en kHz o en MHz. La cantidad representa lo mismo pero con un número más bajo. 
Si usas los splines en realidad te aumenta el valor en lugar de bajar. 
Creo que por ese tema algunas rutinas hacen cuentas con el índice y luego calcula la frecuencia en las unidades que corresponde. 
Pero no sé para qué usas luego esos números, así que tal vez mi observación no tenga sentido.


----------



## seaarg (Ago 9, 2013)

Tal vez estoy interpretando mal los calculos que se hacen con la FFT, tengo entendido lo siguiente:

- Tengo un vector de 1024 samples en dominio de tiempo
- Paso esos samples por FFT
- Tengo ese vector, pero en dominio de frecuencia
- Aplico calculo de magnitud para cada elemento del vector
- Obtengo el indice del vector donde se encuentre la mayor magnitud, salteando el x[0] que es DC
- Calculo el bin size con la Fs de esta forma: bin = Fs / 1024
- Multiplico el indice de mayor magnitud por este bin size para obtener la F fundamental de los samples.

Entonces, para Fs = 40000000 en esa cantidad de samples el bin es = 39062,5

Supongamos que la mayor magnitud de la FFT esta en el indice 1, entonces el resultado del calculo de frecuencia fundamental va a dar 39062,5 HZ (y de ahi para arriba, asi que algo esta mal!! no podria obtener una F menor a 40khz, para esta cantidad de samples, a este sample rate tan grande)

Llego a esa conclusion porque en diversas pruebas, con diversos set de samples, estoy obteniendo esos indicios, pero estoy casi seguro de que hay algo mal o en las presuposiciones, o en las rutinas que hice.

Consegui un programita en vb 6.0 que hace el grafico de la FFT de una senoidal generada por el mismo (con funcion sin()) y el mismo devuelve en que indice esta la max amplitud...

... entonces, le reemplace la generacion artificial de la senoidal por un conjunto de mis samples:

Resulta que este si me da el indice siempre en los primeros lugares (no como mi rutina que siempre me daba en 256 o similares) y, como sample esas pruebas a 20us/div, hice la siguiente cuenta:

bin = 5000000 / 1024 = 4882,81

La aplique y parece ir funcionando, me da valores en HZ similares a la senoidal que yo genero con la placa de sonido. Por supuesto, la resolucion es de 4.8khz entre un indice y el siguiente, por lo tanto, la resolucion es pesima pero al menos es logica.

Tal vez aqui tenga que releer algo que explico palurdo mucho mas arriba, sobre la frecuencia fundamental que cae entre dos indices enteros.

Es curioso tambien como, en este programita esta la opcion de usar 256, 512, 1024, 2048 y 4096 samples. En 1024 (yo tengo 896 muestras y completo con numeros 127 que seria el GND) anda bastante bien, pero a medida que aumento la frecuencia de la senoidal sampleada, tengo que subir la cantidad de samples. Sin embargo, si uso la opcion de 4096 samples para todas las pruebas, en la mayoria de las situaciones la magnitud maxima de la FFT se encuentra en el indice 0 (como si fuera DC)

En fin, me voy a volver loco jeje. Lo que voy a hacer ahora es estudiar otras implementaciones de la FFT, en php por ejemplo asi puedo ver si estoy cometiendo un error y donde.

Me encantaria poder interpretar los simbolos matematicos de la descripcion de la FFT. Solo conozco el de sumatoria  entonces no puedo hacer un algoritmo en codigo directamente a partir de la formula. (ay, la ignorancia!!!)

Edito: Me encontre un error en mi traspaso de codigo desde C a vb.net y ahora en vez de tener la mayor amplitud en la mitad de la FFT siempre, la tengo siempre en el indice de la FFT = 1 (descartando el 0 que es DC)

Tambien adapte ese programita vb 6.0 que comente que funciona, y le puse que trabaje con 8192 samples (896 reales entre -127 y 127 y padding a 0). Desde frecuencias en el generador de algunos cientos de hz hasta los 20khz me da una medida que parece ser buena, para una Fs de 5.000.000

Obviamente, este programa que encontre una un algoritmo un poco distinto para hacer la FFT. Podria adaptarlo pero el codigo C que tenia me parecia mas elegante (este codigo de vb 6 funciona pero es bastante chancho a mi gusto)

Quisiera aclarar que, para el calculo de frecuencia deje de intentar aplicar max y min locales como me decia palurdo, y me emperre en sacar adelante FFT, por la razon de que luego, usando la misma FFT puedo sacar otros valores como THD, etc. Y aclarar que la funcion que me brindo palurdo quedo super firme para vmax y vmin. Esos estan saliendo perfectos.

Sigo intentando buscarle la vuelta. Les comentare resultados!


----------



## seaarg (Ago 10, 2013)

Vuelvo con mas novedades. Pude implementar la FFT 

Consegui un programita hecho en vb.net que hace la FFT, tiene aplicacion de ventana de hamming, hanning y blackman (o ninguna) y funcionaba bastante bien, asi que agarre el codigo fuente y lo adapte a mi programa.

Estoy usando la barbaridad de 131072 muestras (896 de verdad, el resto cero) y es demasiado rapido, las procesa sin problemas!

Se eligio semejante numero para poder disponer de un bin size de 40000000 / 131072 = 305hz

Lo proximo seria agregarle un pequeño codigo que cambie la cantidad de samples que se le envian a la FFT segun el sample rate del ADC, porque es inutil mandar tantos cuando sampleo mas lento que 40mhz. De esa forma, seria mas rapido aun.

El algoritmo es bien feo en cuanto a codigo de programacion (uso de goto)  pero lo importante es que funciona muy bien. Lo adjunto aqui por si alguien lo quiere investigar, yo solamente le tuve que adaptarlo a mi señal, sampling rate variable, etc y hacer el calculo de magnitud sqrt(re_^2+im^2) para poder obtener la frecuencia fundamental.

Luego de unas cuantas adaptaciones, esto es lo que logre visualizar:



Ahi se ve la señal de los 2 canales, la frecuencia y periodo (aun por adaptar a algo legible) de los mismos, los niveles de señal obtenidos por splines, la ventana de grafico de FFT para el canal A (verde) mostrando el espectro, y deje tambien ahi la ventanita del programa generador de señal, para comparar lo que dice que esta generando con lo que arrojan de resultado las matematicas.

En fin, parece que todo va bien! ahora ya me puse ambicioso y le quiero sacar mas informacion a la FFT, como THD, ruido, etc.

Si alguien quiere el programa completo del osciloscopio, digame y con gusto lo subo aqui. Con explicaciones del codigo mediante, otra persona puede hacerse su placa de captura y utilizar este soft (que aun le falta mucho! esta verde, pero pinta bien)_


----------



## Chico3001 (Ago 14, 2013)

y en que quedo el proyecto?? ya funciono???


----------



## Sebastian1989 (Ago 14, 2013)

Excelente proyecto, me están dando ganas de armarme uno, voy a buscar por eBay algunos ADCs rápidos.

PD: Con solo un "poco" mas de trabajo lo podrías vender por medio millón de dolares, solo necesitas ampliar el ancho de banda a 62GHz y tomar 160GS/s. Vídeo de alguien "jugando" con un osciloscopio de medio millón de dolares


----------



## seaarg (Ago 14, 2013)

Chico3001 dijo:


> y en que quedo el proyecto?? ya funciono???



Si que funciono! luego de aprender mucho, renegar mucho y recibir mucha ayuda, ahora tengo una herramienta bastante util para proximos diseños 

Aca va una imagen de una cuadrada generada por un cristal de 8mhz y dividida con un contador. La señal es de 2mhz y esta siendo sampleada a 40mhz, utilizando la interpolacion SPLINES de palurdo (no me canso de agradecerte! jeje) Tambien se ve la ventanita de FFT donde veo la fundamental y las distintas armonicas. Dicha ventana la deje chica para que entre en esta captura de pantalla pero como tengo 2 monitores, normalmente arrastro esa ventana al 2do monitor y veo la FFT en full screen con todo su detalle.



En esta otra, vemos 2 senoidales en canal A y B, y ademas lo que se ve como "fantasma" por detras son 2 senoidales de la referencia. Le adicione al programa poder guardar una captura para despues cargarla de fondo como referencia visual (cuando comparo señales, por ejemplo, para ver si un filtro pasa bajos filtra mas o menos segun valores de componentes)



Y por ultimo, pongo una foto del bicho en su gabinete temporal-definitivo jeje, son tapas de lectora de CD recicladas, ahora tengo que hacerle el frente nomas 



Mejore varios aspectos del programa y ahora anda mucho mas rapido tambien!

Caracteristicas generales:

2 Canales + trigger externo
sample rate max: 40mhz
Escalas de tiempo: Desde 50ns/DIV a 1 S/DIV
Memoria: 896 bytes x2

Analogica:
Entrada de señal 1:1 = 10vpp (100vpp con el selector de la punta) ITEM A MEJORAR
Ancho de banda: 75-130mhz (no probado)
Ancho de banda probado: 4mhz (con señal cuadrada)
Divisores: 1:1, 2:1
Multiplicadores: 1x2, 1x10
AC/DC

Trigger:
Externo, Canal A, Canal B. Nivel seleccionable con mouse (5vpp a mejorar) Trigger automatico, single shot, slope + y -

Software:
Interpolacion SPLINES y Lineal para bases de tiempo 1us/DIV hacia abajo. Para arriba es captura real
Guarda y carga a memoria un DUMP de la señal
Guarda y carga a memoria una referencia de señal para mostrar de background
Grafico FFT del power spectrum de la señal
Mediciones (por el momento): Vmax, Vmin, VPP, Periodo, Frecuencia fundamental (a futuro, THD, etc)
Cursor con mouse que informa VPP, Periodo y Frecuencia del area seleccionada en pantalla.

En fin, estoy muy contento. Lo que andaba haciendo ahora es ver la forma de pasar esto a esquematicos, no conozco mucho de programas para esto y hacerlo en livewire seria un infierno, sumado a que no existen todos los integrados ni en ese ni en proteus. Empece a hacer una especie de esquematico con photoshop, pegando las imagenes de los integrados en los datasheets, para unir con lineas pero ahi esta, pendiente.

El PCB esta adjuntado mas arriba aunque ya ha tenido cambios de componentes y algunas cirugias.

Adjunto aqui tambien el soft, tanto de la PC como el del PIC por si lo quieren chusmear o sacar ideas, o usarlo de base para algun hardware propio.

Esta hecho en vb.NET con framework 4. La comunicacion usb es via HID. Es un soft un tanto complejo y no esta muy comentado para entenderlo facil, pero cualquier duda, aqui estoy para lo que necesiten.


----------



## cosmefulanito04 (Ago 14, 2013)

Muy bueno como va quedando, el tema de medir bien una señal rectangular, está en el flanco (si es muy corto el time rise), tratá de probar con senoidales para ver hasta donde podés medir bien un armónico puro y ver como se comporta el vertical en todo el ancho de banda.


----------



## seaarg (Ago 15, 2013)

Exacto, hasta ahora genero las señales senoidales con la compu y no da mas de 20khz. Ahora que lo tengo en gabinete lo puedo llevar a lo de un tecnico amigo que tiene un generador de señales  para probarle un poco mas, pero creo que su generador no da mas de 200khz. Voy a tener que fabricarme un generador de senoidal por arriba del mhz

De todos modos, aclaro que la señal cuadrada de 2mhz se ve muy parecida en un osciloscopio comercial (con un poco mas de detalle nomas) pero la forma sale asi de fea desde el generador


----------



## cosmefulanito04 (Ago 15, 2013)

seaarg dijo:


> Exacto, hasta ahora genero las señales senoidales con la compu y no da mas de 20khz. Ahora que lo tengo en gabinete lo puedo llevar a lo de un tecnico amigo que tiene un generador de señales  para probarle un poco mas, pero creo que su generador no da mas de 200khz. Voy a tener que fabricarme un generador de senoidal por arriba del mhz
> 
> De todos modos, aclaro que la señal cuadrada de 2mhz se ve muy parecida en un osciloscopio comercial (con un poco mas de detalle nomas) pero la forma sale asi de fea desde el generador



Podrías probar con una cuadrada lo más zarpada que puedas (hasta 20MHz sería lo máximo que podés medir) y quedarte solo con la fundamental usando un buen filtro, hacelo de varios ordenes pasivo, no importa, ya que tenés un OCR para medir y lo podés usar de patrón para comprobar la amplitud de la senoidal.

Cuadrada --> Capacitor (para eliminar DC) --> Filtro pasabjos 1er orden -> Filtro pasabjos 1er orden .... etc

Lo malo de hacer eso es que atenúas bastante la salida, pero en este caso te resulta indistinto. Eso si, trata de usar componentes que tengan el mejor comportamiento posible a esas frecuencias.

¿Qué contadores usás para hacer el divisor de frecuencia?

*Editado:*

Vos usas un cristal, de ahí podés sacar la senoidal.


----------



## seaarg (Ago 15, 2013)

Ah muy buena idea. El generador de cuadrada que uso es un cristal + inversor. Pongo la sonda en la entrada del inversor y tengo la senoidal del cristal.

El contador es 74HC4040

Hoy lo lleve de un tecnico amigo y su generador de cuadrada, triangular, rampa y senoidal dio como maximo a 200khz. El osciloscopio respondio perfecto a todas ellas 

Eso si, como su computadora daba 4.2v en el usb, tuve que poner la entrada del osci en AC porque sino no andaba el trigger (error mio de diseño del trigger, solo anda cuando los integrados estan alimentados con al menos 4.8v como es en mi PC)


----------



## cosmefulanito04 (Ago 15, 2013)

seaarg dijo:


> ...
> El contador es 74HC4040
> ...



En base solo a ese contador, como mencioné en este _tema_ (y trajo cierta controversia ), el problema con esa rectangular/cuadrada es que el Tr es muy corto (7nS típicos trabajando a 5v) para tu ancho de banda y posiblemente también para el ancho de banda del OCR patrón que estás usando (por eso se ve "feo"), ya que al menos necesitarías un ancho de banda mayor a 100MHz.

Si vas a trabajar con señales rectangulares, asegurate de tener Tr de 50 nS o más para medirlo bien.


----------

