nonlinFiltAll(inSig, inhw, trct, wfc1, wfc2, wfc3)

function fs = nonlinFiltAll(inSig, inhw, trct, wfc1, wfc2, wfc3)

%----------------------------------------------------
% Compute filtered signal using nonlinear target filter (combination of
% three filers : see [Malyarenko, et al RCMS 2006])
%
% INPUT: 
% inSig - input TOF signal (resampled to have constant point density per
% peak)
% inhw - HWHM of input signal
% trct - ADC precision level for artifact truncation
% wfc1,2,3 - filter coefficients for 3 filters (columns) to be geometrically
% averaged
%
% OUTPUT:
% fs - filtered signal
%
% USAGE:
% fs = nonlinFiltAll(inSig, inhw, trct, wfc1, wfc2, wfc3);
% Dependency: "filter" function from "signal" toolbox

sins = size(inSig);
if sins(2) > sins(1) inSig=inSig'; end;

swfc1 = size(wfc1);
if swfc1(2) > swfc1(1) wfc1=wfc1'; end;

swfc2 = size(wfc2);
if swfc2(2) > swfc2(1) wfc2=wfc2'; end;

swfc3 = size(wfc3);
if swfc3(2) > swfc3(1) wfc3=wfc3'; end;

Nw = max([length(wfc1),length(wfc2),length(wfc3)]);

fs1 = filter(wfc1, 1, inSig); 

Nw1=length(wfc1);
if (Nw1/2==fix(Nw1/2)) sh1 = -Nw1/2; % for filter symmetry (around maximum)
else sh1 = -(Nw1-1)/2;
end;
fs1 = circShift(fs1, sh1);

fs2 = filter(wfc2, 1, inSig); 

Nw2=length(wfc2);
if (Nw2/2==fix(Nw2/2)) sh1 = -Nw2/2; % for filter symmetry (around maximum)
else sh1 = -(Nw2-1)/2;
end;
fs2 = circShift(fs2, sh1);

fs3 = filter(wfc3, 1, inSig); 

Nw3=length(wfc3);
if (Nw3/2==fix(Nw3/2)) sh1 = -Nw3/2; % for filter symmetry (around maximum)
else sh1 = -(Nw3-1)/2;
end;
fs3 = circShift(fs3, sh1);

fs = sign(fs1.*fs2.*fs3).*((abs(fs1.*fs2.*fs3)).^(1/3.)); 
%thrs=-max(fs)/trct; 
i1 = find((fs1.*fs2<0) | (fs1.*fs3<0) | (fs2.*fs3<0) | (fs.*fs1<0) | (fs.*fs2<0) | (fs.*fs3<0) | (fs <0) | (fs1 < 0) | (fs2 < 0) | (fs3 < 0));
fs(i1) = fs(i1)/(max(fs)*trct); % artifact truncation
%% may use noise level fro artifacts instead

l1 = length(fs);
%% correct for filter shift into early points 11/14/07
fs((end-floor(Nw/2)):end) = inSig((end-floor(Nw/2)):end); 
fs = fs(1:l1);
return;