# 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 =
```
R
```bw_file_path <- system.file("extdata", "testset1.nc", package = 'tagtools', mustWork = TRUE)
```

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

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