Issues with creating a frequency weighting filter for a iso-standard

Hello! I want to mention in the beginning that i am new to signal processing with MATLAB.

I am currently working on weighting a vibration simulation and some acceleration data. The standard I am using for the weighting gives the filter functions using the Laplace variable. The three functions that are used are a high pass & low pass filter, and lastly an a-v transition filter (pic included)

For this the band limiting i use butterworth filters, which i think should be equal to the first function. For the transfer function i use the bilinear tool to get filter coefficients using the info from the standard, and then i filter the signal using all three filters.

It is worth mentioning that the simulation i use gives a displacement value (time series) where i use FFT and thereafter differentiate twice in fourier domain to get a acceleration value. After that i weigh the signal with the weighting filter below. Am i doing this in the correct order? Lastly i convert to dB and plot and this is the result:

I am pretty sure that the filter fails at weighing correctly according to the standard. For example, the weighting factor that is multiplied with the acceleration at f = 100 Hz should be 0.160, and in the plot they are almost the same.

Here is the code:

function filtered_signal=Wh_weighting(acceleration_signal,Fs)
%Purpose: to weight an acceleration signal with Wh weighting, for hand arm
%vibrations

%Input: The time series acceleration signal and the sampling frequency of
%the signal

%Output: The filtered signal using Wh filtering according to ISO 8041 and
%Iso 5349

%Constants from SS_EN_ISO_8041_1_2017 page 12 table 3 Wh
f1=10^.8;
w1=2*pi*f1;
f2=10^3.1;
w2=2*pi*f2;
f3=100/2/pi;
f4=f3;
Q1=1/sqrt(2);
Q2=Q1;
Q4=0.64;
w3=2*pi*f3;
w4=2*pi*f4;

%The nyquist limit
nyqusit_limit=Fs/2;

%A second order highpass filter for band limiting of the signal, b1 and a1
%are coefficients of the highpass filter transfer function in z domain.
[b1,a1]=butter(2,f1/nyqusit_limit,'high');

%A second order lowpass filter for band limiting of the signal.  b2 and a2
%are coefficients of the lowpass filter transfer function in z domain.
[b2,a2]=butter(2,f2/nyqusit_limit,'low');%low pass filter

%Frequency weighting filter
%This filter uses the coeffcients described in ISO 5349-1 (commented is ISO 8041) equation (3) in S
%domain and transforms its coeeficients to z-domain using the bilinear
%transform. So A3 and B3 are s-domain coeffcients and a3 and b3 are
%z-domain coefficients.
%B3=[1/w3 1];
B3 = [2*pi*f4^2 f3*(2*pi*f4)^2];
%A3=[1/w4/w4 1/Q4/w4 1];
A3 = [f3 2*pi*f4*f3/Q4 f3*(2*pi*f4)^2];
[b3,a3]=bilinear(B3,A3,Fs);%a-v transition filter

%Here the signal is filtered several times by different filter to give an
%overall filtering.

%the input signal is first filtered with the low pass filter
filtered_signal=filter(b2,a2,acceleration_signal);

% the output of line 44 is then filtered again with the highpass filter
filtered_signal=filter(b1,a1,filtered_signal);

%lastly the signal is filtered with the frequency weighting filter
filtered_signal=filter(b3,a3,filtered_signal);
%filtered_signal=filter(b3,a3,acceleration_time_series_signal);

% rms(filtered_signal) gives the declaration value a_hw for a mono axial
% vibration


end

Am i implementing this totaly wrong maybe?

Thank you in advance for any help received :)

[Star Strider]

I cannot easily follow the code, or determine what relationship have to ‘Fs’.

To see if the filter Bode plots are what you expect them to be, use the freqz function —

f1=10^.8;
w1=2*pi*f1;
f2=10^3.1;
w2=2*pi*f2;
f3=100/2/pi;
f4=f3;
Q1=1/sqrt(2);
Q2=Q1;
Q4=0.64;
w3=2*pi*f3;
w4=2*pi*f4;

Fs = 65536;                                             % Using Plots For Reference, Assuming Nyquist Freq = 250

