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ímetrost
: tiempo, medido en segundosFs
: 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)
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)
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
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
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')
Consulte también
findpeaks
| FVTool | designfilt
| grpdelay
| periodogram