Estimación de retardo de tiempo de las señales del osciloscopio utilizando correlación cruzada

12

Grabé 2 señales de un oscopio. Se ven así: ingrese la descripción de la imagen aquí

Quiero medir el retraso de tiempo entre ellos en Matlab. Cada señal tiene 2000 muestras con una frecuencia de muestreo de 2001000.5.

Los datos están en un archivo csv. Esto es lo que tengo hasta ahora.
Borré los datos de tiempo del archivo csv para que solo los niveles de voltaje estén en el archivo csv.

x1 = csvread('C://scope1.csv');
x2 = csvread('C://scope2.csv');  
cc = xcorr(x1,x2);
plot(cc);  

Esto da este resultado: ingrese la descripción de la imagen aquí

Por lo que he leído, necesito tomar la correlación cruzada de estas señales y esto debería darme un pico en relación con el retraso de tiempo. Sin embargo, cuando tomo la correlación cruzada de estas señales obtengo un pico en 2000 que sé que no es correcto. ¿Qué debo hacer con estas señales antes de cruzarlas? Solo busco alguna dirección.

EDITAR: después de eliminar el desplazamiento DC, este es el resultado que ahora obtengo:
ingrese la descripción de la imagen aquí

¿Hay alguna manera de limpiar esto para obtener un retraso de tiempo más definido?

EDITAR 2: Aquí están los archivos:
http://dl.dropbox.com/u/10147354/scope1col.csv
http://dl.dropbox.com/u/10147354/scope2col.csv

Nick Sinas
fuente
¿Cómo, exactamente, estás haciendo la correlación cruzada? En respuesta a su pregunta directa, no debería necesitar hacer nada a sus señales antes de la correlación cruzada, aunque en algunos casos el filtrado primero ayuda a eliminar el ruido que puede distorsionar los resultados.
Jim Clay
1
Publique el código que ha utilizado y, lo que es más importante, una gráfica de la señal de correlación cruzada. Algunas herramientas / bibliotecas colocan la puntuación (lag = 0) en el centro del gráfico; No recuerdo si Matlab hace eso.
pichenettes
@pichenettes: publicación actualizada
Nick Sinas
@JimClay: publicación actualizada
Nick Sinas
@NickS. Si sus señales están perfectamente alineadas, obtendrá un pico en el medio de su diagrama de cc. Entonces, el pico en 2000 significa que no hay demora. Ahora digamos que tiene un retraso de 10 muestras, lo que significa que la señal2 está a 10 muestras de la señal1. Esto solo moverá su pico en el cc de 2000 a 2010 (o 1990). Entonces su retraso de tiempo corresponde a su ubicación pico real, MENOS 2000.
Spacey

Respuestas:

11

@NickS

Dado que no es seguro que la segunda señal en los gráficos sea de hecho una versión únicamente retrasada de la primera, deben intentarse otros métodos además de la correlación cruzada clásica. Esto se debe a que la correlación cruzada (CC) es simplemente un estimador de máxima verosimilitud si sus señales son versiones retrasadas entre sí. En este caso, claramente no lo son, por no decir nada sobre la no estacionariedad de ellos tampoco.

En este caso, creo que lo que puede funcionar es una estimación del tiempo de la energía significativa de las señales. De acuerdo, 'significativo' puede o no puede ser algo subjetivo, pero creo que al observar sus señales desde un punto de vista estadístico, podremos cuantificar 'significativo' e ir desde allí.

Con este fin, hice lo siguiente:

PASO 1: Calcule los sobres de señal:

Este paso es simple, ya que se calcula el valor absoluto de salida de la Transformada de Hilbert de cada una de sus señales. Existen otros métodos para calcular sobres, pero esto es bastante sencillo. Este método esencialmente calcula la forma analítica de su señal, en otras palabras, la representación fasorial. Cuando tomas el valor absoluto, estás destruyendo la fase y solo después de la energía.

Además, dado que buscamos una estimación de retraso de tiempo de la energía de sus señales, este enfoque está garantizado.

ingrese la descripción de la imagen aquí

PASO 2: Elimine el ruido con filtros mediales no lineales que conservan los bordes:

