[filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type, k)

function [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type, k)
% % whole_Body_Filter: Calculates the filter coefficients for whole body vibrations.  
% % 
% % Syntax;
% % 
% % [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type, k);
% % 
% % ********************************************************************
% %
% % Description
% % 
% % [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs,
% % type, k);
% % 
% % Returns the filter coefficients, filter types, and scaling factors for 
% % whole body vibrations according to ISO 2631-1:1997(E) given the
% % sampling rate Fs (Hz), the posture type, and the optional variable k. 
% % 
% % The input variables and output variablses are described in more detail
% % in the respective sections below.  
% % 
% % ********************************************************************
% %
% % Input Variables
% % 
% % Fs (Hz) is the sampling frequency.  default is Fs=5000;
% % 
% % type is an integer for the whole body vibration posture.
% %       For the selected posture types the corresponding transfer 
% %       functions, and scaling factors are collected.  The transfer 
% %       functions can be used to compute the filter coefficients 
% %       (see whole_body_time_filter for more information). 
% %       default is type=1;
% % 
% %      The table below relates the type integer to each posture.
% % 
% %      % type=1;    Standing
% %      % type=2;    Seated (Health)
% %      % type=3;    Seated (Comfort)
% %      % type=4;    Laying on back (recumbent) k=1 Pelvis, k=2 Head
% %      % type=5;    Rotational on supporting seat
% %      % type=6;    Rotational on seated backrest
% %      % type=7;    Rotational on feet
% %      % type=8;    Motion sickness
% % 
% % 
% %    There are six filter types which are used to calculate the transfer 
% %    functions for the seven postures.  The seven posture types and six 
% %    filter types are easily confused.   The posture types may use more 
% %    than one filter type.  
% %    
% %    The case structure in this program documents which postures
% %    correspond to the filters.  For example: Posture number 7 
% %    (Motion sickness) uses filter number 4, Wf.   
% %    
% %    Table 1:  ISO 2631 has 6 filter types
% %    Wc Wd We Wf Wk Wl
% %    1  2  3  4  5  6 
% % 
% % 
% % k    ISO 2631 uses the variable k for multiple purposes.  
% %      To avoid overloading the definition of k, in these programs, 
% %      k is only used for posture type 4, more specifically the 
% %      recumbent posture z-axis recording position.
% % 
% %      k specifies the pelvis or head if the recumbent posture was used.
% %      k is 1 for Pelvis vibration recordings 
% %      k is 2 for Head vibration recordings
% %      default is k=1;
% % 
% % 
% % ********************************************************************
% % 
% % Output Variables
% % 
% % filtercoeffs is a cell array of cell arrays of cell arrays of arrays
% %              of filter coefficients.   
% % 
% %              filtercoeffs{filter_type}{transfer_functions}{1}=B;
% %              filtercoeffs{filter_type}{transfer_functions}{2}=A;
% % 
% %              There are three transfer_functions: the Band limiting 
% %              filter (high pass and low pass), Acceleration-Velocity 
% %              transition, and the Upward step. 
% %              
% % 
% % filter_types There are six filter types
% % 
% % Wc Wd We Wf Wk Wl
% % 1  2  3  4  5  6 
% % 
% % scaling_factors are the factors namely: 0, 1, sqrt(2), or sqrt(3) which 
% %             are used to scale a single axis measurement to a multiple 
% %             axis estimation of overall acceleration.  
% % 
% % ********************************************************************
% 
% 
% 
% Example='1';
% 
% % Standing Posture 
% % Using a sampling rate of 5000 Hz 
% 
% Fs=5000;  type=1;  
% [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type);
% 
% 
% 
% Example='2';
% % Seated Posture (Health)
% % Using a sampling rate of 5000 Hz 
% 
% Fs=5000;  type=2;  
% [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type);
% 
% 
% 
% Example='3';
% % Seated Posture (Comfort)
% % Using a sampling rate of 5000 Hz 
% 
% Fs=5000;  type=3;  
% [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type);
% 
% 
% 
% Example='4';
% % Laying on Back Posture recording at Pelvis location
% % Using a sampling rate of 5000 Hz 
% 
% Fs=5000;  type=4; k=1;   
% [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type);
% 
% 
% 
% Example='5';
% % Laying on Back Posture recording at Head location
% % Using a sampling rate of 5000 Hz 
% 
% Fs=5000;  type=4; k=2;   
% [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type);
% 
% 
% 
% Example='6';
% % Rotational on supporting seat Posture 
% % Using a sampling rate of 5000 Hz 
% 
% Fs=5000;  type=5;  
% [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type);
% 
% 
% 
% Example='7';
% % Rotational on seated backrest
% % Using a sampling rate of 5000 Hz 
% 
% Fs=5000;  type=6;  
% [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type);
% 
% 
% 
% Example='8';
% % Rotational on feet
% % Using a sampling rate of 5000 Hz 
% 
% Fs=5000;  type=7;  
% [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type);
% 
% 
% 
% Example='9';
% % Motion Sickness Posture 
% % Using a sampling rate of 100 Hz 
% 
% Fs=100;  type=8; 
% [filtercoeffs, filter_types, scaling_factors]=whole_Body_Filter(Fs, type);
% 
% 
% % ********************************************************************
% %
% % References For whole Body Vibrations
% %
% % ISO 2631-1 1997  Mechanical vibration and shock-Evaluation of human
% %                  exposure to whole-body vibration-Part 1:
% %                  General Requirements
% %
% % ********************************************************************
% %
% % Program Written by Edward L. Zechmann
% %
% %     date    July        2007
% % 
% % modified 3  September   2008
% % 
% % modified 16 December    2008    Updated comments
% %                               
% % modified 21 January     2009    Split the seated posture into two cases 
% %                                 Seated (Health) and Seated (Comfort). 
% %                                 Only documentation needed adjustment.
% %                               
% % modified 17 January     2011    Fixed a bug in the splitting the seated
% %                                 posture into two cases.  The scaling 
% %                                 factors for the two cases were switched.
% % 
% % ********************************************************************
% % 
% % Please Feel Free to Modify This Program
% % 
% % See Also:   whole_body_time_filter, combine_accel_directions_wb, 
% %             config_wb_accels, Vibs_calc_whole_body
% % 

