Theoretical description
What is it?
Locomotion in many animals involves the oscillatory movement of limbs, fins or wings. The idea behind dominant stroke frequency is that animals tend to prefer a particular oscillation rate during routine transport. This rate may be determined by the size, weight and spring constant of anatomical structures, i.e., coinciding with a resonance or pendulum frequency, making it an efficient rate at which to locomote. Although the precise DSF may vary by individual, DSF is often broadly similar across animals with similar size and locomotion style, leading to general allometric relationships between body size and DSF. For example, stroke cycle frequencies of swimming specialist seabirds and marine mammals have been seen to be proportional to mass^0.29 (Sato et al., 2007). However, there can also be significant variation within and across species, e.g., some individuals tend to adopt faster or slower locomotion rates, known as the pace of life syndrome. In addition it is quite common that animals have a different dsf during the descent than the ascent phase, as an animal that encounters large forces (due to their buoyancy; Miller et al., 2004; MacAyeal et al., 2011) may respond by stroking faster.
How is it measured?
Locomotion generally involves the actuation of some of the larger muscles in the body. This produces specific acceleration oscillations as well as small cyclical changes in body orientation, which can be sensed potentially throughout the body. The presence of both acceleration and orientation changes at the locomotion rate means that any of the usual fine-scale motion sensors can be used to measure the DSF. Magnetometers and gyroscopes sense orientation, while accelerometers sense both orientation and specific acceleration, but any of these sensors in a surface-attached tag will produce components that vary at the locomotory rate. The magnitude of these locomotion signals depends on the position of the tag on the body but the oscillation rate, and therefore the DSF, does not.
Although the DSF can be estimated by finding each locomotory cycle in the motion sensor data (i.e., in the accelerometer, magnetometer or gyroscope signals), this would be a lot of work. A simpler approach is to compute a spectral average of the sensor data using Fast Fourier Transforms (FFTs). Locomotory motions are often largest in one or two of the animal’s three cardinal axes. For example, in cetacean swimming, locomotion signals are largest in the dorso-ventral and longitudinal axes. If the axis with the strongest locomotion signal is known, the spectral average can be computed on just this axis. Alternatively, the spectral averages of each axis can be combined, or can be inspected individually for a peak corresponding to the DSF. In either case, the DSF peak is often readily identified.
What is it good for?
The DSF may be a useful parameter in itself, e.g., to evaluate activity budgets and to identify intervals when the animal exceeds its normal locomotory rate (perhaps due to the need for additional force or speed). The DSF can also be useful to guide the design of high-pass or low-pass filters to separate postural components from specific acceleration components in on-animal acceleration data, e.g., in the calculation of ODBA.
Technical description
Usage of tool
The dsf() function calculates the spectrum level of the signal (X), i.e., the amount of power per 1 Hz band. Accelerometer and magnetometer signals have three axes each, so the spectrum level is estimated for each axis. Then, dsf() sums the spectral power in the three axes, and returns the frequency with the maximum power as an estimate of the dominant stroke frequency (fpk, i.e., the peak frequency in the sum of the power spectrum for each axis in Hz). This function also provides the quality of the peak (q) measured by the peak power divided by the mean power of the spectra. This is a dimensionless number which is large if a clear spectral peak exists.
To calculate the spectrum level, the length of the Fast Fourier Transform (FFT), and therefore the frequency resolution, can be specified by the user (Nfft). If not specified, the default value is the power of two closest to 20*fs, i.e., an analysis block length of about 20 seconds and a frequency resolution of about 0.05 Hz. A shorter FFT may be required if movement behavior is highly variable. In contrast, a longer FFT may work well if propulsion is continuous and stereotyped. fs is the sampling rate of the sensor data in Hz (samples per second). This is only required if X is not a sensor structure.
Before computing the spectrum level at each frequency of a signal, dsf() performs a low-pass filter on the signal. This prevents high-frequency transients, e.g., in foraging, from dominating the spectra. The dsf() function allows the user to specify the cut-off frequency in Hz of this low-pass filter (fc). If not specified by the user, fc is set to 5 Hz by default. It is important to note that the interval length over which dsf() will be computed should be at least Nfft/fs seconds, i.e., 20 s for the default FFT length.
Deep-Dive on Functionality
For a better understanding of how this works, let us plot the power spectra of the accelerometer and magnetometer data in the animal frame to visualize the dominant stroke cycle frequency. This will also help us to visualize not only the dominant stroke cycle frequency but how this differs from the low-frequency postural changes in orientation and how we can set an appropriate cut-off frequency for the complementary filter comp_filt(), to appropriately separate stroking, flapping or striding from postural signals. We will keep using the example of the beaked whale stroking during the ascent phase of a deep dive from the dsf tutorial.
bw_file_path = xxx
bw = load_nc(bw_file_path) ;
Aseg = crop_to(bw.A,[35*60 42*60]) ;
Mseg = crop_to(bw.M,[35*60 42*60]) ;
fs = Aseg.sampling_rate ;
Nfft = round(20*fs) ; % default dsf FFT length
% force Nfft to the nearest power of 2
Nfft = 2^round(log(Nfft)/log(2)) ;
[Sa,fa] = spectrum_level(Aseg.data,Nfft,fs,Nfft,floor(Nfft/2)) ;
[Sm,fm] = spectrum_level(Mseg.data,Nfft,fs,Nfft,floor(Nfft/2)) ;
bw_file_path <- system.file("extdata", "testset1.nc", package = 'tagtools', mustWork = TRUE)
bw <- load_nc(bw_file_path)
Aseg <- crop_to(X = bw$A, tcues = c(35*60, 42*60))
Mseg <- crop_to(X = bw$M, tcues = c(35*60, 42*60))
fs <- Aseg$sampling_rate
Nfft <- round(20*fs); # default dsf FFT length force Nfft to the nearest power of 2
Nfft <- 2^round(log(Nfft)/log(2))
Sa <- spectrum_level(x=Aseg$data,nfft=Nfft,sampling_rate=fs,w= Nfft,nov=floor(Nfft/2))$SL
fa <-spectrum_level(x=Aseg$data,nfft=Nfft,sampling_rate=fs,w= Nfft,nov=floor(Nfft/2))$f
Sm <-spectrum_level(x=Mseg$data,nfft=Nfft,sampling_rate=fs,w=Nfft,nov=floor(Nfft/2))$SL
fm <- spectrum_level(x=Mseg$data,nfft=Nfft,sampling_rate=fs,w=Nfft,nov=floor(Nfft/2))$f
#[Sa,fa] <- spectrum_level(x=Aseg$data,nfft=Nfft,sampling_rate=fs,w= Nfft,nov=floor(Nfft/2));
#[Sm,fm] <- spectrum_level(x=Mseg$data,nfft=Nfft,sampling_rate=fs,w=Nfft,nov=floor(Nfft/2));
Here, we plot a spectrogram with 512 length FFTs (because the sampling frequency fs is 25 Hz). spectrum_level() divides the data up into as many blocks of 512 samples as it can, computes the FFT of each of these, and averages the power of the results. It returns the spectral average in decibels i.e., 20*log10(signal size), which is good for sound that has a high dynamic range, but is not always so useful for movement data. We can reverse the decibels by computing 10.^(Sa/20) in MATLAB (note that the operator .^ means compute 10 to the power of each element in Sa/20 and return a vector of the same size as Sa). The units of the spectrum will then be in m/s2/√Hz for acceleration and µT/√Hz for magnetometers.
What is the frequency and time resolution of this spectrogram? Hint: one is 512 divided by the sampling rate and the other is the inverse
The time resolution is 512/25 ≈ 20 seconds, and the frequency resolution is 25/512 ≈ 0.05 Hz.
You can plot the spectral averages for Aseg and Mseg on the same figure:
figure ;
ax1 = subplot(2,1,1) ;
plot(fa,10.^(Sa/20)) ;
xlabel('Frequency, Hz') ;
ylabel('Spectrum level, m/s2') ;
ax2 = subplot(2,1,2) ;
plot(fm,10.^(Sm/20)) ;
xlabel('Frequency, Hz') ;
ylabel('Spectrum level, µT') ;
legend({'x-axis','y-axis','z-axis'}) ;
linkaxes([ax1 ax2],'x') ;
plot(x = fa, y = (10.^(Sa/20)), xlab = "Frequency, Hz", ylab = "Spectrum level, m/s2")
plot(x = fm, y = (10.^(Sm/20)), xlab = "Frequency, Hz", ylab = "Spectrum level, m/s2")
Sa <- data.frame(Sa, fa = fa)
names(Sa) <- c('x-axis','y-axis','z-axis', 'fa')
Sa <- tidyr::pivot_longer(Sa,
cols = c('x-axis','y-axis','z-axis'),
names_to = 'axis',
values_to = 'spec_lev')
Sm <- data.frame(Sm, fm = fm)
names(Sm) <- c('x-axis','y-axis','z-axis', 'fm')
Sm <- tidyr::pivot_longer(Sm,
cols = c('x-axis','y-axis','z-axis'),
names_to = 'axis',
values_to = 'spec_lev')
# make first plot (don't display yet)
ax1 <- gf_path(10^(spec_lev/20) ~ fa,
color = ~axis,
data = Sa,
xlab = 'Frequency, Hz',
ylab = 'Spectrum level, m/s2')
# make second plot (don't display yet)
ax2 <- gf_path(10^(spec_lev/20) ~ fm,
data = Sm,
color = ~axis,
xlab = 'Frequency, Hz',
ylab = expression(paste('Spectrum level, ', mu, 'T', sep = '')))
# note: the interactive version won't work with the mu (special character) so alternative
aax2 <- gf_path(10^(spec_lev/20) ~ fm,
data = Sm,
color = ~axis,
xlab = 'Frequency, Hz',
ylab = 'Spectrum level, uT')
# display both together
grid.arrange(ax1, ax2, ncol = 1)
INTERACTIVE VERSION
iax1 <- ggplotly(ax1)
iax2 <- ggplotly(aax2)
plotly::subplot(iax1, iax2,
nrows = 2,
titleX = TRUE, titleY = TRUE,
shareX = TRUE)
Expand to show figure …
Coming soon …
You could zoom in to get a closer look at the figures.
In this plot, the longitudinal axis is plotted in blue, the dorso-ventral axis in yellow and the lateral axis in red. How do the spectra of Aseg and Mseg compare? (note: ignore any frequencies of the first peak: these are from slow changes in orientation). Looking first at the axis of each sensor with the strongest fundamental signal: is it the same axis for both sensors? If so, why?
Yes, the longitudinal axis. These data are from a swimming whale. The caudal propulsion of whales generates an oscillatory pitching rotation which shows up most clearly in the longitudinal axes of the accelerometer and magnetometer (remember that both the accelerometer and magnetometer are sensitive to orientation changes). Propulsion also produces specific acceleration (i.e., acceleration due to actual movement as opposed to gravity) and this shows up most clearly in the dorso-ventral axis in whale swimming. Because only the accelerometer is sensitive to specific acceleration, the dorso-ventral axis shows a clear peak in the top plot, albeit somewhat lower than in the longitudinal axis.
Do you see a 2nd harmonic (i.e., a smaller peak at twice the DSF) in either sensor? If not both, why do you think this might be? What spectral components (if any) are in the other axes?
There is a small second harmonic peak in the accelerometer data but not in the magnetometer data. In the accelerometer it is most apparent in the longitudinal axis.
Swimming animals that use caudal propulsion can produce one or two bursts of thrust per tailbeat cycle, e.g., by using a powerful down-stroke followed by a lazy up-stroke, or by using powerful strokes in both directions. The resulting thrust is forward-directed irrespective of whether it is produced once or twice per tailbeat.
The spectrum of the accelerometer longitudinal axis therefore tells us whether the whale is swimming with a lazy half-stroke or with two powerful half-strokes: the presence of some second harmonic indicates that two bursts of thrust are being made per tailbeat (although these may not be equally strong). In this example, the magnitude of the magnetometer data in the second harmonic is almost null. This is because the magnetometer measures orientation changes but not specific acceleration such as thrust or heave. The body rotations in caudal propulsion that are measured by a magnetometer are close to sinusoidal and so only produce a first harmonic (or fundamental frequency).
We could also plot the sum of the spectral power in the three axes for the accelerometer and magnetometer data:
va = sum(10.^(Sa/10),2) ; % sum spectral power in the three axes
vm = sum(10.^(Sm/10),2) ; % sum spectral power in the three axes
figure ;
ax1 = subplot(2,1,1) ;
plot(fa,va) ;
xlabel('Frequency, Hz') ;
ylabel('Spectrum level, m/s2') ;
ax2 = subplot(2,1,2) ;
plot(fm,vm) ;
xlabel('Frequency, Hz') ;
ylabel('Spectrum level, µT') ;
linkaxes([ax1 ax2],'x') ;
va <- sum(10.^(Sa/10),2); #sum spectral power in the three axes
vm <- sum(10.^(Sm/10),2); # sum spectral power in the three axes
plot(x = fa, y = va, xlab = "Frequency, Hz", ylab = "Spectrum level, m/s2")
plot(x = fm, y = vm, xlab = "Frequency, Hz", ylab = "Spectrum level, m/s2")
Expand to show figure …
Coming soon …
You could zoom in to get a closer look at the figure. This output will be the same as when we run the dsf() function.
Verification
If the sampling rate of the sensor(s) used for computing the DSF is well above twice the highest stroke frequency of the animal, and the other assumptions are met, dsf will be well estimated. To verify that you have run dsf()correctly, you could plot the accelerometer or magnetometer signals at a time when the animal is locomoting. Look for complete oscillation cycles and count the seconds that it takes for the animal to perform a number of full cycles, for example 10 strokes/strides. The DSF should then be reasonably close to 10/t where t is your second count. This will give you an idea if the output of dsf() is correct, but it may be that you do not get the exact DSF as there could be periods when the animal locomotes at a higher or lower rate.
Plot Aseg, zoom, for example, in the period from 70 to 110 seconds. Using the function ginput(2) mark the start and end of the period containing 10 full strokes. Then calculate how many strokes the animal performs in a second.
figure ;
plott(Aseg.data(:,1),fs) ;
ylabel('Acc x-axis') ;
xlabel('Time, seconds') ;
xlim([70 100]) ;
section = ginput(2) ;
dsf_manual = 10/(section(2)-section(1)) ;
figure
plot(x = Aseg$data[,1], y = fs, xlab = "Acc x-axis", xlim=c(70,100))
section <- locator(n = 2) # section$x will be the 2 x values and section$y will be the two y values
dsf_manual=1/(diff(section$x))
Expand to show figure …
Coming soon …
Caveats
DSF is often broadly similar across animals with similar size and locomotion style. However, there can also be significant variation within and across species, e.g., some individuals tend to adopt faster or slower locomotion rates. In addition it is quite common that diving animals have a different dsf during the descent than the ascent phase, as an animal that encounters large forces (due to their tissue weight or buoyancy) may respond by stroking faster. The same may be true of wingbeat rates in birds with faster stroking needed on ascent than during horizontal or descending flight.
It is essential that the sampling rate of the sensor(s) used for computing the DSF is well above twice the highest stroke frequency of the animal, otherwise aliasing can occur. Signals alias when the sampling rate is less than twice the highest oscillation frequency in the signal. An aliased signal may still give an apparent DSF peak but the peak frequency will be incorrect. In addition, for the dsf() function to work well and give an obvious result you should aim for estimating DSF during periods of locomotion and avoid having long periods where the animal is performing other activities, such as grooming or making social displays, without actually stroking, flapping or walking. Also you may need to avoid periods when the animal is maneuvering rapidly, e.g., during foraging, as this can mask the spectral peak due to locomotion.
References
Miller, P. J. O., Johnson, M. P., Tyack, P. L., and 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
Sato, K., Watanuki, Y., Takahashi, A., Miller, P. J. ., Tanaka, H., Kawabe, R., Ponganis, P. J., Handrich, Y., Akamatsu, T., Watanabe, Y., Mitani, Y., Costa, D. P., Bost, C.-A., Aoki, K., Amano, M., Trathan, P., Shapiro, A., & Naito, Y. (2007). Stroke frequency, but not swimming speed, is related to body size in free-ranging seabirds, pinnipeds and cetaceans. Proceedings of the Royal Society B: Biological Sciences, 274(1609), 471–477. https://doi.org/10.1098/rspb.2006.0005
MacAyeal, L.C., Riskin, D.K., Swartz, S.M., Breuer, K.S. (2011). Climbing flight performance and load carrying in lesser dog-faced fruit bats (Cynopterus brachyotis). Journal of Experimental Biology,214 (5): 786–793. https://doi.org/10.1242/jeb.050195