Estimate Vertical Speed and Forward Speed

Usage

In this vignette you will learn to use depth_rate() and ocdr() to gain insight about how to estimate vertical and forward speed, respectively.

v=depth_rate(P, sampling_rate = NULL, fc= NULL)

v <- depth_rate(P, sampling_rate = NULL, fc= NULL)

Inputs

Outputs

Notes

s=ocdr(P, A, sampling_rate = NULL, fc= NULL, plim = NULL)

s <- ocdr(P, A, sampling_rate = NULL, fc= NULL, plim = NULL)

Inputs

Ouputs

Notes

Theoretical Description

Some tags have a speed sensor but most of them not, however, vertical and forward speed can still be estimated. Here we will present two different estimates that can be obtained with only depth data or a combination of depth and acceleration data.

Technical Description

In this vignette, you will learn how to estimate vertical velocity by differentiating a depth or altitude time series using the depth_rate() function. A low-pass symmetric filter reduces the sensor noise that is amplified by the differentiation.

You will also learn how to estimate forward speed of a flying or diving animal by first computing the altitude or depth-rate (i.e., the first differential of the pressure in meters) and then correcting for the pitch angle using the ocdr() function. This is called the Orientation Corrected Depth Rate.

Load and Visualize the Test Data Set

Load the test dataset testset1, which belongs to the data recorded from a suction cup tag attached to the back of a beaked whale [Note]This dataset has already been converted from a source file that was offloaded from the tag into a NetCDF file. In doing so, some metadata was gleaned from the file and added to the data. Other metadata was added by hand..

bw_file_path = 
bw = load_nc(bw_file_path)

bw_file_path <- system.file("extdata", "testset1.nc", package = 'tagtools', mustWork = TRUE)
bw <- load_nc(bw_file_path)

Then use plott() to inspect it:

plott(bw.P, bw.Aa, bw.Ma)
subplot(311)
ylabel('Depth')
subplot(312)
ylabel('Acc')
subplot(313)
ylabel('Mag')
legend({'x-axis','y-axis','z-axis'})

plott(X = list(Depth = bw$P, Acc = bw$Aa, Mag = bw$Ma))

Expand to show figure …

• MATLAB Output
• R Output

MATLAB Output

R Output

If your depth data needs to be corrected for off animal periods, outliers, pressure offsets and temperature effects, please go to and perform the tutorial on depth checks and corrections before proceeding to estimate the vertical or the forward speed.

Estimate vertical speed

We are going to estimate vertical speed using the depth_rate() function. This function estimates vertical speed as the differential of depth. This differentiation amplifies any present depth sensor noise, and so depth_rate() function smooths the depth rate signal by default (with a smoothing filter cut off frequency of 0.2 Hz, this is a smoothing time constant of 5 seconds) to remove this noise.

Run depth_rate() on bw.P to estimate the vertical speed:

v = depth_rate(bw.P)

v <- depth_rate(bw$P)

Plot vertical speed along with the pitch angle and depth to get a better insight. First you need to compute the pitch from the accelerometer data on the animal frame using a2pr() function. Go to the animal orientation tutorial for a greater detail on how to estimate pitch.

pitch = a2pr (bw.A.data) % estimate pitch angle in radians
fs = P.sampling_rate

figure
plott(P,fs,pitch*180/pi,fs,v,fs)
%add some labels 
subplot(311)
ylabel('Depth m')
subplot(312)
ylabel('Pitch degrees')
subplot(313)
ylabel('Vertical speed, m/s')

pitch <- a2pr (bw$A$data) # estimate pitch angle in radians
fs <- bw$A$sampling_rate

plott(X = list(`Depth (m)` = bw$P, `Pitch (degrees)` = pitch$p*180/pi, `Vertical speed, m/s`=v),
      fsx = fs, interactive = TRUE)

Expand to show figure …

• MATLAB Output
• R Output

MATLAB Output

R Output

Note that pitch is not a smooth signal but contains fluke-stroking information. In order to smooth the pitch signal have a look at the animal orientation tutorial.

