function [yf, ww]  = bl_filt(y, half_width)
% function [yf, ww] = bl_filt(y, half_width)
%
% Blackman filter y.  The filter slides over the whole array,
% using only the non-NaN values but yielding results for all
% elements, so long as the number of consecutive points does not exceed
% twice the half-width.
%
% ww is the filter fraction with data; it is 1 if there are no gaps within
% the filter width, 0.5 if half the filter integral covered gaps,
% NaN (like yf) if there were no good points within the filter range.
%
% Note: the interpolation becomes shaky for small values of the filter
% fraction; either NaNs or some other method of interpolation should
% be used for filter fractions less than something like 0.3-0.5.
%
% For now, this is set up to work only with vectors.
%
% E.F. 94/12/10  (R/V Knorr, WHP I8S)
% JH - force output size to match input size
%

if (nargin ~= 2)
  help bl_filt
  error('wrong number of parameters');
end

[nrows, ncols] = size(y);


%--- generate weights array for given width ---
nf = half_width * 2 + 1;
ny = length(y);

x = (-half_width:half_width)*pi/half_width;
w = 0.42 + 0.5 * cos(x) + 0.08 * cos(2*x);

bad_mask = isnan(y);
good_mask = ~bad_mask;

% Replace NaNs with zeros.
ib = find(bad_mask);
if ~isempty(ib),
   y(ib) = good_mask(ib);
end

ytop = conv(y, w);
ybot = conv(good_mask, w);

% If a gap
ibad = find(ybot == 0);
if (~isempty(ibad))
   ybot(ibad) = NaN * ybot(ibad);
end

% Select the middle ny points of the convolution.
ii = half_width + (1:ny);

yf = ytop(ii) ./ ybot(ii);

if nargout > 1,
   ww = ybot(ii) / sum(w);
end


%swap directions if output was swapped.  Note that the directions disagree
%   only if 2*half_width < length(y)

[mrows, mcols] = size(yf);
if (mrows == ncols)
   yf = yf.';
   if (nargout > 1)
      ww = ww.';
   end
end

