psech.m

linear time－frequency toolbox

function [g]=psech(L,tfr)
%PSECH  Sampled, periodized hyperbolic secant.
%   Usage: g=psech(L);
%          g=psech(L,tfr);
% 
%   Input parameters:
%        L : Length of vector.
%      tfr : ratio between time and frequency support.
%   Output parameters:
%        g : The periodized Gaussian(s).
%
%   PSECH(L,tfr) computes samples of a periodized hyperbolic secant.
%   The function returns a regular sampling of the periodization
%   of the function sech(pi*x)
%
%   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 FFT of g.
%
%   PSECH(L) does the same setting tfr=1.
%
%   If tfr is a vector, PSECH will return a matrix, where each column
%   is a PSECH function with the corresponding value of tfr.
%
%   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 with
%   a rectangular lattice, then the optimal window is generated by
%   g=psech(L,a/b);
%
%   SEE ALSO:  PGAUSS, PBSPLINE, PHERM
%
%   REFERENCES:
%     A. J. E. M. Janssen and T. Strohmer. Hyperbolic secants yield Gabor frames.
%     Appl. Comput. Harmon. Anal., 12(2):259-267, 2002.

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

if nargin>2
  error('Too many input parameters.');
end;

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;

% Parse tfr and order.
if sum(1-(size(tfr)==1))>1
  error('tfr must be a scalar or vector');
end;

safe=12;

R=length(tfr);
g=zeros(L,R);
sqrtl=sqrt(L);

for r=1:R

  w=tfr(r);

  % 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(:,r)=g(:,r)+sech(pi*(lr/sqrtl-k*sqrtl)/sqrt(w));
  end;
  
  % Normalize it.
  g(:,r)=g(:,r)*sqrt(pi/(2*sqrt(L*w)));
  
end;