Main Content

Medir similitudes de señales

Este ejemplo muestra cómo medir similitudes de señales. Le ayudará a responder a preguntas como las siguientes: ¿Cómo comparo señales con distintas longitudes o distintas tasas de muestreo? ¿Cómo averiguo si hay una señal o solo ruido en una medición? ¿Están relacionadas dos señales? ¿Cómo mido un retardo entre dos señales (y cómo las alineo)? ¿Cómo comparo el contenido frecuencial de dos señales? También pueden encontrarse similitudes en diferentes secciones de una señal para determinar si una señal es periódica.

Comparar señales con distintas tasas de muestreo

Piense en una base de datos de señales de audio y en una aplicación de concordancia de patrones en la que deba identificar una canción a medida que se reproduce. Los datos suelen almacenarse con una tasa de muestreo baja para ocupar menos memoria.

load relatedsig.mat

figure
ax(1) = subplot(3,1,1);
plot((0:numel(T1)-1)/Fs1,T1,'k')
ylabel('Template 1')
grid on
ax(2) = subplot(3,1,2); 
plot((0:numel(T2)-1)/Fs2,T2,'r')
ylabel('Template 2')
grid on
ax(3) = subplot(3,1,3); 
plot((0:numel(S)-1)/Fs,S)
ylabel('Signal')
grid on
xlabel('Time (secs)')
linkaxes(ax(1:3),'x')
axis([0 1.61 -4 4])

Figure contains 3 axes objects. Axes object 1 contains an object of type line. Axes object 2 contains an object of type line. Axes object 3 contains an object of type line.

La primera subgráfica y la segunda muestran las señales de plantilla de la base de datos. La tercera subgráfica muestra la señal que deseamos buscar en la base de datos. Solo observando la serie de tiempo, la señal no parece coincidir con ninguna de las dos plantillas. Un examen más detallado revela que, en realidad, las señales tienen diferentes longitudes y tasas de muestreo.

[Fs1 Fs2 Fs]
ans = 1×3

        4096        4096        8192

Las diferentes longitudes impiden calcular la diferencia entre dos señales, pero esto se puede solucionar fácilmente extrayendo la parte común de las señales. Además, no siempre es necesario ecualizar longitudes. Se puede realizar correlación cruzada entre señales con longitudes diferentes, pero es fundamental asegurarse de que tienen tasas de muestreo idénticas. La forma más segura de hacerlo consiste en remuestrear la señal con una tasa de muestreo más baja. La función resample aplica un filtro FIR (paso bajo) anti-aliasing a la señal durante el proceso de remuestreo.

[P1,Q1] = rat(Fs/Fs1);          % Rational fraction approximation
[P2,Q2] = rat(Fs/Fs2);          % Rational fraction approximation
T1 = resample(T1,P1,Q1);        % Change sampling rate by rational factor
T2 = resample(T2,P2,Q2);        % Change sampling rate by rational factor

Encontrar una señal en una medición

Ahora podemos correlacionar en forma cruzada la señal S con las plantillas T1 y T2 utilizando la función xcorr para determinar si hay una coincidencia.

[C1,lag1] = xcorr(T1,S);        
[C2,lag2] = xcorr(T2,S);        

figure
ax(1) = subplot(2,1,1); 
plot(lag1/Fs,C1,'k')
ylabel('Amplitude')
grid on
title('Cross-correlation between Template 1 and Signal')
ax(2) = subplot(2,1,2); 
plot(lag2/Fs,C2,'r')
ylabel('Amplitude') 
grid on
title('Cross-correlation between Template 2 and Signal')
xlabel('Time(secs)') 
axis(ax(1:2),[-1.5 1.5 -700 700 ])

Figure contains 2 axes objects. Axes object 1 with title Cross-correlation between Template 1 and Signal contains an object of type line. Axes object 2 with title Cross-correlation between Template 2 and Signal contains an object of type line.

La primera subgráfica indica que la señal y la plantilla 1 están menos correlacionadas, mientras que el pico alto de la segunda subgráfica indica que la señal está presente en la segunda plantilla.

