Complementary Filter and Applications in Animal Movement

Usage

For the full theoretical description, refer to complementary filter – technical explanation.

Xf=comp_filt(X, sampling_rate = NULL, fc)

Xf<-comp_filt(X, sampling_rate = NULL, fc)

Inputs

Ouputs

Theoretical Description

Locomotion in many animals involves the oscillatory movement of limbs, fins or wings. These oscillations tend to occur at the dominant stroke frequency (dsf) and are overlaid on the constantly changing postural orientation, which tend to occur at a lower rate, as the animal manoeuvres. To separate these two signals that occur at a different rate, we can use a complementary filter using the comp_filt() function. For the full theoretical description, refer to complemetary filter – technical explanation.

Technical Description

In this vignette, you will learn how to use complementary filters, using the comp_filt() function, to separate intervals of movements into distinct frequency bands, i.e., animal orientation (low frequency), gait (mid frequency or dsf) and strikes/flinches (high frequency), in order to gain insight about the locomotion and biomechanical strategies of different aerial, terrestrial and marine animals. For the full technical description, refer to complementary filter – technical explanation.

Verification

To verify that you have run comp_filt()correctly, i.e., that you have set the filter cut off frequencies correctly, to seperate orientation from gait from strikes you will need to plot each signal separately and check that they look as they should. For the full verification, refer to complementary filter – technical explanation.

Caveats

As with any filter, if the cut-off frequency to separate orientation from gait locomotion is too low, the magnitude of the gait locomotion will be overstimated and the animal orientation will be mistakenly identified as strokes. On the other hand is this cut-off frequency is too high, the magnitude of the gait locomotion will be understimated and the stroking will be visible in the animal orientation, i.e., low-frequency signal. For the full verification, refer to complementary filter – technical explanation.

Load and Visualize the Test Data Set

