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 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.
Consulte también
alignsignals
| cpsd
| finddelay
| findpeaks
| mscohere
| xcov
| xcorr