Medir similitudes de señales

Comparar señales con diferentes tasas de muestreo

Considere una base de datos de señales de audio y 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

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 (s)")
linkaxes(ax(1:3),"x")
axis([0 1.61 -4 4])

La primera y la segunda subgráfica 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 una 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) antialiasing 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 sample rate by rational factor
T2 = resample(T2,P2,Q2);        % Change sample rate by rational factor

Encontrar una señal en una medición

Ahora se puede 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(s)") 
axis(ax(1:2),[-1.5 1.5 -700 700])

La primera subgráfica indica que la señal S y la plantilla T1 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 plantilla T2 presenta una señal S de 499 muestras, según indica SampleDiff. Esta información puede utilizarse para alinear las señales.

Medir el retardo entre señales y alinearlas

Considere una situación en la que recopila datos de distintos sensores que registran las vibraciones provocadas por vehículos a ambos extremos de un puente. Cuando analice las señales, puede que necesite alinearlas. Suponga que tiene 3 sensores que funcionan con las mismas tasas de muestreo y que miden las 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")

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 se puede utilizar para alinear las 3 señales desplazándolas en el tiempo. También podemos alinear las señales utilizando la función alignsignals, retrasando 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")

Comparar el contenido frecuencial de 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;         % Sample 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 (s)")
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)")

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

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

Considere 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

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

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

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 cíclico diario, con temperaturas más bajas durante la noche que aumentan durante el día.