Applying a Bessel Filter to Acceleration Data

Dear Matlab Community,

I am trying to learn how to apply a bessel filter to force data from a vertical jump, as the force data is similar to an inverted square wave when airborne. After reading the matlab pages on the besself commands, I cannot wrap my head around how to apply the filter to the data, similar to how I would with a butterworth filter. I am hoping that you all can help me figure this part out, or explain the inaccuracies of my understanding (if they exist) so that i can get this sorted.

Here's what my code looks like for a typical lowpass Butterworth filter:

sf = 1000; % sampling rate

cf = 50/(sf/2); % normalize cutoff frequency for low pass filter, 50 Hz cutoff

[b,a] = butter(4,cf,'low'); % determine filter coefficients (obtain coefficients)

Fz = filtfilt(b,a,Fz); % overwrite Fz raw data with filtered data

Here's what I have thus far for the Bessel filter:

sf = sampling rate

cf = 50(sf/2); % normalize cutoff frequency for low pass filter, 50 Hz cutoff

[b,a] = besself(4,cf,'low'); % determine filter coefficients (obtain coefficients)

What is minning now is how to apply the filter to the Fz data. I cannot locate any examples for this... All I can find are examples of plotting the frequency response.

Thanks for your assistance!

Answer 1

[Chunru]

Do you think you really need bessel filter (for linear group delay)? Otherwise, other filter types give better attenuation/ripple performance.

Unlike butter and others, the besself function does not support the design of digital Bessel filters and you have to do the conversion yourself from analog to digital (which than destroy the constant group delay property).

fs = 1024;                      % sampling rate
fc = 50;                        % cutt off frequency
% Bessel filter is analog filter
% Specify cut-off as rad/second
w0 = 2*pi*fc;
[z,p,k] = besself(6, w0,'low'); % zero, pole and gain form

% Convert to digital fileter
[zd,pd,kd] = bilinear(z,p,k,fs);   % z-domain zero/pole/gain
[sos,g] = zp2sos(zd,pd,kd);        % convert to second order section 

figure;
%freqz(b, a, 2048, fs);
freqz(sos, 2048, fs);

% filter signal
x = randn(5000,1);
y = filtfilt(sos, g, x);

plot(x); hold on;
plot(y);

[Daniel]

Thank you so much. Perfet answer.