Issues with creating a frequency weighting filter for a iso-standard
25 views (last 30 days)
Show older comments
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 :)
12 Comments
Star Strider
on 17 Apr 2023
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.
Answers (1)
Star Strider
on 17 Apr 2023
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.
.
4 Comments
Star Strider
on 18 Apr 2023
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)
[Hbn,Hbd] = numden(Hb)
Hbn = double(sym2poly(Hbn))
Hbd = double(sym2poly(Hbd))
[B1,A1] = bilinear(Hbn,Hbd,Fs)
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)
[Hwn,Hwd] = numden(Hw)
Hwn = double(sym2poly(Hwn))
Hwd = double(sym2poly(Hwd))
[B3,A3] = bilinear(Hwn,Hwd,Fs)
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.
.
See Also
Categories
Find more on Bartlett in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!