function [g,tfr]=pgauss(L,varargin)
%PGAUSS  Sampled, periodized Gaussian
%   Usage: g=pgauss(L);
%          g=pgauss(L,tfr);
%          g=pgauss(L,...);
%          [g,tfr]=pgauss( ... );
% 
%   Input parameters:
%         L    : Length of vector.
%         tfr  : ratio between time and frequency support.
%
%   Output parameters:
%         g    : The periodized Gaussian.
%
%   `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/tfr))$.
%
%   The $l^2$ norm of the returned Gaussian is equal to 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*.
%
%   `[g,tfr] = pgauss( ... )` will additionally return the time-to-frequency
%   support ratio. This is useful if you did not specify it (i.e. used the
%   `'width'` or `'bw'` flag).
%
%   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)`.
%
%   In addition to the `'width'` flag, `pgauss` understands the following
%   flags at the end of the list of input parameters:
%
%     'fs',fs     Use a sampling rate of *fs* Hz as unit for specifying the
%                 width, bandwidth, centre frequency and delay of the
%                 Gaussian. Default is *fs=[]* which indicates to measure
%                 everything in samples.
%
%     'width',s   Set the width of the Gaussian such that it has an
%                 effective support of *s* samples. This means that
%                 approx. 96% of the energy or 79% of the area
%                 under the graph is contained within *s* samples. 
%                 This corresponds to -6dB or to width at the
%                 half of the height. 
%                 This is equivalent to calling `pgauss(L,pi*s^2/4L*log(2))`.
%
%     'atheight',ah  Used only in conjuction with 'width'. Forces the 
%                    Gaussian to width *s* at the *ah* fraction of the
%                    height.
%
%     'bw',bw     As for the `'width'` argument, but specifies the width
%                 in the frequency domain. The bandwidth is measured in 
%                 normalized frequencies, unless the `'fs'` value is given.
%
%     'cf',cf     Set the centre frequency of the Gaussian to *fc*.  
%
%     'wp'        Output is whole point even. This is the default.
%
%     'hp'        Output is half point even, as most Matlab filter
%                 routines.
%
%     'delay',d   Delay the output by *d*. Default is zero delay.
%
%   In addition to these parameteres, `pgauss` accepts any of the flags
%   from |normalize|. The output will be normalized as specified.
%
%   If this function is used to generate a window for a Gabor frame, then
%   the window giving the smallest frame bound ratio is generated by
%   `pgauss(L,a*M/L)`.
%
%   Examples:
%   ---------
%
%   This example creates a Gaussian function, and demonstrates that it is
%   its own Discrete Fourier Transform:::
%
%     g=pgauss(128);
%
%     % Test of DFT invariance: Should be close to zero.
%     norm(g-dft(g))
% 
%   The next plot shows the Gaussian in the time domain:::
%
%     plot(fftshift(pgauss(128)));
% 
%   The next plot shows the Gaussian in the frequency domain on a log scale:::
%
%     magresp(pgauss(128),'dynrange',100);
%     
%   The next plot shows the Gaussian in the time-frequency plane:::
% 
%     sgram(pgauss(128),'tc','nf','lin');
%
%   See also:  dgtlength, psech, firwin, pbspline, normalize
%
%   Demos:  demo_pgauss
%
%   References: mazh93

% AUTHOR : Peter L. Søndergaard.

%   First reference on this found in mazh93 eq. 63

if nargin<1
  error('Too few input parameters.');
end;

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;

% Define initial value for flags and key/value pairs.
definput.import={'normalize'};

definput.flags.centering={'wp','hp'};
definput.flags.delay={'nodelay','delay'};
definput.flags.width={'tfr','width','bw'};

definput.keyvals.tfr=1;
definput.keyvals.delay=0;
definput.keyvals.width=0;
definput.keyvals.fs=[];
definput.keyvals.cf=0;
definput.keyvals.bw=0;
definput.keyvals.atheight=[];

[flags,keyvals,tfr]=ltfatarghelper({'tfr'},definput,varargin);

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

if ~isempty(keyvals.atheight) && ~flags.do_width
    error(['%s: Param. ''atheight'' must be used together with param.',...
           ' ''width''. '],upper(mfilename));
end

if isempty(keyvals.atheight)
    keyvals.atheight = 0.5;
end

if keyvals.atheight >= 1 || keyvals.atheight <=0
    error('%s: Param. ''atheight'' must be in the range ]0,1[.',...
           upper(mfilename));
end


fs=keyvals.fs;

if flags.do_wp
  cent=0;
else
  cent=0.5;
end;

if isempty(fs)
  
  if flags.do_width
     tfr=pi/(4*log(1/keyvals.atheight))*keyvals.width^2/L; 
  end;
  
  if flags.do_bw
     tfr=L/(keyvals.bw*L/2)^2;
  end;
  
  delay_s=keyvals.delay;
  cf_s   =keyvals.cf;
else
  
  if flags.do_width
    tfr=(keyvals.width*fs)^2/L;
  end;

  if flags.do_bw
    tfr=L/(keyvals.bw/fs*L)^2;
  end;
  
  delay_s=keyvals.delay*fs;
  cf_s   =keyvals.cf/fs*L;
end;

g=comp_pgauss(L,tfr,cent-delay_s,cf_s);

g=normalize(g,flags.norm);

