BS 6841 frequency weighting functions

Dear expertises,

I am now trying to code the weighting function in BS 6841:1987 - Whole body vibration Evaluation. However, I cannot code correctly and don't know what is wrong. I am, therefore, really hope that you can figure out my problem. Thank you very much for your help.

Here is my code.

clear all
clc
f1=0.4; f2=100;q1=0.71; q2=0.55;  f3=16; q3=0.9; f4=16; q4=0.95; f5=2.5; f6=4; k=0.4; %Wb
w1=2*pi*f1; w2=2*pi*f2; w3=2*pi*f3; w4=2*pi*f4; w5=2*pi*f5; w6=2*pi*f6;

N=61440; % Number of data points
Fs=100; % Sampling rate
f=linspace(0,Fs/2,N/2);
Hb=(f.^2./((f.^2+w1/q1*f+w1^2))).*(w2^2./(f.^2+w2/q1*f+w2^2));
Hw=((f+w3).*(f.^2+w5/q3*f+w5^2))./((f.^2+w4/q2*f+w4^2).*(f.^2+w6/q4*f+w6^2))*2*pi*k*f4^2*f6^2/(f3*f5^2);

Wb=Hb.*Hw;
loglog(f,Wb)