B3 = [2*pi*f4^2 f3*(2*pi*f4)^2];
%A3=[1/w4/w4 1/Q4/w4 1];
A3 = [f3 2*pi*f4*f3/Q4 f3*(2*pi*f4)^2];
[b3,a3]=bilinear(B3,A3,Fs);%a-v transition filter

figure
freqz(B3, A3, 2^16, Fs)

Do the same analysis for , and use the filtfilt function, not filter, to do the actual filtering.

I am not posting this as an nswer becausse I am not certain that I understand the question.

EDIT — (16 Apr 2023 at 21:06)

Corrected ‘Fs’, and re-ran code.

.

[Elias Jakobsson]

Thank you for your answer :)

I notice now that my post might be a bit confusing since it is alot of information at a time, sorry about that. I tried to keep it as simple as possible. I don't think there is any point though to upload the main code since it is very long and has alot of different parts. The only thing i do there related to this function is plot, and fft.

Fs is the sampling frequency and is set to Fs = 65536. I missed to mention that.

So essentially the question i have is if i am using the signal processeing functions correctly according to the filter functions given in the standard, like f.e is bilinear the correct way to create the a-v transition? As you might have noticed i am quite new to this, so i apologize in advance :/

The weighting factors that are to be multiplied with the input signal should look like this (from the standard):

Is this something i can validate with freqz?

Also, when i use filtfilt instead of filter, the result is quite odd:

The one thing i notice from the start is that the filtered signal should never be larger in value than the input signal if i am not mistaken.

I am not really sure what to do now, maybe i am approaching this in the wrong way?

[Star Strider]

The bliniear approach appears to be appropriate and the implementation correct, however I fail to understand the rest of it.

Most filters do not amplify (some active filters may amplify, although that is not common practice) so the filtered signals have magnitudes that are never greater than the filter passbands.

I find it interesting that what appear to be the orange magnitude functions end at different frequencies (perhaps 105 to 175 Hz) while the frequency of the blue signals go out to 250 Hz. That indicates to me that the filters are not created with the same sampling frequency as the signal (assuming that the frequency upper limits are the Nyquist frequencies of the various filters and signals).

I am not certain that your approach is wrong, simply that the filters and signals may be inconsistent in the sampling frequencies used to create the filters. The sampling frequencies supplied to the filters need to be consistent with the sampling frequencies of the signals they are to work with for the filters to operate correctly.

[Elias Jakobsson]

Hm that is interesting, i am pretty sure i use the same variable Fs for the entirety of my code, especially when running the ffts. There might be an issue on that aspect that i am missing. I can supply a snippet of the relevant code incase you can find something:

%....
% Simulation fft calculation, this is  a simple fft, I used my own fft code
[fft1,freq1]=FFTwelch(sim_solution(:,6,1),Fs); %sim_solution is the output from the ode-solver
% sim_solution(:,6,1) is a certain velocity (6) for a certain applied
% force (1).

% Here i differentiate once in fourier domain to get the acceleration:
fft_sim_db1 = mag2db(abs(fft1.*(2*pi*freq1))./1e-6);

% This line is for the filtering, differentiate first
diffed_sol1 = abs(fft1.*(2*pi*freq1));

% Use the weighting function
fft1_weighted= Wh_weighting(diffed_sol1,Fs);

% Convert to dB
fft_sim_weighted_db1 = mag2db(fft1_weighted./1e-6);


%%Plotting

axs=[0 250 0 150 ];
close all
figure
subplot(3,1,1)

plot(freq1,fft_sim_db1)
hold on
plot(freq1,fft_sim_weighted_db1)

title("Simulation")
ylabel("Vibration acceleration level [dB]")
xlabel("Frequency [Hz]")
legend("Unfiltered","Filtered")
axis(axs)

%The rest is for the measurement data, so it is left out