if (nargin < 1 || isempty(Fs)) || ~isnumeric(Fs) 
    Fs=5000;
end

if (nargin < 2 || isempty(type)) || ~isnumeric(type) 
    type=1;
end

if (nargin < 3 || isempty(k)) || ~isnumeric(k) 
    k=1;
end



filtercoeffs={};
filter_types=[];
scaling_factors=[];

% From ISO 2631

switch type
    
    case 1
        % for standing people
        kx=1;   %Wd;     x-axis
        ky=1;   %Wd;     y-axis
        kz=1;   %Wk;     z-axis
        
        filter_types=[2 5];
        scaling_factors=[kx ky 0; 0 0 kz];
        
    case 2
        % Scaling factors for Health
        kx=1.4;   %Wd;     %x-axis
        ky=1.4;   %Wd;     %y-axis
        kz=1;     %Wk;     %z-axis

        filter_types=[2 5];
        scaling_factors=[kx ky 0; 0 0 kz];
        
    case 3
        % for seated persons
        % Scaling factors for Comfort
        
        kx=1;   %Wd;     %x-axis
        ky=1;   %Wd;     %y-axis
        kz=1;   %Wk;     %z-axis

        filter_types=[2 5];
        scaling_factors=[kx ky 0; 0 0 kz];
        
    case 4
        % Laying on back (recumbent)
        kx=1;   %Wd;     %x-axis (under pelvis)
        ky=1;   %Wd;     %y-axis (under pelvis)     
        kz=1;   %Wk;     %z-axis (under pelvis)
        kz=1;   %Wl;     %z-axis (under head) k=1;
        
        if nargin < 3
            k = menu('Is the Vibration Measurement under the Pelvis or Head?','Pelvis','Head');
        end
        
        if isequal(k, 1)
            filter_types=[2 5];
        else
            filter_types=[2 6];
        end

        scaling_factors=[kx ky 0; 0 0 kz];
        
    case 5
        %rotational on supporting seat
        kx=0.63;    %We;     %x-axis
        ky=0.4;     %We;     %y-axis
        kz=0.2;     %We;     %z-axis
        
        filter_types=[3];
        scaling_factors=[kx ky kz];
        
    case 6
        %rotational on seated backrest
        kx=0.8;     %Wc;     %x-axis
        ky=0.5;     %Wd;     %y-axis
        kz=0.4;     %Wd;     %z-axis
        
        filter_types=[1 2];
        scaling_factors=[kx 0 0; 0 ky kz];
        
    case 7 
        %rotational on feet
        kx=0.25;    %Wk;     %x-axis
        ky=0.25;    %Wk;     %y-axis
        kz=0.4;     %Wk;     %z-axis
        
        filter_types=[5];
        scaling_factors=[kx ky kz];
        
    case 8 
        % Motion sickness
        kz=1;%        Wf;

        filter_types=[4];
        scaling_factors=[kz kz kz];
        
    otherwise
        kx=1;   %Wd;     %x-axis
        ky=1;   %Wd;     %y-axis
        kz=1;   %Wk;     %z-axis
        
        filter_types=[2 5];
        scaling_factors=[kx ky 0; 0 0 kz];
        