Este es un paso importante. El objetivo aquí es suavizar sus envolturas de energía, pero sin destruir ni suavizar sus bordes y tiempos de subida rápidos. En realidad, hay un campo completo dedicado a esto, pero para nuestros propósitos aquí, simplemente podemos usar un filtro Medial no lineal fácil de implementar . (Filtrado medio). Esta es una técnica poderosa porque, a diferencia del filtrado medio , el filtrado medial no anulará sus bordes, pero al mismo tiempo 'suavizará' su señal sin una degradación significativa de los bordes importantes, ya que en ningún momento se realiza ninguna aritmética en su señal (siempre que la longitud de la ventana sea impar). Para nuestro caso aquí, seleccioné un filtro medial de tamaño de ventana de 25 muestras:

ingrese la descripción de la imagen aquí

PASO 3: Eliminar tiempo: Construir funciones de estimación de densidad del núcleo gaussiano:

¿Qué ocurriría si miraras la trama anterior de lado en lugar de la forma normal? Matemáticamente hablando, eso significa, ¿qué obtendrías si proyectaras cada muestra de nuestras señales sin ruido en el eje de amplitud y? Al hacerlo, lograremos eliminar el tiempo, por así decirlo, y podremos estudiar únicamente las estadísticas de la señal.

Intuitivamente, ¿qué sale de la figura de arriba? Si bien la energía del ruido es baja, tiene la ventaja de que es más "popular". En contraste, mientras que la envolvente de señal que tiene energía es más enérgica que el ruido, está fragmentada a través de umbrales. ¿Qué pasa si consideramos la 'popularidad' como una medida de energía? Esto es lo que haremos con la implementación (mi crudo) de una función de densidad del núcleo (KDE), con un núcleo gaussiano.

Para hacer esto, se toma cada muestra y se construye una función gaussiana utilizando su valor como media, y se selecciona un ancho de banda preestablecido (varianza) a priori. Establecer la varianza de su gaussiano es un parámetro importante, pero puede establecerlo según las estadísticas de ruido según su aplicación y las señales típicas. (Solo tengo tus 2 archivos para activar). Si luego construimos la Estimación de KDE, obtenemos la siguiente gráfica:

ingrese la descripción de la imagen aquí

Puede pensar en el KDE como una forma continua de un histograma, por decirlo así, y la varianza como su ancho de bin. Sin embargo, tiene la ventaja de garantizar un PDF sin problemas en el que luego podemos realizar el primer y segundo cálculo derivado. Ahora que tenemos los KDE gaussianos, podemos ver dónde las muestras de ruido alcanzan su mayor popularidad. Recuerde que el eje x aquí representa las proyecciones de nuestros datos en el espacio de amplitud. Por lo tanto, podemos ver en qué umbrales es más "enérgico" el ruido, y esos nos dicen qué umbrales debemos evitar.

En la segunda gráfica, se toma la primera derivada de los KDE gaussianos , y seleccionamos la abscisa de la primera muestra después de la primera derivada después del pico de la mezcla de gaussianos para alcanzar un cierto valor cercano a cero. (O primer cruce por cero). Podemos usar este método y estar 'seguros' porque nuestro KDE se construyó con gaussianos suaves de ancho de banda razonable, y se tomó la primera derivada de esta función suave y sin ruido. (Por lo general, las primeras derivadas pueden ser problemáticas en cualquier cosa que no sean señales SNR altas ya que aumentan el ruido).

Las líneas negras muestran entonces en qué umbrales sería prudente 'segmentar' la imagen, para evitar todo el ruido de fondo. Si luego aplicamos a nuestras señales originales, alcanzamos los siguientes gráficos, con las líneas negras que indican el inicio de la energía de nuestras señales:

ingrese la descripción de la imagen aquí

Esto produce así un muestras.δt=241

Espero que esto haya ayudado.

Spacey
fuente
Guau. Muchas gracias. Todas estas son técnicas nuevas para mí que comenzaré a investigar. ¿Hay alguna forma de que pueda echar un vistazo al código matlab que usaste?
Nick Sinas
Así que tengo los pasos 1 y 2 realizados en Matlab, y mis resultados coinciden con los suyos, pero tengo problemas con el paso 3. ¿Qué funciones usaste?
Nick Sinas
@NickS. Pregunte ... y recibirá, envíeme un correo electrónico y puedo enviarle el código que utilicé.
Spacey
@Mohammed ¿Podría publicar su código para estimar el retraso de tiempo? Le he enviado un correo electrónico con respecto a este asunto, así que por favor ayuda
6

