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 SETNORM. 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, setnorm
%
% 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-2023 Peter L. Soendergaard <peter@sonderport.dk> and others.
% This file is part of LTFAT version 2.6.0
%
% 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={'setnorm'};
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=setnorm(g,flags.norm);