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.

MATLAB
v=depth_rate(P, sampling_rate = NULL, fc= NULL)
R
v <- depth_rate(P, sampling_rate = NULL, fc= NULL)
Inputs
  • P: The depth or altitude vector (a regularly sampled time series) or depth or altitude sensors list, sampled at sampling_rate Hz. p can have any units.
  • sampling_rate: The sampling rate of the sensor data in Hz (samples per second). This is only needed if p is not a sensor structure.
  • fc: (optional). Specifies the cut-off frequency or frequencies of the smoothing filter in Hz. If fc is not given, a default value is used of 0.2 Hz (i.e., a smoothing time constant of 5 seconds).
Outputs
  • v: vertical velocity with the same sampling rate as P. v has the same dimensions as P. The unit of v depends on the unit of P. If P is in meters, v is in meters/second.
Notes
  • The low-pass filter is a symmetric finite impulse response (FIR) filter with length 4fs/fc. The FIR filter removes the group delay of the filters.
MATLAB
s=ocdr(P, A, sampling_rate = NULL, fc= NULL, plim = NULL)
R
s <- ocdr(P, A, sampling_rate = NULL, fc= NULL, plim = NULL)
Inputs
  • P: The depth or altitude vector (a regularly sampled time series) or depth or altitude sensors list in meters, sampled at sampling_rate Hz.
  • A: The nx3 acceleration matrix with columns [ax ay az] or acceleration sensor list (e.g., from readtag.R). Acceleration can be in any consistent unit, e.g., g or m/s^2. A must have the same number of rows as P.
  • sampling_rate: The sampling rate of the sensor data in Hz (samples per second). This is only needed if p is not a sensor structure.
  • fc: (optional). (optional) Specifies the cut-off frequency of a low-pass filter to apply to P after computing depth-rate and to A before computing pitch. The filter cut-off frequency is in Hz. The filter length is 4*sampling_rate/fc. Filtering adds no group delay. If fc is empty or not given, the default value of 0.2 Hz (i.e., a 5 second time constant) is used.
  • plim: (optional) Specifies the minimum pitch angle in radians at which speed can be computed. Errors in speed estimation using this method increase strongly at low pitch angles. To avoid estimates with poor accuracy being used in later analyses, speed estimates at low pitch angles are replaced by NaN (not-a-number). The default threshold for this is 20 degrees.
Ouputs
  • s: is the forward speed estimate in m/s. s has the same size as P.
Notes
  • This function assumes that a positive pitch angle is an anti-clockwise rotation around the y-axis (lateral axis) follwiing the [north,east,up] navigation frame and a [forward,right,up] local frame. Thus, a descending animal willhave a negative pitch angle.
  • The low-pass filter is a symmetric finite impulse response (FIR) filter with length 4fs/fc. The FIR filter removes the group delay of the filters.

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..

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

Then use plott() to inspect it:

MATLAB
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'})
R
plott(X = list(Depth = bw$P, Acc = bw$Aa, Mag = bw$Ma))
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:

MATLAB
v = depth_rate(bw.P)
R
v <- depth_rate(bw$P)

Vertical speed will be an accurate estimation of the actual forward speed of the animal during periods of high pitch angle either during descent and ascent diving phases. If no steep periods exist (magnitude greater than 45 deg), vertical speed may not be an accurate estimation of the actual forward speed of the animal.

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.

MATLAB
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')
R
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)
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.

MATLAB
fc = 0.05 % cut-off frequency of the low-pass filter
v2 = depth_rate(bw.P, fc)
R
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.

MATLAB
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')
R
plott(X = list(`ver_speed_fc0.2` = v, `ver_speed_fc0.05` = v2), fsx = bw$A$sampling_rate)
MATLAB Output
R Output

It seems that v works well during the whole time except during the ascent part of the deep dive where a higher smoothing seems to do better. Thus, we should keep v2 which uses fc=0.05.

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:

MATLAB
fc = 0.05 % set the cut-off filter frequency. 
s = ocdr(bw.P,bw.A,fc)
R
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.

MATLAB
%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')
R
#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)
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:

MATLAB
plim=30*pi/180 % plim should be specified in radians
s2=ocdr(bw.P,bw.A,fc,plim)
R
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.

MATLAB
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');
R
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'))
MATLAB Output
R Output

Common Questions

Vertical speed and orientation corrected depth rate will be very similar during periods of high pitch angle either during descent and ascent diving phases. If no steep periods (greater than 30 degrees absolute pitch angles) orientation corrected depth rate will not be computed and vertical speed will not be an accurate estimation of the actual forward speed of the animal.

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