Filtration Process for human step detection using inertial measuring unit

Question

Muhammad Hammad Malik 2020-12-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/699260-filtration-process-for-human-step-detection-using-inertial-measuring-unit

编辑： Muhammad Hammad Malik 2021-1-15

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

14 个评论
显示 12更早的评论隐藏 12更早的评论

Mathieu NOE 2020-12-27

在 MATLAB Online 中打开

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 

Mathieu NOE 2021-1-13

hello

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

Muhammad Hammad Malik 2021-1-14

编辑：Muhammad Hammad Malik 2021-1-15

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.

请先登录，再进行评论。

请先登录，再回答此问题。

Filtration Process for human step detection using inertial measuring unit

14 个评论
显示 12更早的评论隐藏 12更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Filtration Process for human step detection using inertial measuring unit

14 个评论 显示 12更早的评论隐藏 12更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

14 个评论
显示 12更早的评论隐藏 12更早的评论