Hay algunos problemas al hacer esto con autocorrelación

  1. Enorme desplazamiento DC (ya arreglado)
  2. Ventana de tiempo: el xcorr () de Matlab tiene la convención molesta de esencialmente "poner a cero" la señal en ambos extremos a medida que desliza el retraso de tiempo. Es decir, la ventana de datos es una función del desfase temporal. Esto creará una forma triangular para cualquier señal estacionaria (incluidas las ondas sinusoidales). Las mejores opciones son elegir una ventana de correlación para que el tamaño de la ventana más el lapso de tiempo máximo se ajusten a su ventana de datos totales, o para normalizar la correlación cruzada por el número de muestras no rellenadas.
  3. Las dos señales no me parecen particularmente correlacionadas. La forma es algo similar, pero el espaciado específico de picos y caídas es bastante diferente, por lo que dudo que incluso una autocorrelación adecuada produzca mucha información aquí.

Un enfoque mucho más simple sería usar un detector de umbral para encontrar los puntos de partida y simplemente usar la diferencia entre estos puntos como el retraso.

Hilmar
fuente
4

Como indicaron las pichenettes, en este caso un pico en el medio de la salida indica un retraso de 0. El desplazamiento del pico desde el punto medio es su retraso de tiempo.

EDITAR: Me preocupa que la correlación sea casi un triángulo perfecto. Eso me indica que la correlación cruzada no está haciendo normalización de potencia. Eso da un sesgo injusto a los retrasos más pequeños sobre los retrasos más grandes. Modificaría su llamada de xcorr a "cc = xcorr (x1, x2, 'imparcial');".

Eso no es una solución perfecta, porque los resultados de gran retraso ahora son más inestables que los resultados de bajo retraso porque se basan en menos datos. Un pico grande en las extremidades puede ser falso por la misma razón que puede obtener el 100% de cara y sin colas con solo unos pocos lanzamientos de monedas, mientras que es extremadamente improbable que ocurra en muchos lanzamientos.

Jim Clay
fuente
¿Significa que las señales no se retrasan?
Nick Sinas el
No estoy seguro, ¿dónde está el pico? Puedo ver que está cerca del medio, pero no está claro que en realidad esté en el medio. Además, hay un problema de normalización de energía que abordaré en una edición de mi respuesta.
Jim Clay
El parámetro 'imparcial' definitivamente lo hace ver mejor. más de lo que esperaría Seguiré investigando esto. Gracias.
Nick Sinas
@JimClay Quizás Nick S, está correlacionando los sobres de sus señales y no las señales reales, (¿Nick es esto cierto?). Esto produciría (aproximadamente) esta forma triangular que imagino.
Spacey
2
@NickS. El comentario de Mohammad me hizo ver y darme cuenta de que tienes una gran compensación de CC que está arruinando tus resultados. Reste la media de ambas señales y luego ejecute xcorr sobre ellas. Lo intentaría sin la opción "imparcial" primero.
Jim Clay
4

Como lo han señalado los demás, y parece que te has dado cuenta en base a tu última edición de la pregunta, no parece que la correlación cruzada te dé una buena estimación del retraso de tiempo para los conjuntos de datos mostrados. La correlación mide la similitud en la forma entre dos series de tiempo al deslizar una sobre la otra por un rango de retrasos de tiempo y calcular un producto interno entre las dos series en cada retraso. El resultado tendrá una gran magnitud cuando las dos series sean cualitativamente similares, o estén "correlacionadas" entre sí. Esto es similar a cómo un producto interno de dos vectores es más grande cuando los dos vectores apuntan en la misma dirección.

El problema con los datos que mostró es que (al menos para los fragmentos que podemos ver) no parece haber mucha similitud en la forma. No hay demora que pueda aplicar a una de las señales para que se parezca a la otra, que es exactamente lo que está haciendo al calcular su correlación cruzada.

Sin embargo, hay casos en los que la correlación cruzada es útil. Supongamos que su segunda señal fue realmente una versión del original en función del tiempo, incluso con algún ruido adicional agregado:

