www.gusucode.com > signal 案例源码程序 matlab代码 > signal/TakeDerivativesOfASignalExample.m

    %% Take Derivatives of a Signal
% You want to differentiate a signal without increasing the noise power.
% MATLAB(R)'s function |diff| amplifies the noise, and the resulting
% inaccuracy worsens for higher derivatives. To fix this problem, use a
% differentiator filter instead.
%
% Analyze the displacement of a building floor during an earthquake. Find
% the speed and acceleration as functions of time.

%%
% Load the file |earthquake|. The file contains the following variables:
% 
% * |drift|: Floor displacement, measured in centimeters
% * |t|: Time, measured in seconds
% * |Fs|: Sample rate, equal to 1 kHz

load(fullfile(matlabroot,'examples','signal','earthquake.mat')) 

%%
% Use |pwelch| to display an estimate of the power spectrum of the signal.
% Note how most of the signal energy is contained in frequencies below 100
% Hz.

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

%%
% Use |designfilt| to design an FIR differentiator of order 50. To include
% most of the signal energy, specify a passband frequency of 100 Hz and a
% stopband frequency of 120 Hz. Inspect the filter with |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)

%%
% Differentiate the drift to find the speed. Divide the derivative by |dt|,
% the time interval between consecutive samples, to set the correct units.

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

vdrift = filter(d,drift)/dt;

%%
% The filtered signal is delayed. Use |grpdelay| to determine that the
% delay is half the filter order. Compensate for it by discarding samples.

delay = mean(grpdelay(d))

tt = t(1:end-delay);
vd = vdrift;
vd(1:delay) = [];

%%
% The output also includes a transient whose length equals the filter
% order, or twice the group delay. |delay| samples were discarded above.
% Discard |delay| more to eliminate the transient.

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

%%
% Plot the drift and the drift speed. Use |findpeaks| to verify that the
% maxima and minima of the drift correspond to the zero crossings of its
% derivative.

[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

%%
% Differentiate the drift speed to find the acceleration. The lag is twice
% as long. Discard twice as many samples to compensate for the delay, and
% the same number to eliminate the transient. Plot the speed and
% acceleration.

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

%%
% Compute the acceleration using |diff|. Add zeros to compensate for the
% change in array size. Compare the result to that obtained with the
% filter. Notice the amount of high-frequency noise.

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