end
       
% Constants for Wc Wd We Wf Wk Wl
% This is an array containing all of the filter specification constants
% According to ISO 2631-1-1997
% The large numbers 1000000000 approximate infinity.  Later in the program
% booleans are used to detect the large numbers and use the exact formula 
% to calculate the filter coefficients.  
% 
f1=[    0.4         0.4         0.4         0.08        0.4     0.4         ];
f2=[    100         100         100         0.63        100     100         ];
f3=[    8           2           1           1000000000	12.5	1000000000  ];
f4=[    8           2           1           0.25        12.5	1000000000  ];
Q4=[    0.63        0.63        0.63        0.86        0.63	0           ];
f5=[    1000000000	1000000000	1000000000	0.0625      2.37	3.75        ];
Q5=[    0           0           0           0.8         0.91	0.91        ];
f6=[    1000000000	1000000000	1000000000	0.1         3.35	5.32        ];
Q6=[    0           0           0           0.8         0.91	0.91        ];


for e1=1:length(filter_types);

    % band pass filter implemented by a low pass and a high pass filter
    e2=1; % High Pass Filter
    
    ft=filter_types(e1);
    
    w1=2*pi*f1(ft);
    
    num1=[1 0 0];
    den1=[1 2^(0.5)*w1 w1^2];

    [B, A] = bilinear(num1, den1, Fs);
    filtercoeffs{e1}{e2}{1}=B;
    filtercoeffs{e1}{e2}{2}=A;
    
    e2=2; % Low Pass Filter
    
    w2=2*pi*f2(ft);
    
    num1=[1];
    den1=[ (1/w2)^2 2^(0.5)/w2 1];
    
    [B, A] = bilinear(num1, den1, Fs);
    
    filtercoeffs{e1}{e2}{1}=B;
    filtercoeffs{e1}{e2}{2}=A;
    
    e2=3; % acceleration-velocity transition filter
    
    w3=2*pi*f3(ft);
    w4=2*pi*f4(ft);
    Q44=Q4(ft);
    
    if f3(ft) > 20
        num1=[1];
    else
        num1=[1/w3 1];
    end
    
    if f4(ft) > 20
        den1=[1];
    else
        den1=[ (1/w4)^2 1/(Q44*w4) 1];
    end
    
    if length(den1) < 3
        B=1;
        A=1;
    else
        [B, A] = bilinear(num1, den1, Fs);
    end
    
    filtercoeffs{e1}{e2}{1}=B;
    filtercoeffs{e1}{e2}{2}=A;
    
    e2=4; % upward step filter
    
    w5=2*pi*f5(ft);
    w6=2*pi*f6(ft);
    Q55=Q5(ft);
    Q66=Q6(ft);
    
    if f5(ft) > 20
        num1=[1];
        den1=[1];
    else
        num1=(w5/w6)^2*[ (1/w5)^2 1/(Q55*w5) 1];
        den1=[ (1/w6)^2 1/(Q66*w6) 1];
    end
    
    if length(den1) < 3
        B=1;
        A=1;
    else
        [B, A] = bilinear(num1, den1, Fs);
    end
    
    filtercoeffs{e1}{e2}{1}=B;
    filtercoeffs{e1}{e2}{2}=A;
    
end