function gout=blfilter(winname,fsupp,varargin)
%BLFILTER  Construct a band-limited filter
%   Usage:  g=blfilter(winname,fsupp,fc);
%           g=blfilter(winname,fsupp,fc,...);
%
%   Input parameters:
%      winname  : Name of prototype
%      fsupp    : Support length of the prototype
%
%   `blfilter(winname,fsupp)` constructs a band-limited filter. The parameter
%   *winname* specifies the shape of the frequency response. The name must be
%   one of the shapes accepted by |firwin|. The support of the frequency
%   response measured in normalized frequencies is specified by *fsupp*.
%
%   `blfilter(winname,fsupp,fc)` constructs a filter with a centre
%   frequency of *fc* measured in normalized frequencies.
%
%   If one of the inputs is a vector, the output will be a cell array
%   with one entry in the cell array for each element in the vector. If
%   more input are vectors, they must have the same size and shape and the
%   the filters will be generated by stepping through the vectors. This
%   is a quick way to create filters for |filterbank| and |ufilterbank|.
%
%   `blfilter` accepts the following optional parameters:
%
%     'fs',fs     If the sampling frequency *fs* is specified then the support
%                 *fsupp* and the centre frequency *fc* is specified in Hz.
%
%     'complex'   Make the filter complex valued if the centre frequency
%                 is non-zero. This is the default.
%
%     'real'      Make the filter real-valued if the centre frequency
%                 is non-zero.
%
%     'delay',d   Set the delay of the filter. Default value is zero.
%
%     'scal',s    Scale the filter by the constant *s*. This can be
%                 useful to equalize channels in a filter bank.
%
%     'pedantic'  Force window frequency offset (g.foff) to a subsample 
%                 precision by a subsample shift of the firwin output.
%                 
%
%   It is possible to normalize the transfer function of the filter by
%   passing any of the flags from the |normalize| function. The default
%   normalization is `'energy'`.
%
%   The filter can be used in the |pfilt| routine to filter a signal, or
%   in can be placed in a cell-array for use with |filterbank| or
%   |ufilterbank|.
%
%   Output format:
%   --------------
%
%   The output *g* from `blfilter` is a structure. This type of structure can
%   be used to describe any bandlimited filter defined in terms of its
%   transfer function. The structure contains the following fields:
%
%     g.H     This is an anonymous function taking the transform length *L* as
%             input and producing the bandlimited transfer function in the
%             form of a vector.
%
%     g.foff  This is an anonymous function taking the transform length *L* as
%             input and procing the frequency offset of *H* as an integer. The
%             offset is the value of the lowest frequency of *H* measured in
%             frequency samples. *foff* is used to position the bandlimited
%             tranfer function stored in *H* correctly when multiplying in the
%             frequency domain.
%
%     g.delay  This is the desired delay of the filter measured in samples.
%
%     g.realonly
%             This is an integer with value *1* if the filter defined a
%             real-valued filter. In this case, the bandlimited transfer
%             function *H* will be mirrored from the positive frequencies to
%             the negative frequencies. If the filter is a natural lowpass
%             filter correctly centered around *0*, `realonly` does not need
%             to be *1*.
%
%     g.fs    The intended sampling frequency. This is an optional parameter
%             that is **only** used for plotting and visualization.
%
%   See also: firfilter, firwin, pfilt, filterbank

% Define initial value for flags and key/value pairs.
definput.import={'normalize'};
definput.importdefaults={'energy'};
definput.keyvals.delay=0;
definput.keyvals.fc=0;
definput.keyvals.fs=[];
definput.keyvals.scal=1;
definput.keyvals.min_win=1;
definput.flags.pedantic = {'nopedantic','pedantic'};
definput.flags.real={'complex','real'};

[flags,kv]=ltfatarghelper({'fc'},definput,varargin,'blfilter');

if flags.do_pedantic
    fc_offset = @(L,fc) L/2*fc-round(L/2*fc);
else
    fc_offset = @(L,fc) 0;
end

[fsupp,kv.fc,kv.delay,kv.scal]=scalardistribute(fsupp,kv.fc,kv.delay,kv.scal);

if ~isempty(kv.fs)
    fsupp=fsupp/kv.fs*2;
    kv.fc=kv.fc/kv.fs*2;
else
   %If fs is not specified, allow fsupp to be in range 0-2
   assert(all(fsupp>0) && all(fsupp<=2),...
          '%s: Filter support should be in range ]0-2].',...
          upper(mfilename));
end;

% Sanitize
kv.fc=modcent(kv.fc,2);

Nfilt=numel(fsupp);
gout=cell(1,Nfilt);

if ischar(winname)
    wn = {winname};
elseif iscell(winname)
    wn = winname;
else
    error('%s: Incorrect format of winname.',upper(mfilename));
end


for ii=1:Nfilt
    g=struct();
    
    if flags.do_1 || flags.do_area 
        g.H=@(L)    fftshift(firwin(wn{1},max([kv.min_win,...
                                    round(L/2*fsupp(ii))]),wn{2:end},...
                                    flags.norm,...
                                    'shift',-fc_offset(L,kv.fc(ii))))...
                                    *kv.scal(ii)*L;        
    end;
    
    if  flags.do_2 || flags.do_energy
        g.H=@(L)    fftshift(firwin(wn{1},max([kv.min_win,...
                                    round(L/2*fsupp(ii))]),wn{2:end},...
                                    flags.norm,...
                                    'shift',-fc_offset(L,kv.fc(ii))))...
                                    *kv.scal(ii)*sqrt(L);                
    end;
        
    if flags.do_inf || flags.do_peak
        g.H=@(L)    fftshift(firwin(wn{1},max([kv.min_win,...
                                    round(L/2*fsupp(ii))]),wn{2:end},...
                                    'shift',-fc_offset(L,kv.fc(ii))))...
                                    *kv.scal(ii);        
        
    end;
        
    g.foff=@(L) round(L/2*kv.fc(ii))-floor(max([kv.min_win,round(L/2*fsupp(ii))])/2);
    g.realonly=flags.do_real;
    g.delay=kv.delay(ii);
    g.fs=kv.fs;
    gout{ii}=g;
end;

if Nfilt==1
    gout=g;
end;

