Filtration Process for human step detection using inertial measuring unit

Currently i am working on pedestrian step detection using inertial measuring unit (accelerometer), i need to do filtraton of my signal in preprocessing level. could anyone suggest which one will be the best filtration algorithm to get good accuracy.

here i have attched the normalized signal. looking for the kind response. Thanks

[Mathieu NOE]

hello

maybe you should first look at what frequencies must be kept / filtered out

this little code will do fft (with averaging if you specify it) and time / frequency analysis (spectrogram in short)

then , in a second step , you have to think which filter you need to apply and at which sections of your signal

[Muhammad Hammad Malik]

thanks for your response. I am trying to learn it. Could to just tell e how i can decide which frequencies i shoud filtered out. is there any method or just a random selection. Thanks

[Mathieu NOE]

maybe I should look first at your data , if you can share it

tx

[Mathieu NOE]

problem solved ?

[Muhammad Hammad Malik]

like the one you mentioned in your comment, i also want to do like that. want to apply filter using window size,i used fir filter but still could not get the exact result.

see attached .mat file. i wnt to design a filter with window size 2, and max cuttoff freq of 5 hz

i am also trying to do it in python, if you know about that

[Mathieu NOE]

hello

what kind of fitering (with fir) do you want to implement ? what is the target ?

FYI, below a code to read the accel data (here I picked the 3rd accel data) and do averaged fft , spectrogram , filtering and compare spectrum before and after filtering

hope it helps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load('IMU_Hand.mat');
% IMU_ULISS =   struct with fields:
%                     Mag: [3×36588 double]
%                    Gyro: [3×36588 double]
%                     Acc: [3×36588 double]
%                    time: [1×36588 double]
%                    size: 36588
%                step_idx: [1×220 double]
%     step_instant_target: [1×36588 double]

time = IMU_ULISS.time;
accel = IMU_ULISS.Acc;
signal = accel(3,:);
samples = length(signal);
dt = (time(end)-time(1))/(samples-1);
Fs = 1/dt;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1024;    % 
OVERLAP = 0.75;

% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 

% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;

%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 4;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);

% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end


figure(1),plot(freq,sensor_spectrum_dB,'b');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  

% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(fsg);
    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
    my_title = ('Spectrogram (dB (A))');
else
    my_title = ('Spectrogram (dB (L))');
end

% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');

%%%%%%%%%%%%%%%%%%%%%%
%% bandpass filter section %%%%%%
f_low = 0.25;
f_high = 5;
N = 2;
% signal_filtered = zeros(size(signal));
[b,a] = butter(N,2/Fs*[f_low f_high]);
% now let's filter the signal
signal_filtered = filter(b,a,signal);
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq,fft_spectrum_filtered] = myfft_peak(signal_filtered, Fs,NFFT,OVERLAP);
sensor_spectrum_filtered_dB = 20*log10(fft_spectrum_filtered);% convert to dB scale (ref = 1)

figure(3),semilogx(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before filter','after filter');
xlabel('Frequency (Hz)');ylabel(' dB')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end

function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);

samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    s_tmp((1:samples)) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);
signal = signal(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
    fft_spectrum = fft_spectrum(select);
    freq_vector = (select - 1)*Fs/nfft;
end 

[Muhammad Hammad Malik]

thanks for your code, i will look into it.

i want to detect features as many as possible from hand data to do step detection. i have extracted some features and now i want to train my algorithm to detect automatically my step instants using unsupervised ML. could you guide bit about this.

[Mathieu NOE]

hello

hve you already coded something ?

[Muhammad Hammad Malik]

yes, i have something in python but now i am moving to matlab

[Muhammad Hammad Malik]

i tried your but in result i am not getting complete signal instead a very short signal. my sampling freq is 200, and i want to use it after normalizing the data. see attached images of your filter, i couldnot understand why getting just this one, i need whole signal, i have also added normalized signal.kindly check.

norm

yours filter result

below is what i have done in python but result is not good.

""""""

from scipy.signal import kaiserord, lfilter, firwin, freqz

sample_rate = 200.0

# The Nyquist rate of the signal.

nyq_rate = sample_rate / 2.0

order=2

width = 1.2/nyq_rate

ripple_db = 10.0

N, beta = kaiserord(ripple_db, width)

cutoff_hz = 0.03

taps = firwin(N, cutoff_hz/nyq_rate, order, window=('kaiser', beta))

d = lfilter(taps, 1.0, n)

plt.plot(d, label='Filtered signal')

plt.xlabel('time (seconds)')

plt.grid(True)

plt.axis('tight')

plt.legend(loc='upper left')

plt.show()

"""""

[Mathieu NOE]

hello

again, before jumping on the code itself - what is your primary goal on this data (you have accel , gyro and mag infos)

which ones are of interest and what are you trying to do with the filtering ? once we understand the targets it's easier to come up with a code

tx

[Muhammad Hammad Malik]

i want to work both with accel and gyro. i my goal is to do step detection. for that first i need to normalize my accx,y,z and then apply filtration using window size and cut off frequency to smooth the signal and then will calculate features from this filter. This is the first thing at the moment i want to do

[Mathieu NOE]

hello

so basically you are looking at time events that would show a "bump" corresponding to a step ?

[Muhammad Hammad Malik]

yes. so for that first i need to do preprocessing. i used L2 normalized for normaliation of the data and now i want to apply filter to smooth the data, after that will extract features from that for further step detection.

Look at above picture, this is what i want. this filtration done on normalized acceleration of time series data.

this is done with 3hz cutoff freq but when i change 0.03hz cutoff i am still getting the same signal. i used lowpass FIR filter.