[~,I] = max(abs(C2));
SampleDiff = lag2(I)
SampleDiff = 499
timeDiff = SampleDiff/Fs
timeDiff = 0.0609

El pico de la correlación cruzada implica que la señal está presente en la plantilla T2 después de 61 ms. En otras palabras, la señal T2 presenta un desfase positivo de 499 muestras con respecto a la señal S, según indica SampleDiff. Esta información puede utilizarse para alinear las señales.

Medir un retardo entre señales y alinearlas

Considere una situación en la que recopila datos de distintos sensores que graban vibraciones provocadas por vehículos a ambos extremos de un puente. Cuando analice las señales, puede que necesite alinearlas. Suponga que tiene tres sensores que funcionan con las mismas tasas de muestreo y que miden señales que provoca el mismo evento.

figure
ax(1) = subplot(3,1,1);
plot(s1)
ylabel('s1')
grid on
ax(2) = subplot(3,1,2); 
plot(s2,'k')
ylabel('s2')
grid on
ax(3) = subplot(3,1,3); 
plot(s3,'r')
ylabel('s3')
grid on
xlabel('Samples')
linkaxes(ax,'xy')

Figure contains 3 axes objects. Axes object 1 contains an object of type line. Axes object 2 contains an object of type line. Axes object 3 contains an object of type line.

También podemos utilizar la función finddelay para encontrar el retardo entre dos señales.

t21 = finddelay(s1,s2)
t21 = -350
t31 = finddelay(s1,s3)
t31 = 150

t21 indica que s2 presenta un desfase negativo con respecto a s1 de 350 muestras y t31 indica que s3 presenta un desfase positivo con respecto a s1 de 150 muestras. Ahora, esta información puede utilizarse para alinear las tres señales desplazando las señales en el tiempo. También podemos alinear las señales directamente utilizando la función alignsignals, que alinea dos señales retardando la primera señal.

s1 = alignsignals(s1,s3);
s2 = alignsignals(s2,s3);

figure
ax(1) = subplot(3,1,1);
plot(s1)
grid on 
title('s1')
axis tight
ax(2) = subplot(3,1,2);
plot(s2)
grid on 
title('s2')
axis tight
ax(3) = subplot(3,1,3); 
plot(s3)
grid on 
title('s3')
axis tight
linkaxes(ax,'xy')

Figure contains 3 axes objects. Axes object 1 with title s1 contains an object of type line. Axes object 2 with title s2 contains an object of type line. Axes object 3 with title s3 contains an object of type line.

Comparar el contenido frecuencial de las señales

Un espectro de potencia muestra la potencia presente en cada frecuencia. La coherencia espectral identifica la correlación frecuencia-dominio entre señales. Los valores de coherencia que tienden a 0 indican que los componentes de frecuencia correspondientes no están correlacionados, mientras que los valores que tienden a 1 indican que los componentes de frecuencia correspondientes están correlacionados. Considere dos señales y sus espectros de potencia respectivos.

Fs = FsSig;         % Sampling Rate

[P1,f1] = periodogram(sig1,[],[],Fs,'power');
[P2,f2] = periodogram(sig2,[],[],Fs,'power');

figure
t = (0:numel(sig1)-1)/Fs;
subplot(2,2,1)
plot(t,sig1,'k')
ylabel('s1')
grid on
title('Time Series')
subplot(2,2,3)
plot(t,sig2)
ylabel('s2')
grid on
xlabel('Time (secs)')
subplot(2,2,2)
plot(f1,P1,'k')
ylabel('P1')
grid on
axis tight
title('Power Spectrum')
subplot(2,2,4)
plot(f2,P2)
ylabel('P2')
grid on
axis tight
xlabel('Frequency (Hz)')

Figure contains 4 axes objects. Axes object 1 with title Time Series contains an object of type line. Axes object 2 contains an object of type line. Axes object 3 with title Power Spectrum contains an object of type line. Axes object 4 contains an object of type line.

La función mscohere calcula la coherencia espectral entre las dos señales. Confirma que sig1 y sig2 tienen dos componentes correlacionados en torno a 35 Hz y 165 Hz. En las frecuencias en las que la coherencia espectral es alta, la fase relativa entre los componentes correlacionados se puede estimar con la fase de espectro cruzado.

