function [g,tfr]=psech(L,p2,p3,p4)
%PSECH  Sampled, periodized hyperbolic secant
%   Usage: g=psech(L);
%          g=psech(L,tfr);
%          g=psech(L,s,'samples);
%          [g,tfr]=psech( ... );
%
%   Input parameters:
%      L   : Length of vector.
%      tfr : ratio between time and frequency support.
%   Output parameters:
%      g   : The periodized hyperbolic cosine.
%
%   `psech(L,tfr)` computes samples of a periodized hyperbolic secant.
%   The function returns a regular sampling of the periodization
%   of the function $sech(\pi\cdot x)
%
%   The returned function has norm 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*.
%
%   `psech(L)` does the same setting $tfr=1$.
%
%   `psech(L,s,'samples')` returns a hyperbolic secant with an effective
%   support of *s* samples. This means that approx. 96% of the energy or 74%
%   or the area under the graph is contained within *s* samples. This is
%   equivalent to `psech(L,s^2/L)`.
%
%   `[g,tfr] = psech( ... )` additionally returns the time-to-frequency
%   support ratio. This is useful if you did not specify it (i.e. used
%   the `'samples'` input format).
%
%   The function is whole-point even.  This implies that `fft(psech(L,tfr))`
%   is real for any *L* and *tfr*.
%
%   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
%   `psech(L,a*M/L)`.
%
%   Examples:
%   ---------
%
%   This example creates a `psech` function, and demonstrates that it is
%   its own Discrete Fourier Transform:::
%
%     g=psech(128);
%
%     % Test of DFT invariance: Should be close to zero.
%     norm(g-dft(g))
% 
%   The next plot shows the `psech` in the time domain compared to the Gaussian:::
%
%     plot((1:128)',fftshift(pgauss(128)),...
%          (1:128)',fftshift(psech(128)));
%     legend('pgauss','psech');
% 
%   The next plot shows the `psech` in the frequency domain on a log
%   scale compared to the Gaussian:::
%
%     hold all;
%     magresp(pgauss(128),'dynrange',100);
%     magresp(psech(128),'dynrange',100);
%     legend('pgauss','psech');
%     
%   The next plot shows `psech` in the time-frequency plane:::
% 
%     sgram(psech(128),'tc','nf','lin');
%
%   See also:  pgauss, pbspline, pherm
%
%   References: jast02-1

complainif_argnonotinrange(nargin,1,4,mfilename);

if nargin==1
  tfr=1;
end;

if size(L,1)>1 || size(L,2)>1
  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;

safe=12;

g=zeros(L,1);
sqrtl=sqrt(L);

w=tfr;

% Outside the interval [-safe,safe] then sech(pi*x) is numerically zero.
nk=ceil(safe/sqrt(L/sqrt(w)));

lr=(0:L-1).';
for k=-nk:nk  
  g=g+sech(pi*(lr/sqrtl-k*sqrtl)/sqrt(w));
end;

% Normalize it.
g=g*sqrt(pi/(2*sqrt(L*w)));