If you find that a higher smoothing is needed, set the cut-off frequency of the low-pass filter smaller (i.e., fc=0.05, which is a smoothing time constant of 20 seconds), so that higher frequencies are blocked by the filter.

fc = 0.05 % cut-off frequency of the low-pass filter
v2 = depth_rate(bw.P, fc)

fc <- 0.05 # cut-off frequency of the low-pass filter
v2 <- depth_rate(p=bw$P, fs=fs,fc=fc)

Plot v and v2 to see the effect of a higher smoothing in the vertical speed estimation and which one would be better to use.

figure
plot((1:length(v))/fs,v)
hold on 
plot((1:length(v2))/fs,v2, 'red')
ylabel('Vertical speed, m/s')
legend('fc = 0.2','fc = 0.05')

plott(X = list(`ver_speed_fc0.2` = v, `ver_speed_fc0.05` = v2), fsx = bw$A$sampling_rate)

Expand to show figure …

• MATLAB Output
• R Output

MATLAB Output

R Output

Estimate forward speed

We are going to estimate vertical speed using the ocdr() function. This function estimates the forward speed by first differentiating the depth data and then correcting it for the pitch angle. This is why this is also called the orientation corrected depth rate. This function assumes that (i) the animal moves in the direction of its longitudinal axis (which is normally the case), and (ii) that a descending animal will have a negative pitch angle while an ascending animal will have a positive pitch angle.

Similar to the depth_rate() function, ocdr() smooths the depth rate by default (with a smoothing low-pass filter cut off frequency of 0.2 Hz, this is a smoothing time constant of 5 seconds) to remove the depth sensor noise. In addition ocdr also smooths the accelerometer signal by default (with a smoothing low-pass filter cut off frequency of 0.2 Hz, this is a smoothing time constant of 5 seconds)  before estimating the pitch angle (see documentation of a2pr() function for a greater detail and better understanding).

Note that we found out that a cut-off frequency of 0.05 was a better smoothing when running depth_rate(). Thus we should use fc = 0.05 when running ocdr() on bw.P and bw.A to estimate the forward speed or orientation corrected depth rate:

fc = 0.05 % set the cut-off filter frequency. 
s = ocdr(bw.P,bw.A,fc)

fc <- 0.05 # set the cut-off filter frequency. 
s <- ocdr(p=bw$P$data,A= bw$A$data, sampling_rate = bw$A$sampling_rate,fc=fc, plim = NULL)

Let’s plot forward speed along with the smooth pitch angle and depth to get a better insight. First you need to compute the pitch from the accelerometer data on the animal frame using a2pr() function. Go to the animal orientation tutorial for a greater detail on how to estimate pitch.

%Estimate smooth pitch
nf = round(4*fs/fc) ;
Alf = fir_nodelay(A,nf,fc/(fs/2)) ;% estimate the low frequency acceleration
smoothpitch = a2pr(Alf); 

%plot depth, smoothpich in degrees and ocdr in m/s
figure
plott(P,fs,smoothpitch*180/pi,fs,s,fs)
%add some labels 
subplot(311)
ylabel('Depth m')
subplot(312)
ylabel('Pitch deg')
subplot(313)
ylabel('ocdr m/s')

#Estimate smooth pitch
nf <- round(4*fs/fc) ;
Alf <- fir_nodelay(x=bw$A$data,n=nf,fc=fc/(fs/2),qual="low"); # estimate the low frequency acceleration
smoothpitch <- a2pr(Alf); 

#plot depth, smoothpich in degrees and ocdr in m/s

plott(X = list(`Depth (m)` = bw$P, `Pitch (degrees)` = smoothpitch$p*180/pi, `ocdr, m/s`=s),fsx = fs, interactive = TRUE)

Expand to show figure …

• MATLAB Output
• R Output

MATLAB Output

R Output