[Cxy,f] = mscohere(sig1,sig2,[],[],[],Fs);
Pxy = cpsd(sig1,sig2,[],[],[],Fs);
phase = -angle(Pxy)/pi*180;
[pks,locs] = findpeaks(Cxy,'MinPeakHeight',0.75);

figure
subplot(2,1,1)
plot(f,Cxy)
title('Coherence Estimate')
grid on
hgca = gca;
hgca.XTick = f(locs);
hgca.YTick = 0.75;
axis([0 200 0 1])
subplot(2,1,2)
plot(f,phase)
title('Cross-spectrum Phase (deg)')
grid on
hgca = gca;
hgca.XTick = f(locs); 
hgca.YTick = round(phase(locs));
xlabel('Frequency (Hz)')
axis([0 200 -180 180])

Figure contains 2 axes objects. Axes object 1 with title Coherence Estimate contains an object of type line. Axes object 2 with title Cross-spectrum Phase (deg) contains an object of type line.

El desfase entre los componentes de 35 Hz se acerca a -90 grados y el desfase entre los componentes de 165 Hz, a -60 grados.

Encontrar periodicidades en una señal

Piense en un conjunto de mediciones de la temperatura en un edificio de oficinas durante la temporada de invierno. Se realizaron mediciones cada 30 minutos durante unas 16,5 semanas.

load officetemp.mat  

Fs = 1/(60*30);                 % Sample rate is 1 sample every 30 minutes
days = (0:length(temp)-1)/(Fs*60*60*24); 

figure
plot(days,temp)
title('Temperature Data')
xlabel('Time (days)')
ylabel('Temperature (Fahrenheit)')
grid on

Figure contains an axes object. The axes object with title Temperature Data contains an object of type line.

Con temperaturas en valores bajos del tramo de los 70 grados Fahrenheit, tiene que eliminar la media para analizar pequeñas fluctuaciones en la señal. La función xcov elimina la media de la señal antes de calcular la correlación cruzada. Devuelve la covarianza cruzada. Limite el desfase máximo al 50% de la señal para obtener una buena estimación de la covarianza cruzada.

maxlags = numel(temp)*0.5;
[xc,lag] = xcov(temp,maxlags);         

[~,df] = findpeaks(xc,'MinPeakDistance',5*2*24);
[~,mf] = findpeaks(xc);

figure
plot(lag/(2*24),xc,'k',...
     lag(df)/(2*24),xc(df),'kv','MarkerFaceColor','r')
grid on
xlim([-15 15])
xlabel('Time (days)')
title('Auto-covariance')

Figure contains an axes object. The axes object with title Auto-covariance contains 2 objects of type line.

Observe las fluctuaciones dominantes y menores en la covarianza automática. Los picos dominantes y menores parecen ser equidistantes. Para verificar si lo son, calcule y represente la diferencia entre las ubicaciones de los picos posteriores.

cycle1 = diff(df)/(2*24);
cycle2 = diff(mf)/(2*24);

subplot(2,1,1)
plot(cycle1)
ylabel('Days')
grid on
title('Dominant peak distance')
subplot(2,1,2)
plot(cycle2,'r')
ylabel('Days')
grid on
title('Minor peak distance')

Figure contains 2 axes objects. Axes object 1 with title Dominant peak distance contains an object of type line. Axes object 2 with title Minor peak distance contains an object of type line.

mean(cycle1)
ans = 7
mean(cycle2)
ans = 1

Los picos menores indican 7 ciclos/semana y los picos dominantes, 1 ciclo/semana. Esto tiene sentido, dado que los datos proceden de un edificio con temperatura controlada siguiendo un calendario de 7 días. El primer ciclo de 7 días indica que hay un comportamiento semanal cíclico de la temperatura del edificio, con temperaturas más bajas durante los fines de semana que vuelven a la normalidad durante los días hábiles. El comportamiento del ciclo de 1 día indica que también hay un comportamiento diario cíclico, con temperaturas más bajas durante la noche que aumentan durante el día.

Consulte también

| | | | | |