[yh, B, A, errors]=hand_arm_time_fil(y, Fs, B, A, settling_time, resample_filter)

function [yh, B, A, errors]=hand_arm_time_fil(y, Fs, B, A, settling_time, resample_filter)
% % hand_arm_time_fil: Calculates the hand arm vibrations filter coefficients and returns the filtered time record
% %
% % Syntax;
% %
% % [yh, B, A, errors]=hand_arm_time_fil(y, Fs, B, A, settling_time, resample_filter)
% %
% % ********************************************************************
% %
% % Description
% %
% % Calculates the hand arm vibrations filter coefficients and returns the
% % filtered time record.
% %
% %
% % There are two options for the downsampling filters to optimize
% % performance for continuous signals or for impulsive signals. 
% % For continuous noise the time domain does not have significant 
% % impulses; however, for impulsive time records there are often very
% % large impulses with distinctive peaks.  
% % 
% % There are two antialiasing filters and interpolation schemes available.
% % The first program is the built-in Matlab "resample" progam which
% % uses a Kaiser window fir filter for antialising and uses an unknown 
% % interpolation method.  The second program available for downsampling 
% % is bessel_down_sample which uses a Bessel filter for antialiasing 
% % and uses interp with the cubic spline option for interpolation.  
% % 
% % The resample function has good antialising up to the Nyquist frequency;
% % however, it has significant ringing effect when there are impulses.  
% % The bessel_down_sample function has good antialising; however, there is
% % excessive attenuation near the Nyquist frequency.  
% % The bessel_down_sample function experiences no ringing due to impulses
% % so it is very useful for peak estimation.  
% % 
% % 
% % The input and output variables are described in more detail in the
% % sections below respectively.
% %
% %
% % ********************************************************************
% %
% % Input Variables
% %
% % y=randn(10000, 10)  % multichannel input time record in (Pa).
% %                     % Processsing assumes that y has more channels
% %                     % than time record samples.
% %                     % default is y=randn(10000, 10);
% %
% % Fs=50000;           % (Hz) sampling rate in Hz.
% %                     % default is 50000 Hz.
% %
% % B is an array of feedforward filter coefficients for the hand arm vibration
% % 
% % A is an array of feedback filter coefficients for the hand arm vibration
% %                     % frequency weighting filter.
% %                     % default is B=1; A=1; The program recalculates
% %                     % B and A if the input values are not numeric,
% %                     % are empty, or have a length less than 4.
% %
% % settling_time=0.1;  % (seconds) is the time it takes the filter to
% %                     % settle (seconds).
% %                     % default is settling_time=0.1;
% %
% % resample_filter=1;  % type of filter to use when resamling
% %                     % 1 resample function Kaiser window fir filter
% %                     % 2 Bessel filter
% %                     % otherwise resample function Kaiser window fir
% %                     % filter
% %                     % default is resample_filter=1; (Kaiser window)
% %
% % ********************************************************************
% %
% % Output Variables
% %
% % yh is the time record which has been filtered by the hand arm
% %             vibration frequency filter.
% % 
% % B is an array of feedforward filter coefficients for hand arm
% %             vibration frequency weighting filter.
% % 
% % A is an array of feedback filter coefficients for hand arm
% %             vibration frequency weighting filter.
% %
% % errors indicates whether the resampled sampling rate is within the
% %             tolerance of 5000 Hz.  Within tolerance 0 is output.
% %             outside of tolerance 1 is output.
% %
% % **********************************************************************
%
%
% Example='1';
%
% % An exaple with 1 channels of vibrations data
%
% y=randn(1, 50000);
% Fs=10000;
% B=1;
% A=1;
% settling_time=1;
%
% [yh, B, A, errors]=hand_arm_time_fil(y, Fs, B, A, settling_time);
%
% % Plot the results
% t_vibs=1/Fs*(0:(50000-1));
% subplot(2,1,1);
% plot(t_vibs, y);
% title('Unweighted Time Record');
% hold on;
% subplot(2,1,2);
% plot(t_vibs, yh);
% title('Weighted Time Record');
%
% % ********************************************************************
% %
% % This program requires the Matlab Signal Processing Toolbox
% %
% % 
% % List of Dependent Subprograms for 
% % hand_arm_time_fil
% % 
% % FEX ID# is the File ID on the Matlab Central File Exchange
% % 
% % 
% % Program Name   Author   FEX ID#
% %  1) bessel_antialias		Edward L. Zechmann			
% %  2) bessel_digital		Edward L. Zechmann			
% %  3) bessel_down_sample		Edward L. Zechmann			
% %  4) convert_double		Edward L. Zechmann			
% %  5) filter_settling_data3		Edward L. Zechmann			
% %  6) geospace		Edward L. Zechmann			
% %  7) get_p_q2		Edward L. Zechmann			
% %  8) hand_arm_fil		Edward L. Zechmann			
% %  9) LMSloc		Alexandros Leontitsis		801	
% % 10) match_height_and_slopes2		Edward L. Zechmann			
% % 11) remove_filter_settling_data		Edward L. Zechmann			
% % 12) resample_interp3		Edward L. Zechmann			
% % 13) rms_val		Edward L. Zechmann			
% % 14) sub_mean2		Edward L. Zechmann			
% %
% %
% % ********************************************************************
% %
% % Program Modified by Edward L. Zechmann
% %
% %  created    June        2007
% %
% % modified  3 September   2008    Updated comments
% %
% % modified 10 December    2008    Upgraded to include
% %                                 filter settling and resampling.
% %
% % modified 16 December    2008    Use convolution to simplify filter
% %                                 coefficients (B1, B2, A1, A2) into
% %                                 arrays B and A.
% %
% %                                 Finished modifications to support
% %                                 filter settling and resampling.
% %
% % modified 19 December    2008    Updated comments
% %
% % modified  5 January     2009    Added sub_mean to remove running
% %                                 average using a a number of averages
% %                                 per second at one-half the lowest
% %                                 frequency band of interest(1 Hz).
% %
% % modified  9 October     2009    Updated Comments
% % 
% % modified  9 July        2010    Added an option to resample using a
% %                                 Bessel Filter.  Updated comments.
% %
% % modified  4 August      2010    Updated Comments
% % 
% %
% %
% %
% % ********************************************************************
% %
% % Please Feel Free to Modify This Program
% %
% % See Also: hand_arm_fil, vibs_calc_hand_arm
% %