Errors in speed estimation using this method increase strongly at low pitch angles. However, to avoid estimates with poor accuracy being used in later analyses, speed estimates at low pitch angles are replaced automatically by NaN (not-a-number) in ocdr(). The default ocdr plim is 20 degrees. Note that ocdr was not plotted for absolute pitch angles smaller than 20 degrees.

However, based on previous publications we should set the pitch limit (plim) to 30 degrees (see Miller et al., 2004; Simon et al., 2012) for better estimates. You could do so by setting plim = 30*pi/180:

plim=30*pi/180 % plim should be specified in radians
s2=ocdr(bw.P,bw.A,fc,plim)

plim <- 30*pi/180 # plim should be specified in radians
s2 <- ocdr(p=bw$P,A=bw$A,sampling_rate=fs,fc=fc,plim=plim)

Compare forward and vertical speed

Vertical speed is the is the differential of the depth rate and forward speed or orientation corrected depth rate is the differential of the depth rate in meters corrected for the pitch angle.

Let’s plot vertical speed, forward speed and also the smooth pitch angle. Note that forward speed is not plotted for pitch angles smaller than 30 degrees. To plot the smooth pitch signal we have added a few lines here. However, for a full description on how to do it, please refer to the animal orientation tutorial.

figure
ax1 = subplot(3,1,1);
plot((1:length(bw.P.data))/fs/60,bw.P.data)
axis ij
grid on
ylabel('Depth m')

ax2 = subplot(3,1,2);
plot((1:length(smoothpitch))/fs/60,smoothpitch*180/pi)
grid on
ylabel('Pitch deg')

ax3 = subplot(3,1,3);
plot((1:length(v2))/fs/60,v2)
hold on 
plot((1:length(s))/fs/60,s2,'r')
grid on
ylabel('speed m/s')
legend(' vertical speed','forward speed')
xlabel('Time (min)')

linkaxes([ax1 ax2 ax3], 'x');

orig_graph_settings <- par(# get 3 rows and 1 column of subplots
    mfrow = c(3,1), 
    # set margins
    mar = c(4.1, 4.1, 1.1, 1.1))
# ax1
plot(x = c(1:length(bw$P$data)) / bw$P$sampling_rate / 60,
     y = bw$P$data,
     ylim = rev(range(bw$P$data, na.rm = TRUE)),
     ylab = 'Depth (m)',
     xlab = '',
     type = 'l')
# equivalent to grid on isn't available; in ggformula/ggplot you change the "theme" the graph is using and in base R graphics you can use the function grid() but have to specify nx and ny, number of cells in each dimension
# ax2
plot(x = c(1:length(smoothpitch$p)) / fs / 60,
     y = smoothpitch$p*180/pi,
     ylab = 'Smoothed Pitch (deg)',
     xlab = '',
     type = 'l')
# ax3
plot(x = c(1:length(v2))/fs/60,
     y = v2,
     type = 'l',
     ylab = 'Speed (m/sec)',
     xlab = 'Time (min)')
lines(x = c(1:length(s2))/fs/60,
      y = s2,
      col = 'darkred',
      type = 'l')
legend("topright", bty = 'n', # for no box around legend
       col = c('black', 'darkred'),
       lty = 1, # for solid lines in legend
       lwd = 1, # width of lines in legend
       legend = c('vertical speed','forward speed'))

Expand to show figure …

• MATLAB Output
• R Output

MATLAB Output

R Output

Common Questions

Further Resources

For additional information on some of the methods presented here, please refer to the depth checks and corrections and animal orientation tutorials.

References

Miller, P. J. O., Johnson, M. P., Tyack, P. L., & Terray, E. A. (2004). Swimming gaits, passive drag and buoyancy of diving sperm whales Physeter macrocephalus. Journal of Experimental Biology, 207, 1953–1967. https://doi.org/10.1242/jeb.00993

Simon, M., Johnson, M., & Madsen, P. T. (2012). Keeping momentum with a mouthful of water: behavior and kinematics of humpback whale lunge feeding. The Journal of Experimental Biology, 215(Pt 21), 3786–3798. https://doi.org/10.1242/jeb.071092