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:
%     S. Mallat and Z. Zhang. Matching pursuits with time-frequency
%     dictionaries. IEEE Trans. Signal Process., 41(12):3397-3415, 1993.
%     
%     
%
%   Url: http://ltfat.github.io/doc/fourier/pgauss.html

% Copyright (C) 2005-2016 Peter L. Soendergaard <peter@sonderport.dk>.
% This file is part of LTFAT version 2.1.2
%
% 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 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);