if nargin < 1 || isempty(y) || ~isnumeric(y)
    y=randn(10, 10000);
end

if nargin < 2 || isempty(Fs) || ~isnumeric(Fs)
    Fs=10000;
end

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

if nargin < 4 || isempty(A) || ~isnumeric(A)
    A=1;
end

% Hand-arm vibrations require a longer filter settling time than for sound
% because the important frequencies are at least an order of magnitude
% lower.
if (nargin < 5 || isempty(settling_time)) || ~isnumeric(settling_time)
    settling_time=1;
end

if (nargin < 6 || isempty(resample_filter)) || ~isnumeric(resample_filter)
    resample_filter=1;
end

if ~isequal(resample_filter, 2)
    resample_filter=1;
end


% Remove the running average from the signal.
% The number of averages per second can be about one-tenth the lowest
% frequency to resolve.  For hand-arm vibrations, the lowest frequency to
% resolve according to ISO 5349-1 is 1.0 Hz
[y]=sub_mean2(y, Fs, 0);



% % ********************************************************************
% %
% set the flag for indicating whether to resample
flag0=0;

max_freq=60000;
min_freq=5000;

if Fs <= max_freq
    p=1;
    q=1;
    Fsn=Fs;
    errors=0;
else
    [Fsn, p, q, errors]=get_p_q2(Fs, max_freq, min_freq);
end

if ~isequal(Fsn, Fs)
    flag0=1;
end



% % ********************************************************************
% %

[num_samples, num_channels]=size(y);

% set the flag for indicating whether to transpose
flag1=0;

% Transpose the data so that the filters act along the time records
% not along the channels.
if num_samples < num_channels
    flag1=1;
    y=y';
    [num_samples, num_channels]=size(y);
end


if (nargin < 4) || (logical(length(B) < 4) ||  logical(length(A) < 4)) || logical(~isnumeric(B) || ~isnumeric(A))

    [B, A]=hand_arm_fil(Fsn);

end