Load the test dataset Brydes_test, which belongs to the data recorded from a suction cup tag attached to the back of a Bryde’s whale published in Izadi et al. (2022) [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.. This dataset is built into the tagtools package, so you can access it using system.file.

beb_file_path = 
beb = load_nc(beb_file_path)

beb_file_path <- system.file("extdata", "XXX", package = 'tagtools', mustWork = TRUE)
beb <- load_nc(beb_file_path)

Then use plott() to inspect it:

plott(beb.P, beb.Aa, beb.Ma)
legend({'x-axis','y-axis','z-axis'})

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

This dataset contains few consecutive surface lunges. We want to describe the biomechanical locomotion of this animal, and so we will separate the animal orientation (low frequency) from the stride, flap or fluke motion (mid frequency) from the rapid muscle movements and unsteady flow (high-frequency).

For tags not glued, attached with suction cups, it could happen that the tag is not properly attached (2 out of 4 suction cups attached) or that 2 out of 4 suction cups get detached earlier than the other 2 and so the tag vibrates. Stroking periods of the animal deployment should resemble a smooth sine wave in the acceleration data. However, if this are spiky such as in the figure below this may be an indication that the tag is vibrating, and we should try to remove that noise. If you need to do this on your data, please check the Remove high frequency noise due to tag vibration tutorial before you continue.

Complementary filter

To separate slow changes in orientation from stride, flap, or fluke motion dynamics in locomotion, and from rapid muscle movements we need to choose a filter cut-off frequency below the stride, flap, or fluke stroke frequency. We can estimate the dominant stroke frequency using dsf(). Go to the dominant stroke frequency tutorial for a greater detail on how to estimate it and use the  dsf() function.

Run dsf() on beb.Aa to get the mean stroking rate:

dsfa = dsf(beb.Aa) % estimated stroking rate in Hz from the accelerometer data

dsfa <- dsf(beb$Aa)$fpk   # estimated stroking rate in Hz

dsfa is 0.35 Hz (from Izadi et al., 2022).

A good starting choice for the filter cut-off frequency of the complementary filter to separate low and mid frequency acceleration movements is 70% of the dominant stroking rate, but cut-off frequencies of 50 to 70% the dsf could be good options. Call this value fc1. Depending on your reserach question you should use a different cut-off frequency. In fact, Izadi et al. (2022) used a cut-off frequency of 0.17 which is approximately 50% of the dsf during lunges.

To further separate the mid from the high frequency accelerometer data we need to set a second cut-off frequency filter, fc2. To get all the information regarding stroking in the mid frequency acceleration signal, fc2 should be at least twice the stroking rate of interest. In this case we will set fc2 to 0.8Hz being approximately twice the stroking rate during surface lunges.

Run a complementary filter on Aa to separate the slow, mid and fast time-scales. Recall that Aa is stored under bw as beb.Aa.

fc = 'YourValueHere'% your value for fc in Hz, a number, without quotes
fc1 = 0.25 %this is approximately 70% of the dsf (0.35 Hz)
fc2 = 0.8 % this is approximately twice the dsf (0.35 Hz)
Af  = comp_filt(beb.Aa, fc1,fc2)

fc <- 'YourValueHere' # your value for fc in Hz, a number, without quotes
fc1 <- 0.25 #this is approximately 70% of the dsf (0.35 Hz)
fc2 <- 0.8 # this is approximately twice the dsf (0.35 Hz)
Af <- comp_filt(beb$Aa, fc = fc1,fc2) 
str(Af, max.level = 1) 

The complementary filter returns a structure containing three data matrices: the low-pass filtered data, mid-pass filtered data and the high-pass filtered data. Each of these is a three-column matrix because the filter is run on each column of the input data. So, it is like you now have three accelerometers in the tag: one is only sensitive to low frequencies, another is only sensitive to mid frequencies, and the other is only sensitive to high frequencies. If you would like to get each matrix out of the cell array, do:

Alow =  Af{1} % low frequency A data
Amid =  Af{2} % mid frequency A data
Ahigh = Af{3} % high frequency A data

Alow <- Af[[1]]     # low frequency A data
Amid <- Af[[2]]     # mid frequency A data
Ahigh <- Af[[3]]    # high frequency A data

The sampling rate of these is the same as for the original data. For simplicity, make a variable sampling_rate equal to the sampling rate and use plott() to plot the three filtered accelerations along with the dive profile:

sampling_rate=beb.P.sampling_rate;
plott(beb.P.data,sampling_rate,Alow,sampling_rate,Ahigh,sampling_rate)
axis ij

sampling_rate <- beb$A$sampling_rate
plott(X = list(`Depth (m)` = beb$P$data,
               `LF Accel` = Alow, 
               `MF Accel` = Amid, 
               `HF Accel` = Ahigh),
      fsx = sampling_rate)

The low and mid versions of acceleration are sometimes called ‘static’ and ‘dynamic’ acceleration. If the filtering worked, you should see that Alow has large, relatively slow changes in acceleration which are mostly to do with orientation. These are missing in Amid and Ahigh. Amid has the mid frequency specific acceleration, having to do with propulsion (stride, fluke or flap), and Ahigh has the high frequency acceleration which has to do with rapid muscle movements such as strikes and flinches, and unsteady flow.

You could estimate animal orientation, pitch and roll of the animal from this low frequency acceleration. If you want to also estimate the heading of the animal you need to get the low frequency magnetometer data in the same way as you did for the accelerometer. Go through the animal orientation technical explanation for greater detail on how to estimate pitch, roll and heading and use the  a2pr() and m2h() functions.

Further Resources

For a more in-depth explanation of the usage and inner workings of the comp_filt() tool, refer to the Technical Explanation of the complementary filter.

References

Izadi, S., Aguilar de Soto, N., Constantine, R., & Johnson, M. (2022). Feeding tactics of resident Bryde’s whales in New Zealand. Marine Mammal Science, 1–14. https://doi.org/10.1111/mms.12918