% subplot(3,1,2)
% plot(freqs(:,1), measurement_fft_db(:,1))
% hold on
% plot(freqs(:,1),measurement_fft_weighted_db(:,1))
% title("Measurement Operator 1")
%  axis(axs)
% ylabel("Vibration acceleration level [dB]")
% xlabel("Frequency [Hz]")
% legend("Unfiltered","Filtered")
% 
% subplot(3,1,3)
% plot(freqs(:,4),measurement_fft_db(:,4))
% hold on
% plot(freqs(:,4), measurement_fft_weighted_db(:,4))
% title("Measurement Operator 2")
% axis(axs)
% ylabel("Vibration acceleration level [dB]")
% xlabel("Frequency [Hz]")
% legend("Unfiltered","Filtered")

%% Function FFTwelch
% This is the function for the fft process

function [FFT_amplitude_spectrum,Frequencies]=FFTwelch(input_signal,Fs)
%Purpose: This function gives the Smallband FFT amplitude spectrum of a signal.
% Input: the time series signal and its sampling frequency.

%This is turns of a warning for an already filtered signal.
warning('off','signal:internal:filteringfcns:ForcedAllpassDesign');% avoid allpass filter warning

%Anti-aliasing: This line of code does a lowpass filter for frequencies
%above the Nyquist limit, Fs/2, in order to avoid aliasing.
filtered_signal=lowpass(input_signal,Fs/2,Fs,"Steepness",.99); %in order to avoid aliasing we filter above the nyquist freq, also to avoid windowing

%Takes the length of the signal.
N=length(filtered_signal);

%Determines the number of averages made in the pwelch() function. The more
%averages the better.
NFFT=N;

% Power spectral density: Pwelch gives the power spectral density of the
% signal. The input is the signal, a rectangular window, which does
% nothing. We use a rectwin so that pwelch does not assume any windows, and
% we do not need windows since we have done anti-aliasing filters.
[Power_spectral_density,Frequencies] = pwelch(filtered_signal,rectwin(N),[],NFFT,Fs);%power spectral density

% Calculating the FFT: The fft is the square root of the psd times 2.
FFT_amplitude_spectrum=sqrt(Power_spectral_density*2);
FFT_amplitude_spectrum(1)=FFT_amplitude_spectrum(1)/2;% we do not want to double the mean(zero amplitude)
FFT_amplitude_spectrum(end)=FFT_amplitude_spectrum(end)/2;
end

Sorry for not giving the entire code but i think it would be unneccessary for this particular problem. Also i should add that there are parts of this code that i did not write myself, hence the confusion from my side.

Thank you for your time btw :)

[Star Strider]

I do not see anything wrong with the posted code. The lowpass call to create ‘filtered_signal’ is appropriate, although I prefer 'ImpulseResponse','iir' to design an efficient elliptic filter.

My concern is that and must be supplied with the same sampling frequency as the signal they are being used to filter. That value of ‘Fs’ needs to be derived from the sampling frequency of each signal in order for the filters to work correctly with it. (If all the signals have the same sampling frequency, then one filter design will work for all of them.)

If the signal sampling frequency — the consecutive sampling inteervals — are not constant, use the resample function to make them constant. (This is required for most signal processing techniques, the only exception I am aware of being the nufft function.) Then use that sampling frequency to design the filters and use those filters to process the resampled signals.

.

[Elias Jakobsson]

Ah i see, that makes total sense! At which point and in which way do you think is appropriate to resample these signals? I attempted to do this directly after the fft function (the code shows only one force and one position, although i do the same process for several others):

[fft1,freq1]=FFTwelch(sim_solution(:,6,1),Fs);

[fft1re,freq1re] = resample(fft1,freq1,Fs); %Error here
%Used this method from "help": 
%[Y,Ty] = resample(X,Tx,Fs) uses interpolation and an anti-aliasing filter to resample the signal
% at a uniform sample rate, Fs, expressed in hertz (Hz).

% Differentiation and weighting
diffed_sol1 = abs(fft1re.*(2*pi*freq1re));
fft1_weighted= Wh_weighting(diffed_sol1,Fs);

And i got the errors:

Error using zeros
Requested 2147483648x1 (16.0GB) array exceeds maximum array size preference (15.9GB). This might cause MATLAB to become unresponsive.