% % ********************************************************************
%
% Downsample the data to a reasonable sampling rate to keep the hand arm
% vibrations filters as stable as possible.
if isequal(flag0, 1)

    % Adding filter settling data
    num_data_pts=num_samples;
    [y11, num_settle_pts]=filter_settling_data3(Fs, y(:, 1), settling_time);

    % Downsampling
    if isequal(resample_filter, 1)
        buf=resample(y11, p, q);
    else
        [buf]=bessel_down_sample(y11, Fs, Fs*p/q, settling_time);
    end

    % removing filter settling data
    [buf]=remove_filter_settling_data(buf, q/p, num_settle_pts, num_data_pts, 0);

    n2=length(buf);
    y2=zeros(n2, num_channels);
    y2(:, 1)=buf;
    clear('buf');
    if num_channels > 1
        for e1=2:num_channels;

            % Adding filter settling data
            [y11, num_settle_pts]=filter_settling_data3(Fs, y(:, e1), settling_time);

            % Downsampling
            if isequal(resample_filter, 1)
                buf=resample(y11, p, q);
            else
                [buf]=bessel_down_sample(y11, Fs, Fs*p/q, settling_time);
            end
            
            % removing filter settling data
            [buf]=remove_filter_settling_data(buf, q/p, num_settle_pts, num_data_pts, 0);
            y2(:, e1)=buf;
        end
    end
else
    y2=y;
end

clear('y');


% % ********************************************************************
%
% Apply the hand-arm weighting filters to each channel of data

yh=zeros(size(y2));

for e1=1:num_channels;
    
    % Adding filter settling data
    num_data_pts=length(y2(:, e1));
    [y11, num_settle_pts]=filter_settling_data3(Fsn, y2(:, e1), settling_time);

    % Apply the hand-arm vibrations filters
    buf = real(filter(B, A, y11 ));

    % Remove the settling data from the time record.
    [buf]=remove_filter_settling_data(buf, 1, num_settle_pts, num_data_pts, 2);    
    yh(:, e1)=buf;
    
end


% % ********************************************************************
%
% Resample if necessary so the output has the same size as the input
if isequal(flag0, 1)

    % Adding filter settling data
    num_data_pts=length(yh(:, 1));
    [y11, num_settle_pts]=filter_settling_data3(Fs, yh(:, 1), settling_time);

    % Upsampling
    if isequal(resample_filter, 1)
        buf=resample(y11, q, p);
    else
        t_in=q/p*1/Fs*(-1+(1:length(y11)));
        [buf]=resample_interp3(y11, t_in, 1/Fs);
        [buf]=bessel_antialias(buf, Fs, Fs/2.5, settling_time);
    end


    % removing filter settling data
    [buf]=remove_filter_settling_data(buf, q/p, num_settle_pts, num_data_pts, 1);


    num_samples2=length(buf);
    yhh=zeros(num_samples2, num_channels);
    yhh(:, 1)=buf;
    clear('buf');
    if num_channels > 1
        for e1=2:num_channels;

            % Adding filter settling data
            num_data_pts=length(yh(:, 1));
            [y11, num_settle_pts]=filter_settling_data3(Fs, yh(:, e1), settling_time);

            % Upsampling
            if isequal(resample_filter, 1)
                buf=resample(y11, q, p);
            else
                t_in=q/p*1/Fs*(-1+(1:length(y11)));
                [buf]=resample_interp3(y11, t_in, 1/Fs);
                [buf]=bessel_antialias(buf, Fs, Fs/2.5, settling_time);
            end

            % removing filter settling data
            [buf]=remove_filter_settling_data(buf, q/p, num_settle_pts, num_data_pts, 1);

            yhh(:, e1)=buf;
        end
    end
else
    yhh=yh;
end

clear('yh');


% % ********************************************************************
%
% Append the data to the output variable
if isequal(flag0, 1)
    if num_samples2 > num_samples
        num_samples2=num_samples;
    end

    yh=yhh(1:num_samples2, 1:num_channels);
else
    yh=yhh(1:num_samples, 1:num_channels);
end


% % ********************************************************************
%
% Output yh must have the same size as the original input matrix y.
% Transpose yh if necessary.
if isequal(flag1, 1)
    yh=yh';
end