How to apply Butterworth filter to an inertial signal metrics?

Hello everyone, I have a CSV file of an inertial signal having 9 columns (three columns each for accelerometer, magnetometer, and gyroscope) and 5000 rows. I want to apply butterworth (BW) filter to clean this file. As per MATLAB documents, we can apply BW filter to array. Please guide me how to filter my inertial signal metrics. Should I need to break my metrics into array or any alternative approac is required to this filtering?

Answer 1

[Mathieu NOE]

hello

this is a demo code (with 3 channels dummy data) that will display time and averaged fft plot

it includes banpass and notch filters (I commented the last one as it was not your request but you can keep it under your arm for more future fun)

simply add one line of code to load your csv file (with csvread , readmatrix , readtable ...and so forth) and you're done !

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% dummy signals
samples  = 5000;
signal = randn(samples,3);
Fs = 1000;
dt = 1/Fs;
time = dt*(0:samples -1)';

%% band pass  filter section %%%%%%
%bandpass filter[Fc1=0.4 Hz   Fc2=45 Hz]
f_low = 40;
f_high = 245;
N_bpf = 2;           % filter order
[b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
signal_filtered = filtfilt(b,a,signal);

% %% notch filter 
% fc_notch = 50;
% w0 = 2*pi*fc_notch/Fs;
% p = 0.95; % something slightly less than 1 (will change the width and the depth of the notch filter).
% 
% % digital notch (IIR)
% num1z=[1 -2*cos(w0) 1];
% den1z=[1 -2*p*cos(w0) p^2];
% 
% % now let's filter the signal 
% signal_filtered = filtfilt(num1z,den1z,signal_filtered);
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 500;    % buffer length
OVERLAP = 0.5; % (50% ) overlap will be used only if NFFT is shorter than signal length

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
subplot(2,1,1),plot(time,signal);grid on
title(['Raw data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend('accelerometer', 'magnetometer', ' gyroscope')
subplot(2,1,2),plot(time,signal_filtered);grid on
title(['Filtered data Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend('accelerometer', 'magnetometer', ' gyroscope')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);
df = freq(2)-freq(1); % frequency resolution 

figure(2),
subplot(2,1,1),plot(freq,20*log10(spectrum_raw));grid on
title(['Raw data Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
legend('accelerometer', 'magnetometer', ' gyroscope')
subplot(2,1,2),plot(freq,20*log10(spectrum_filtered));grid on
title(['Filtered data Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
legend('accelerometer', 'magnetometer', ' gyroscope')

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 percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

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

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

%    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*ones(1,channels));
        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