Error in signal.internal.resample.uniformResample (line 81)
    y = coder.nullcopy(zeros(syTrue));

Error in resample>nonUniformResample (line 281)
[y, h] = signal.internal.resample.uniformResample(xGrid, isDimValSet, Dim, dimIn, varargin{1}, ...

Error in resample (line 214)
            nonUniformResample(isDimValSet,Dim, m, method1, dimIn, xIn, ...

Related documentation

[Star Strider]

The resample function works in the time domain, so you have to use the time domain vectors with it. I am not certain what those are here, however that needs to be done first.

[sr,tr] = resample(s,t,Fs); % ‘s’ & ‘t’ Are The Original Signal & Time Vectors, ‘sr’ & ‘tr’ Are The Resampled Versions

Then, do all the signal processing (fft, filtering, and anything else) on the resampled signals.

.

[Elias Jakobsson]

Ok so the good news is that the simulation looks a bit more accurate than before, but the filter is still looking kinda bad so i am not sure what it could be:

I did the exact same routine with the same weighting function, now with the resampling in the time domain before fft:

%Resampling
[vr_6, tr] = resample(sim_solution(:,6,1),tid,Fs);

%FFT
[fft1,freq1]=FFTwelch(vr_6,Fs);

% Differentiation
diffed_sol1 = abs(fft1.*(2*pi*freq1));

% Weighting
fft1_weighted= Wh_weighting(diffed_sol1,Fs);

%% Weighting func
function filtered_signal=Wh_weighting(acceleration_signal,Fs)
%Purpose: to weight an acceleration signal with Wh weighting, for hand arm
%vibrations

%Input: The time series acceleration signal and the sampling frequency of
%the signal

%Output: The filtered signal using Wh filtering according to ISO 8041 and
%Iso 5349

%Constants from SS_EN_ISO_8041_1_2017 page 12 table 3 Wh
f1=10^.8;
w1=2*pi*f1;
f2=10^3.1;
w2=2*pi*f2;
f3=100/2/pi;
f4=f3;
Q1=1/sqrt(2);
Q2=Q1;
Q4=0.64;
w3=2*pi*f3;
w4=2*pi*f4;

%The nyquist limit
nyquist_limit=Fs/2;

%A second order highpass filter for band limiting of the signal, b1 and a1
%are coefficients of the highpass filter transfer function in z domain.
[b1,a1]=butter(2,f1/nyquist_limit,'high');

%A second order lowpass filter for band limiting of the signal.  b2 and a2
%are coefficients of the lowpass filter transfer function in z domain.
[b2,a2]=butter(2,f2/nyquist_limit,'low');%low pass filter

%Frequency weighting filter
%This filter uses the coeffcients described in ISO 8041 equation (3) in S
%domain and transforms its coeeficients to z-domain using the bilinear
%transform. So A3 and B3 are s-domain coeffcients and a3 and b3 are
%z-domain coefficients.
%B3=[1/w3 1];
B3 = [2*pi*f4^2 f3*(2*pi*f4)^2];
%A3=[1/w4/w4 1/Q4/w4 1];
A3 = [f3 2*pi*f4*f3/Q4 f3*(2*pi*f4)^2];
[b3,a3]=bilinear(B3,A3,Fs);%a-v transition filter
%[numd,dend] = bilinear(num,den,fs) converts the s-domain transfer function 
% specified by numerator num and denominator den to a discrete equivalent.

%Here the signal is filtered several times by different filter to give an
%overall filtering.

%the input signal is first filtered with the low pass filter
filtered_signal=filtfilt(b2,a2,acceleration_signal);

% the output of line 44 is then filtered again with the highpass filter
filtered_signal=filtfilt(b1,a1,filtered_signal);

%lastly the signal is filtered with the frequency weighting filter
filtered_signal=filtfilt(b3,a3,filtered_signal);
%filtered_signal=filtfilt(b3,a3,acceleration_signal);

% rms(filtered_signal) gives the declaration value a_hw for a mono axial
% vibration


end

[Star Strider]

I am not certain what the plots are. It could help to plot the original signal fft, the filter Bode plots, and the frequency domain of the filtered signal. I still do not understand what you are doing with respect to designing the filters in the context of the signals.

[Elias Jakobsson]

Yeah maybe I explained it poorly, I'll try and take it from the beginning. I use an ode solver for my project to solve for different displacements, velocities e.t.c in a dynamic system. From this result i isolate a specific velocity and resample it and use the FFT func, and thereafter differentiating this velocity in the f domain to get the acceleration. Lastly I filter this signal with the weighting function.

The weighting function is an attempt att creating the filter that the standard (that I mentioned in the main post) gives and it seems to be formulated correctly i hope.

I forgot to clarify earlier that the first subplots that i sent contain 1. The simulation(ode) FFT (blue: unfiltered, orange: filtered) 2. &3. FFT of measurement data from an accelerometer.

The plot that I most recently sent contains only the simulation FFT (blue) and it filtered version (orange).

Also I'm not sure I fully understand what these Bode plots of the filter show, is it supposed to correspong to the weighting factors that are multiplied with the acceleration value?

I hope I didn't missunderstand your latest comment, english is not my first language. If you could specify any more info you need to help me, I'll supply it. Thanks :)

[Star Strider]

Your explanation helps. I still have a problem with the filtered signal plots, since in the filter passband the signal should be the same as the unfiltered signal. Only in the stopbands should it be different, and significantly attenuated. That does not appear to be the situation, if I understand the plots correctly.

The Bode plots produced by the freqz function show the magnitude and phase response of the filter at the plotted frequencies. They should be the same as the filter design parameters, essentially being a ‘quality control¹ for the filter design.

[Elias Jakobsson]

This might be a rather dumb question, but how do i know what the appropriate sample rate should be? I have tried different rates, while still having it consistent in the process (resampling etc.), and the filtered fft becomes significantly different while the unfiltered fft stays the same. (The Fs of my measurement is Fs=65536 so it should be the same)

F.e i tried Fs = 10000 and Fs = 100000, and the lower value gives alot more discontinuities in the filtered fft.

Not really sure how to proceed from this sadly..

Answer 1

[Star Strider]

I just saw your latest Comment now.

The sampling frequency can be anything that you determine to be appropriate for the system you are integrating. The highest frequency in the signal (and I do not know what that is) should be less than the Nyquist frequency, creating the sampling frequency to be at least 2.5 times the highest frequency in the signal, and higher if possible.

Since you are creating your data with an ode45 (or similar functions) integration of your system, the easiest way to create consistent sampling intervals (the inverse of the sampling frequency) is to define ‘tspan’ as a vector of more than two elements. For example to create a time vector from 0 to 1000 seconds (or whatever time span you are using) with a sampling frequency of 5000, something like this will work:

Fs = 5000; % Sampling Frequency (Hz)

Tlen = 1000; % Signal Length (sec)

tspan = linspace(0, Tlen*Fs, Tlen*Fs+1)/Fs; % ‘tspan’ Vector

The integration will still be adaptive, however the results will be interpolated and will only be reported at the times corresponding to the multiple-element ‘tspan’ vector. There will likely be no need to do any resampling.

Now that you have the signal and a specific sampling frequency for it, all the filters can be created with that sampling frequency. They should then produce the result you want. The passbands and stopbands can be defined with respect to the Nyquist frequency (and they all must be less than the Nyquist frequency) so for example using this sampling frequency (5000):

Fn = Fs/2; % Nyquist Frequency

BandpassFreqs = [100 250]/Fn; % Normalised Frequencies

I am posting this as an answer because it might actually be the solution to this problem.

.

[Elias Jakobsson]

I have been using this time vector for my project:

Fs=65536; %Sample rate
Ts=1/Fs; %Step
tid=0:Ts:20; %Time vector

It is essentially the same vector. The signal is created with ode45 with this vector:

% Ode 45 options,
options = odeset('RelTol',1e-5,'AbsTol',1e-6);


%for loop that calculates the solution x, x contains position, velocity and
%gass mass parameters for the system
for counter=1:3
    [t,x]=ode45(@(t,x) State_space_new(t,x,K_vector,C_vector,M_vector,f_hand(counter)),tid,x_total_IC,options);
    sim_solution(:,:,counter)=x;
end

Where the state space function calculates the solution to the diff. eq. y' = Ay + f (As mentioned earlier the solution contains displacements and velocities).

What I am curious about though, is what normalized cutoff frequencies i should use for the butterworth filters. Right now i am using the resonance frequencies from the standard to calculate this (in the weighting function):

f1=10^.8;
f2=10^3.1;

%The nyquist limit
nyquist_limit=Fs/2;

[b1,a1]=butter(2,f1/nyquist_limit,'high'); %Highpass
[b2,a2]=butter(2,f2/nyquist_limit,'low');%Lowpass

filtered_signal=filtfilt(b2,a2,acceleration_signal); %Acceleration_signal is the input
filtered_signal=filtfilt(b1,a1,filtered_signal);

Is this the right way to implement the standard? (ref pic below)

These constants are used for the filter functions in the first pic i uploaded.

[Star Strider]

The Control System Toolbox makes defining and realising these filters straightforward.

Returning too the original definitions —

I copied that to here for convenience.

Coding the filters are then straightforward —

format long

Fs=65536; %Sample rate
Fn = Fs/2;
f1=6.310/Fn;                                % Normalise The Frequency By The Nyquist Frequency
w1=2*pi*f1;
f2=1258.9/Fn;                               % Normalise The Frequency By The Nyquist Frequency
w2=2*pi*f2;
f3=15.915/Fn;                               % Normalise The Frequencies By The Nyquist Frequency
f4=f3;
Q1=1/sqrt(2);
Q2=Q1;
Q4=0.64;
w3=2*pi*f3;
w4=2*pi*f4;
K = 1;
Ts = 1/Fs;                                  % Sampling Interval

opts = bodeoptions;                         % For The 'bodeplot' Plots
opts.FreqUnits = 'Hz';

s = tf('s');
Hb = (s^2*pi^2*f2^2) / ((s^2 + 2*pi*f1*s/Q1 + 4*pi^2*f1^2)*(s^2 + 2*pi*f2*s/Q1 + 4*pi^2*f2^2))

Hb =
 
                         0.007284 s^2
  -----------------------------------------------------------
  0.5 s^4 + 0.1715 s^3 + 0.02943 s^2 + 5.01e-05 s + 4.265e-08
 
Continuous-time transfer function.

figure
bodeplot(Hb, opts)
grid

Hbz = c2d(Hb, Ts, 'tustin')

Hbz =
 
  8.479e-13 z^4 - 1.696e-12 z^2 + 8.479e-13
  -----------------------------------------
        z^4 - 4 z^3 + 6 z^2 - 4 z + 1
 
Sample time: 1.5259e-05 seconds
Discrete-time transfer function.

Hbzn = Hbz.Numerator

Hbzn = 1×1 cell array
    {[8.479326400652394e-13 0 -1.695865280130479e-12 0 8.479326400652394e-13]}

Hbzd = Hbz.Denominator

Hbzd = 1×1 cell array
    {[1 -3.999994764868243 5.999984294618431 -3.999984294632135 0.999994764881946]}

Hw = ((s + 2*pi*f3)*2*pi*K*f4^2) / ((s^2 + 2*pi*f4*s/Q2 + 4*pi^2*f4^2)*f3)

Hw =
 
          1.048e-06 s + 3.198e-09
  ---------------------------------------
  0.0003434 s^2 + 1.482e-06 s + 3.198e-09
 
Continuous-time transfer function.

figure
bodeplot(Hw, opts)
grid

Hwz = c2d(Hw, Ts, 'tustin')

Hwz =
 
  2.328e-08 z^2 + 1.084e-15 z - 2.328e-08
  ---------------------------------------
               z^2 - 2 z + 1
 
Sample time: 1.5259e-05 seconds
Discrete-time transfer function.

Hwzn = Hwz.Numerator

Hwzn = 1×1 cell array
    {[2.328234100700516e-08 1.084134789677823e-15 -2.328233992287037e-08]}

Hwzd = Hwz.Denominator

Hwzd = 1×1 cell array
    {[1 -1.999999934147595 0.999999934147597]}

If you have the Control System Toolbox, you can run this code offline. If you do not have it, and if you want to change anything in my code, you can run it here in a Comment, and then once the post is submitted, you can copy the numerator and denominator vectors from it (as you can here, after I post this). You can then delete that Comment if you want to. I did my best to proofread this, however there could still be errors that I did not catch, so please check it.

Use these numerator and denominator vectors inside the square brackets (copy the square brackets as well and ignore the curly braces) as the transfer function vectors for the Signal Processing Toolbox filters. No further processing should be necessary for them. (The 'tustin' method is the same as the 'bililnear' method.)

.

[Elias Jakobsson]

Impressive! Although, I am having trouble impletmenting this in matlab since it produces warnings every different attempt i try.

It does not matter in what way I use these vectors but i always get these warnings:

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  3.700753e-17.

Though i am not certain i am using this toolbox the right way:

Fs=65536; %Sample rate
Fn = Fs/2;
f1=6.310/Fn;                                % Normalise The Frequency By The Nyquist Frequency
w1=2*pi*f1;
f2=1258.9/Fn;                               % Normalise The Frequency By The Nyquist Frequency
w2=2*pi*f2;
f3=15.915/Fn;                               % Normalise The Frequencies By The Nyquist Frequency
f4=f3;
Q1=1/sqrt(2);
Q2=Q1;
Q4=0.64;
w3=2*pi*f3;
w4=2*pi*f4;
K = 1;
Ts = 1/Fs;   

s = tf('s');
Hb = (s^2*pi^2*f2^2) / ((s^2 + 2*pi*f1*s/Q1 + 4*pi^2*f1^2)*(s^2 + 2*pi*f2*s/Q1 + 4*pi^2*f2^2));

Hbz = c2d(Hb, Ts, 'tustin');
Hbzn = Hbz.Numerator;
Hbznv = cell2mat(Hbzn);
Hbzd = Hbz.Denominator;
Hbzdv = cell2mat(Hbzd);

Hw = ((s + 2*pi*f3)*2*pi*K*f4^2) / ((s^2 + 2*pi*f4*s/Q2 + 4*pi^2*f4^2)*f3);

Hwz = c2d(Hw, Ts, 'tustin');

Hwzn = Hwz.Numerator;
Hwznv = cell2mat(Hwzn);

Hwzd = Hwz.Denominator;
Hwzdv = cell2mat(Hwzd);

% I made them into vectors directly with cell2mat
B1 = Hbznv;
A1 = Hbzdv;

B3 = Hwznv;
A3 = Hwzdv;

filtered_signal=filtfilt(B1,A1,acceleration_signal);
filtered_signal=filtfilt(B3,A3,filtered_signal);

I hope this is what you meant. I also tried putting these vectors in the bilinear function but got the same errors, i figured they were not needed if the tustin is the same as bilinear. But either way i get the errors regarding the badly scaled matrices. Its weird since it should basically be the same as before right?

[Star Strider]

The sampling frequency may be too high to produce reliable filters. Is that high a sampling frequency absolutely required, or is it possible to use perhaps a tenth of that?

Although the Control System Toolbox bodeplot plots were encouraging, I am not happy with the way the filters are transferring to the discrete domain, so i tried with the Symbolic Math Toolbox as well to see if I could improve on the result.

I am unfortunately not encouraged —

Fs=65536; %Sample rate
Fn = Fs/2;
f1=6.310/Fn;                                % Normalise The Frequency By The Nyquist Frequency
w1=2*pi*f1;
f2=1258.9/Fn;                               % Normalise The Frequency By The Nyquist Frequency
w2=2*pi*f2;
f3=15.915/Fn;                               % Normalise The Frequencies By The Nyquist Frequency
f4=f3;
Q1=1/sqrt(2);
Q2=Q1;
Q4=0.64;
w3=2*pi*f3;
w4=2*pi*f4;
K = 1;
Ts = 1/Fs;   

acceleration_signal = randn(Fs*5,1);

% s = tf('s');
syms s
Hb = (s^2*pi^2*f2^2) / ((s^2 + 2*pi*f1*s/Q1 + 4*pi^2*f1^2)*(s^2 + 2*pi*f2*s/Q1 + 4*pi^2*f2^2));
Hb = simplify(Hb, 500)

Hb = 

[Hbn,Hbd] = numden(Hb)

Hbn = 

Hbd = 

Hbn = double(sym2poly(Hbn))

Hbn = 1×3

1.0e+36 * 9.9141 0 0

Hbd = double(sym2poly(Hbd))

Hbd = 1×5

1.0e+38 * 6.8056 2.3349 0.4005 0.0007 0.0000

[B1,A1] = bilinear(Hbn,Hbd,Fs)

B1 = 1×5

1.0e-11 * 0.0848 0.0001 -0.1700 0.0004 0.0847

A1 = 1×5

1.0000 -4.0000 6.0000 -4.0000 1.0000

figure
freqz(B1,A1, 2^16, Fs)
set(subplot(2,1,1), 'XLim',[0 2])
set(subplot(2,1,2), 'XLim',[0 2])

Hw = ((s + 2*pi*f3)*2*pi*K*f4^2) / ((s^2 + 2*pi*f4*s/Q2 + 4*pi^2*f4^2)*f3);
Hw = simplify(Hw, 500)

Hw = 

[Hwn,Hwd] = numden(Hw)

Hwn = 

Hwd = 

Hwn = double(sym2poly(Hwn))

Hwn = 1×2

1.0e+33 * 6.4557 0.0197

Hwd = double(sym2poly(Hwd))

Hwd = 1×3

1.0e+36 * 2.1155 0.0091 0.0000

[B3,A3] = bilinear(Hwn,Hwd,Fs)

B3 = 1×3

1.0e-07 * 0.2328 0.0000 -0.2328

A3 = 1×3

1.0000 -2.0000 1.0000

figure
freqz(B3,A3, 2^16, Fs)
set(subplot(2,1,1), 'XLim',[0 2])
set(subplot(2,1,2), 'XLim',[0 2])

return                                  % Stop Here

Hbz = c2d(Hb, Ts, 'tustin');
Hbzn = Hbz.Numerator;
Hbznv = cell2mat(Hbzn);
Hbzd = Hbz.Denominator;
Hbzdv = cell2mat(Hbzd);

Hw = ((s + 2*pi*f3)*2*pi*K*f4^2) / ((s^2 + 2*pi*f4*s/Q2 + 4*pi^2*f4^2)*f3);

Hwz = c2d(Hw, Ts, 'tustin');

Hwzn = Hwz.Numerator;
Hwznv = cell2mat(Hwzn);

Hwzd = Hwz.Denominator;
Hwzdv = cell2mat(Hwzd);

% I made them into vectors directly with cell2mat
B1 = Hbznv;
A1 = Hbzdv;

B3 = Hwznv;
A3 = Hwzdv;

[sos1,g1] = tf2sos(B1,A1)
[sos3,g3] = tf2sos(B3,A3)

figure
freqz(B1,A1, 2^16, Fs)
set(subplot(2,1,1), 'XLim',[0 2])
set(subplot(2,1,2), 'XLim',[0 2])

figure
freqz(sos1, 2^16, Fs)
set(subplot(2,1,1), 'XLim',[0 2])
set(subplot(2,1,2), 'XLim',[0 2])

figure
freqz(B3,A3, 2^16, Fs)

figure
freqz(sos3, 2^16, Fs)

filtered_signal=filtfilt(sos1,g1,acceleration_signal);
filtered_signal=filtfilt(sos3,g3,filtered_signal);

I need to be away for a while. I will return to this later.

.