Main Content

La traducción de esta página está obsoleta. Haga clic aquí para ver la última versión en inglés.

Tomar derivadas de una señal

Quiere diferenciar una señal sin aumentar la potencia de ruido. La función diff de MATLAB® amplifica el ruido y la imprecisión resultante empeora con derivadas más altas. Para solucionar este problema, utilice un filtro diferenciador en su lugar.

Analice el desplazamiento del suelo de un edificio durante un terremoto. Encuentre la velocidad y la aceleración como funciones de tiempo.

Cargue el archivo earthquake. El archivo contiene las siguientes variables:

  • drift: desplazamiento del suelo, medido en centímetros

  • t: tiempo, medido en segundos

  • Fs: tasa de muestreo, igual a 1 kHz

load('earthquake.mat')

Utilice pwelch para mostrar una estimación del espectro de potencia de la señal. Observe cómo la mayoría de la energía de la señal se encuentra en frecuencias inferiores a 100 Hz.

pwelch(drift,[],[],[],Fs)

Figure contains an axes object. The axes object with title Welch Power Spectral Density Estimate contains an object of type line.

Utilice designfilt para diseñar un diferenciador FIR de orden 50. Para incluir la mayoría de la energía de la señal, especifique una frecuencia de banda de paso de 100 Hz y una frecuencia de banda de parada de 120 Hz. Inspeccione el filtro con fvtool.

Nf = 50; 
Fpass = 100; 
Fstop = 120;

d = designfilt('differentiatorfir','FilterOrder',Nf, ...
    'PassbandFrequency',Fpass,'StopbandFrequency',Fstop, ...
    'SampleRate',Fs);

fvtool(d,'MagnitudeDisplay','zero-phase','Fs',Fs)

Figure Filter Visualization Tool - Zero-phase Response contains an axes object and other objects of type uitoolbar, uimenu. The axes object with title Zero-phase Response contains 2 objects of type line.

Diferencie los desvíos para encontrar la velocidad. Divida la derivada entre dt, el intervalo de tiempo entre muestras consecutivas, para establecer las unidades correctas.

dt = t(2)-t(1);

vdrift = filter(d,drift)/dt;

La señal filtrada se retrasa. Utilice grpdelay para determinar que el retraso es la mitad del orden del filtro. Compénselo descartando muestras.

delay = mean(grpdelay(d))
delay = 25
tt = t(1:end-delay);
vd = vdrift;
vd(1:delay) = [];

La salida también incluye un transitorio cuya longitud es igual al orden del filtro o dobla el retraso del grupo. Las muestras delay se descartaron en el paso anterior. Descarte más delay para eliminar el transitorio.

tt(1:delay) = [];
vd(1:delay) = [];

Represente el desvío y la velocidad de desvío. Utilice findpeaks para verificar que los máximos y mínimos del desvío se corresponden con los cruces por cero de su derivada.

[pkp,lcp] = findpeaks(drift);
zcp = zeros(size(lcp));

[pkm,lcm] = findpeaks(-drift);
zcm = zeros(size(lcm));

subplot(2,1,1)
plot(t,drift,t([lcp lcm]),[pkp -pkm],'or')
xlabel('Time (s)')
ylabel('Displacement (cm)')
grid

subplot(2,1,2)
plot(tt,vd,t([lcp lcm]),[zcp zcm],'or')
xlabel('Time (s)')
ylabel('Speed (cm/s)')
grid

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

Diferencie la velocidad de desvío para encontrar la aceleración. El desfase es el doble de largo. Descarte el doble de muestras para compensar el retardo y el mismo número para eliminar el transitorio. Represente la velocidad y la aceleración.

adrift = filter(d,vdrift)/dt;

at = t(1:end-2*delay);
ad = adrift;
ad(1:2*delay) = [];

at(1:2*delay) = [];
ad(1:2*delay) = [];

subplot(2,1,1)
plot(tt,vd)
xlabel('Time (s)')
ylabel('Speed (cm/s)')
grid

subplot(2,1,2)
plot(at,ad)
ax = gca;
ax.YLim = 2000*[-1 1];
xlabel('Time (s)')
ylabel('Acceleration (cm/s^2)')
grid

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

Calcule la aceleración utilizando diff. Añada ceros para compensar el cambio del tamaño del arreglo. Compare el resultado con el obtenido con el filtro. Observe el aumento de ruido de alta frecuencia.

vdiff = diff([drift;0])/dt;
adiff = diff([vdiff;0])/dt;

subplot(2,1,1)
plot(at,ad)
ax = gca;
ax.YLim = 2000*[-1 1];
xlabel('Time (s)')
ylabel('Acceleration (cm/s^2)')
grid
legend('Filter')
title('Acceleration with Differentiation Filter')

subplot(2,1,2)
plot(t,adiff)
ax = gca;
ax.YLim = 2000*[-1 1];
xlabel('Time (s)')
ylabel('Acceleration (cm/s^2)')
grid
legend('diff')

Figure contains 2 axes objects. Axes object 1 with title Acceleration with Differentiation Filter contains an object of type line. This object represents Filter. Axes object 2 contains an object of type line. This object represents diff.

Consulte también

| | | |

Temas relacionados