1
$\begingroup$

I'm working on filtering accelerometer XYZ data. The aim is gait-pattern recognition via finding characteristic frequency points and feeding these to a artificial neural network. The trouble is, and commonly seems to be so, filtering out a massive DC component (a strong peak at 0 Hz) in the power-vs-frequency domain, which overshadows higher frequency components. I'm working with Octave.

In literature, the typical approach is to gather accelerometer XYZ data, smooth some noise jitter with a sliding-window-average, apply a 4th order Butterworth highpass (at 0.9 Hz cutoff) and then take an Euclidean sum of the individual processed components (XYZ that is).

So far, whatever I try, there is always a strong peak centered at 0 Hz. Trying has included: 1) Subtracting mean of raw data (x=x-mean(x), for all channels) 2) Various orders of Butterworth on individual channels before the Euclidean sum, in various settings

What did get rid of the strong 0 Hz peak was applying the Butterworth filter on the Euclidean sum of XYZ channels. The characteristic 2 Hz gait peak remains and some higher order frequencies are amplified.

My question is: is there anything fundamentally wrong in filtering the XYZ data after the Euclidean sum? Or: is this correct? Or: what is happening, why is the 0 Hz peak not getting filtered with individual-channel filtering (that would be a hint at something non-distributive).

Here are some plots: Raw XYZ accelerometer data Raw XYZ accelerometer data http://www.hot.ee/jaaniussikesed/rawdataxyz.png

Data with very slight sliding window average (pretty much raw) Data with very slight sliding window average (pretty much raw) http://www.hot.ee/jaaniussikesed/walk_xyz_nf_fft.png

4th order Butterworth, cutoff at about 1 Hz, on individual XYZ before Euclidean sum 4th order Butterworth, cutoff at about 1 Hz, on individual XYZ before Euclidean sum http://www.hot.ee/jaaniussikesed/walk_xyz_f_fft.png

4th order Butterworth, cutoff at about 1 Hz, taken on a Euclidean sum of XYZ accelerometer data. Notice no 0 Hz peak. 4th order Butterworth, cutoff at about 1 Hz, taken on data after Euclidean sum of XYZ http://www.hot.ee/jaaniussikesed/walk_euc_f_fft.png

4th order Butterworth filter frequency response: 4th order Butterworth filter frequency response http://www.hot.ee/jaaniussikesed/frq_response.png

Some data-data: Samples: 1112 Sample rate: 84 samp/s

Butterworth implementation in Octave: Fc= 0.5; #Cutoff frequency in Hz Wcn=Fc/(srate/2); #Normalized cutoff frequency order=4; #Order of filter [B,A]=butter(order,Wcn,"high"); #Define Butterworth highpass filter

Note1: changing the cutoff frequency between 0.1-1 does not really change the picture. Note2: there have been quite many posts on a very similar problem, but they all have their nuances, as does this I believe.

Thank you for any help or/and advice!

EDIT1: Here is the code (running Octave 4.0.0, packages: control, data-smoothing, optim, signal, struct):

# Load data
load("ID:walk-2016-05-02T22:58:23+02:00.mat");

# Separate XYZ and clip
clipc=300;
ax=datamat(clipc:(end-clipc),1);
ay=datamat(clipc:(end-clipc),2);
az=datamat(clipc:(end-clipc),3);

srate=floor(srate); #sample rate

# Design filter
pkg load all
Fc= 1; #Cutoff frequency in Hz
Wcn=Fc/(srate/2); #Normalized cutoff frequency, srate (sample rate) loaded from file
order=4; #Order of filter
[B,A]=butter(order,Wcn,"high");

# Apply filter
axf=filtfilt(B,A,ax);
ayf=filtfilt(B,A,ay);
azf=filtfilt(B,A,az);

#Euclidean norm
afen=(axf.^2+ayf.^2+azf.^2).^(1/2);

#FFT
N_fft=2^13; #Frequency sample points, empirical
afft=fft(afen,N_fft); #Centered double sided FFT
freqpwr=afft.*conj(afft)/(N_fft*rows(afft)); #Power of each frequency spectrum component
freqpoints=srate*(0:N_fft/2-1)/N_fft; #Bins/points of FFT frequency axis

#Plot power spectral density
plot(freqpoints,freqpwr(1:N_fft/2));
title("One Sided Power Spectral Density");       
xlabel("Frequency (Hz)");
ylabel("Power spectral density");

fprintf("Program paused. Press enter to continue.\n");
pause;
close;

Link to data file.

$\endgroup$
1
$\begingroup$
  1. I notice that after the butterworth filtering the high frequency components are attenuated. That suggests that somehow you are actually low pass filtering. There might be a mistake in your code.

  2. A sliding window averaging is a crude low pass filter. Therefore, given that you have both a low pass filter (window averaging) and a high pass filter, you might try combining both operations by using a band pass filter.

  3. You have to remove the DC offset from each channel before taking the magnitude (euclidean sum). Consider a simple example - the offsets are $[1,2,3]$ and the actual acceleration vector is $[1,1,0]$. The combined magnitude is $\sqrt{(1+3)^2 + (1+2)^2 + (0+3)^2} = \sqrt{34}$. But if the values for two axes were flipped, say $[1,0,1]$, you would get $\sqrt{(1+3)^2 + (0+2)^2 + (3+1)^2} = \sqrt{36}$. Ideally the magnitude should be invariant to rotation of the axes.

$\endgroup$
  • $\begingroup$ +1: The non-DC peaks should definitely not be attenuated after the Butterworth filter. $\endgroup$ – Peter K. Apr 29 '16 at 11:24
  • $\begingroup$ I managed to pre-filter the data to remove the mentioned "gravity bias". But the filtering is still not working out for me. Could you perhaps have a look at the code? I pasted my full code and also added a link to the .mat datafile with the raw (g accel. removed) values. $\endgroup$ – Pythzepf May 2 '16 at 22:34
  • $\begingroup$ @Pythzepf The magnitude is alway positive so it will have a large DC component itself. You could try applying the high pass filter to the magnitude - afen=filtfilt(B,A,afen);. To plot the data try using the log log(freqpwr(1:N_fft/2)), as that will enable you to more easily see frequencies. $\endgroup$ – geometrikal May 2 '16 at 23:48
  • $\begingroup$ Are you only interested in changes in acceleration magnitude? $\endgroup$ – geometrikal May 2 '16 at 23:53
  • $\begingroup$ @geometrikal Thank you for the suggestions! Would there be anything fundamentally wrong with filtering the magnitude instead of the magnitude components (XYZ) separately? Am I losing or producing information? That kind of was my initial question, since filtering the magnitude gets (got) rid of the DC in all considerations. Log is an option, but findpeaks() will overproduce results on a log scale. My interest is in finding features for the artificial neural network. I plan on using both the magnitude and frequency of the significant peaks. $\endgroup$ – Pythzepf May 3 '16 at 11:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.