Remove High Frequency Noise Due to Tag Vibration

Usage

In this tutorial, you will learn how to identify and remove noise in the accelerometer and magnetometer data due to a tag partially attached to the animal, while keeping the biological information: animal orientation (low frequency) and strikes/flinches (high frequency). This tutorial is specific for tags attached to the animal with suction cups or that are attached with a harness, not for glued tags. Sometimes it could happen that the tag is not properly attached (only two suction cups attached) or that two suction cups get detached earlier than the other two and so the tag vibrates. These vibrations can be identified as high frequency noise in a spectrogram and if the frequency of vibration does not coincide with any biological movement this can be removed by applying a low-pass filter.

How to identify a high frequency noise due to a vibrating tag

Lets plot the depth along with the acceleration data. But first as always lets upload that data.

xxx_file_path = 
xxx = load_nc(xxx_file_path)

plott(xxx.P,xxx.Aa)
subplot(211)
ylabel('Depth')
subplot(212)
ylabel('Acc')
legend({'x-axis','y-axis','z-axis'})

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

plott(X = list(Depth = xxx$P, Acc = xxx$Aa))

In this plot, in MATLAB the longitudinal axes is plotted in blue, the dorso-ventral axes in yellow, and the lateral axes in red. While in R, the longitudinal axes is plotted in blue, the dorso-ventral axes in green, and the lateral axes in red.

You could check if this spiky sine waves are throughout all the deployment or only at the end of the deployment. If this is the case, then we recommend that you plot two power spectra as specified below, one for a period of smooth sine waves and another for a period with spiky sine waves, so you could easily see the vibration frequency and where to set the cut off filter frequency of the low-pass filter to remove the noise due to tag vibration.

Low-pass filter design

To know at which frequency this vibration occurs you would need to plot the power spectra of the accelerometer data. This will also help us to visualize not only the tag vibration frequency, but how this one does or does not differ from the dominant stroke cycle frequency and the low frequency postural changes in orientation and how we can set an appropriate cut off-frequency for the low-pass filter, to appropriately remove the noise signal due to the vibration of the tag and keep the stroking, flapping or striding as well as the postural signals.

Select the whole deployment or a segment of the deployment to plot the power spectra. To do this use the function crop_to():

fs = 50; % Sampling frequency (Hz)
% x1 and x2 are the starting and end point of the accelerometer segment to be plotted in seconds.
Aseg = crop_to(Aa, [x1*fs, x2*fs]);

fs <- 50 # Sampling frequency (Hz)
# x1 and x2 are the starting and end point of the accelerometer segment to be plotted in seconds.
Aseg <- crop_to(Aa, sampling_rate = fs, tcues <- c(x1, x2))

Here, we plot a spectrogram with 1024 length Fast Fourier Transforms (FFTs) (because the sampling frequency fs is 50 Hz). spectrum_level() divides the data up into as many blocks of 1024 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 which 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/s²/√Hz for acceleration and µT/√Hz for magnetometers. For a greater detail see the technical explanation for dsf.

Lets plot the spectral averages for Aseg:

Nfft = 1024 % force the length of the Fast Fourier Transform to be 1024 samples 
fs = Aseg.sampling_rate 
[Sa,fa] = spectrum_level(Aseg.data,Nfft,fs,Nfft,floor(Nfft/2));

figure
plot(fa,10.^(Sa/20))
xlabel('Frequency, Hz')
ylabel('Spectrum level, m/s2')
legend({'x-axis','y-axis','z-axis'})

Nfft <- 1024 # force the length of the Fast Fourier Transform to be 1024 samples 
fs <- Aseg$sampling_rate 


[Sa,fa] <- spectrum_level(x=Aseg$data,nfft=Nfft,sampling_rate=fs,w= Nfft,nov=floor(Nfft/2));

plot(x = fa, y = (10.^(Sa/20)), xlab = "Frequency, Hz", ylab = "Spectrum level, m/s2")

In this plot, in Matlab the longitudinal axes is plotted in blue, the dorso-ventral axes in yellow and the lateral axes in red. While in R, the longitudinal axes is plotted in blue, the dorso-ventral axes in green, and the lateral axes in red.

Here, you would see two, three or four peaks. The first peak with the lowest frequency is due to the animal orientation, the second peak is placed at the dominant stroke cycle frequency of the whale. In swimming animals that use caudal propulsion, only one or two bursts of thrust can be produced per cycle (i.e., corresponding to each half-stroke) and so the kinematic waveforms can be characterized by their fundamental (the dominant stroke cycle frequency) and second harmonic (which occurs at twice the stroking rate). It could also be visible another peak at higher frequencies which may be related with the noise due to tag vibration. The cut-off filter to be used in the complementary filter should be set before the start of this last peak, so we remove the whole upper-frequency region of the signal that results from the tag vibrating.

Apply the low-pass filter to the accelerometer and magnetometer data.

The cut off filter to be used in the complementary filter should be set before the start of this last peak, so as we remove the whole signal due to the tag vibrating. Call this value fc. Run a complementary filter on A to separate the biological from the noise data.

fc = 'YourValueHere'% your value for fc in Hz, a number, without quotes
Af = comp_filt(Aa, fc)

fc <- 'YourValueHere' # your value for fc in Hz, a number, without quotes
Af <- comp_filt(beb$Aa, fc = fc) 

The complementary filter returns a structure containing two data matrices: the low-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 two accelerometers in the tag—one is only sensitive to low frequencies (which in this case are related to the biological data of interest) and the other is only sensitive to high frequencies (which in this case refer to the noise signal due to the tag vibration). Here we are only interested in the biological signal so we will ignore the high frequency signal.

Abio = Af{1} % low frequency A data

Abio <- Af[[1]]     # low 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 two filtered accelerations along with the dive profile:

fs = A$sampling_rate 

plott(P.data,fs,Aa.data,fs,Abio.data,fs)
axis ij
legend({'x-axis','y-axis','z-axis'})

fs <- A$sampling_rate 

plott(X = list(`Depth (m)` = xxx$P$data,
               `Accel` = xxx$Aa, 
               `Acc bio` = Abio),
      fsx = fs)

Compare A with Abio, after removing the noise due to the tag vibration, the periods of fluking strokes in Abio should be seen as a smooth sine wave signal, while in A this sine wave will be spiky.

Apply the same filter to the magnetometer data, plot the data and compare it with the unfiltered magnetometer data.

Mf = comp_filt(xx.Ma, fc)
Mbio = Mf{1} % low frequency M data

plott(xx.P.data,fs,xx.Ma.data,fs,xx.Mbio.data,fs)
axis ij
legend({'x-axis','y-axis','z-axis'})

Mf <- comp_filt(xx$Ma, fc = fc) 
Mbio <- Mf[[1]]  # low frequency M data

plott(X = list(`Depth (m)` = xxx$P$data,
               `Magnet` = xx$Ma, 
               `Mag bio` = Mbio),
      fsx = fs)

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.