a = csvread('scope1col.csv');
a = a - mean(a);               % to remove DC offset
b = a(200:end) + sqrt(0.05)*randn(1801,1);
figure; subplot(211); plot(a); subplot(212); plot(b)

ingrese la descripción de la imagen aquí

Ahora no está claro de inmediato que las dos señales estén relacionadas por un retraso de tiempo. Sin embargo, si tomamos la correlación cruzada, obtenemos:

[c,lags] = xcorr(a,b);
igure; plot(lags,c); grid on; xlabel('Lag'); ylabel('Cross-correlation');

ingrese la descripción de la imagen aquí

que muestra un pico en el retraso correcto de 200 muestras. La correlación puede ser una herramienta útil para determinar el retraso de tiempo, cuando se aplica a conjuntos de datos que contienen el tipo correcto de similitud.

Jason R
fuente
¿Alguna idea de qué más podría hacer? ¿Quizás otra técnica que no sea la correlación cruzada o quizás algún tipo de filtro? Gracias.
Nick Sinas
@NickS. También lo he mirado y no son copias retrasadas el uno del otro. Dicho esto, ¿quieres estimar el retraso de la energía ? Creo que eso tendría más sentido en este caso, VS retraso de las señales . Si nos cuenta más sobre el canal / experimento subyacente que está ejecutando, podemos brindarle más información sobre posibles rutas.
Spacey
@ Mohammad Gracias. El canal subyacente es el acero. ¿Alguna idea de cómo estimar el retraso de la energía?
Nick Sinas
@Mohammad, ¿crees que la distorsión de las señales podría ser algún tipo de reverberación que podría limpiarse con filtrado?
Nick Sinas
@NickS. Puede haber algunos trucos de limpieza de reverberación (no sé cómo se llevarían a cabo), pero he improvisado algo simple que será un estimador de energía si quieres echar un vistazo.
Spacey
0

Basado en la sugerencia de Muhammad, intenté hacer un guión de Matlab. Sin embargo, no puedo deducir si construye una distribución gaussiana basada en las variaciones y luego toma una estimación de KDE o si realiza una estimación de KDE con suposición gaussiana.

También es difícil inferir cómo traduce el tiempo de desplazamiento de KDE en el dominio del tiempo. Aquí está mi intento de hacerlo. Cualquier usuario que esté interesado en usar el script es libre y actualiza la versión mejorada si es posible.

%% Initialising data

Ws1 = data1;
Ws2 = data2;
mWs1 = nanmean(Ws1);
mWs2 = nanmean(Ws2);
sdWs1 = nanstd(Ws1);
sdWs2 = nanstd(Ws2);

%% Computing the signal envelopes
Ws1d = Ws1 - mWs1;
Ws2d = Ws2 - mWs2;
h1 = abs(hilbert(Ws1d));
h2 = abs(hilbert(Ws2d));
figure();
subplot(211)
plot([Ws1d, h1])
subplot(212)
plot([Ws2d, h2])

%% Denoise the signal with edge preserving nonlinear medial filtering
w = 25;
mf1 = medfilt1(h1, w);
mf2 = medfilt1(h2, w);
figure();
subplot(211)
plot(mf1)
subplot(212)
plot(mf2)

<%% Remove time: construct the gaussian kernel density estimation functions>
% Using the kde from Matlab central directly on the filtered data
data1 = mf1;
[bw1, den1, xmesh1, cdf1] = kde(data1, 2^14);
der1 = diff(den1);
data2 = mf2;
[bw2, den2, xmesh2, cdf2] = kde(data2, 2^14);
der2 = diff(den2);
figure();
plot([der1, der2]);
legend('Sig1', 'Sig2')

% the other method as explained in Muhammad's post
for i = 1:length(mf1)
gf1(:,i) = mf1(i) + sdWs1*randn(1000,1);
gf2(:,i) = mf2(i) + sdWs2*randn(1000,1);
end
[bwM1, denM1, xmeshM1, cdfM1] = kde(gf1(:,1), 2.^11);
dd1 = diff(denM1);
[bwM2, denM2, xmeshM2, cdfM2] = kde(gf2(:,1), 2.^11);
dd2 = diff(denM2);
figure();
plot([dd1, dd2]);
legend('Sig1', 'Sig2')
AshG
fuente