DC component in accelerometer data - filter before or after XYZ Euclidean sum

Question

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.

Answer 1

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.