function [g,tfr]=pgauss(L,p2,p3,p4)
%PGAUSS  Sampled, periodized Gaussian.
%   Usage: g=pgauss(L);
%          g=pgauss(L,tfr);
%          g=pgauss(L,s,'samples);
%          g=pgauss(L,tfr,cent);
%          g=pgauss(L,s,'samples,cent);
%          [g,tfr]=pgauss( ... );
% 
%   Input parameters:
%        L : Length of vector.
%      tfr : ratio between time and frequency support.
%     cent : Centering.
%   Output parameters:
%        g : The periodized Gaussian(s).
%
%   PGAUSS(L,tfr) computes samples of a periodized Gaussian.
%   The function returns a regular sampling of the periodization
%   of the function exp(-pi*(x.^2)).
%
%   The returned function has norm == 1.
%
%   The parameter tfr determines the ratio between the effective
%   support of g and the effective support of the DFT of g. If tfr>1 then
%   g has a wider support than the DFT of g.
%
%   PGAUSS(L) does the same setting tfr=1.
%
%   PGAUSS(L,s,'samples') returns a Gaussian with an effective support of s
%   samples. This means that approx. 96% of the energy or 79% or the area
%   under the graph is contained within s samples. This is equivalent to
%   PGAUSS(L,s^2/L);
%
%   [g,tfr] = pgauss( ... ) will additionally return the time-to-frequency
%   support ratio. This is usefull if you did not specify it (i.e. used
%   the 'samples' input format).
%
%   The function is whole-point even. This implies that FFT(PGAUSS(L,tfr))
%   is real for any L and tfr. The DFT of g is equal to PGAUSS(L,1/tfr).
%
%   PGAUSS(L,tft,cent) will generate a differently centered Gaussian.
%   Setting cent=0 generates a whole point even function (as above)
%   and setting cent=.5 generates a half point even function, as most
%   other Matlab filter routines.
%
%   If this function is used to generate a window for a Gabor frame, then
%   the window giving the smalles frame bound ratio is generated by
%   PGAUSS(L,a*M/L);
%
%   SEE ALSO:  IIRPAR, PSECH, FIRWIN, PBSPLINE
%
%   EXAMPLES:  EXAMP_PGAUSS
%
%   REFERENCES:
%     S. Mallat and Z. Zhang. Matching pursuits with time-frequency
%     dictionaries. IEEETrans.SignalProcess., 41(12):3397-3415, 1993.
%     

% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% Author : Peter Soendergaard.

%   First reference on this found in mazh93 eq. 63

error(nargchk(1,4,nargin));

if (prod(size(L,1))~=1 || ~isnumeric(L))
  error('L must be a scalar');
end;

if rem(L,1)~=0
  error('L must be an integer.')
end;

switch(nargin)
 case 1
  tfr=1;
  cent=0;
 case 2
  tfr=p2;
  cent=0;
 case 3
  if ischar(p3)
    switch(lower(p3))
     case {'s','samples'}
      tfr=p2^2/L;
     otherwise
      error('Unknown argument %s',p3);
    end;
    cent=0;
  else
    tfr=p2;
    cent=p3;
  end;
 case 4
  tfr=p2^2/L;
  cent=p4;
end;



if (prod(size(tfr,1))~=1 || ~isnumeric(tfr))
  error('tfr must be a scalar.');
end;

if (prod(size(cent,1))~=1 || ~isnumeric(cent))
  error('cent must be a scalar.');
end;

g=comp_pgauss(L,tfr,cent);

