cantight.m

Matlab时频分析工具箱,希望能对大家有所帮助啊

function gamma=cantight(p1,p2,p3,p4)
%CANTIGHT  Canonical tight window.
%   Usage:  gamma=cantight(a,M,L);
%           gamma=cantight(g,a,M);
%           gamma=cantight(g,a,M,L);
%
%   Input parameters:
%         g     : Gabor window.
%         a     : Length of time shift.
%         M     : Number of modulations.
%         L     : Length of window. (optional)
%   Output parameters:
%         gamma : Canonical tight window, column vector.
%
%   CANTIGHT(a,M,L) computes a nice tight IIR window of length L for a
%   lattice with parameters a, M.
%
%   CANTIGHT(g,a,M) computes the canonical tight window of the Gabor frame
%   with window g and parameters a, M.  It is assumed that g is an IIR
%   window, and the length of the returned tight window is the same as the
%   length of g.
%
%   If the length of g is equal to M, then the input window is assumed to
%   be a FIR window. In this case, the canonical tight window also has
%   length g.
%
%   CANTIGHT(g,a,M,L) pads or truncates g to a length L IIR window before
%   calculating the tight window. gamma will also have length L.
%
%   If a>M then an orthonormal window of the Gabor Riesz sequence with
%   window g and parameters a and M will be calculated.
%
%   SEE ALSO:  CANDUAL, FIR2IIR, DGT

% 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 Soendergaard.

wasrow=0;
dofirdual = 0;

error(nargchk(3,4,nargin));

if numel(p1)==1
  % First argument is a scalar.
    
  % Type 1: Nice tight win for sep system.

  a=p1;
  M=p2;
  L=p3;

  assert_squarelat(a,M,1,'CANTIGHT',0);
  [b,N,L]=assert_L(L,L,L,a,M,'CANTIGHT');

  % Create default window, a Gaussian.
  g=comp_pgauss(L,a/b,0);

else
  % First argument assumed to be a vector.
  
  g=p1;
  L=length(g); % might be overwritten later.
  
  % Take care of row/column preferences. Convert to column layout.
  if size(g,1)==1
    wasrow=1;
    % Make g a column vector.
    g=g(:);
  end;
  
  if numel(p2)==1
    % Second argument a scalar.
    
    if nargin<3
      error('Too few input parameters.');
    end;
    
    a=p2;
    M=p3;
    
    assert_squarelat(a,M,1,'CANTIGHT',0);

    Lwindow=length(g);
    Ls=Lwindow;

    if nargin<4
      [b,N,L]=assert_L(Ls,Lwindow,[],a,M,'CANTIGHT');
      
      if Ls<L
        % We are dealing with a FIR window.
        if Ls<=M
          % The dual will also be FIR
          dofirdual=1;
          g=fir2iir(g,L);
        else
          error('The computed tight window is an IIR window. Please specify the exact length L of the window: h=candual(g,a,M,L)');
        end;
      end;

    else
      L=p4;
      [b,N,L]=assert_L(L,Lwindow,L,a,M,'CANTIGHT');
      g=fir2iir(g,L);
    end;
       
  end;
end;

% If input is real, output must be real as well.
inputwasreal = isreal(g);

scale=1;
if a>M
  % Handle the Riesz basis (dual lattice) case.
  % Swap a and M, and scale differently.
  scale=sqrt(a/M);
  tmp=a;
  a=M;
  M=a;
end;

% Do the actual computation
gf=comp_wfac(g,a,M);
gammaf=comp_cantight_fac(gf,L,a,M)*scale;
gamma=comp_iwfac(gammaf,L,a,M);

if dofirdual
  gamma=iir2fir(gamma,Ls);
end;

if inputwasreal
  gamma=real(gamma);
end;
   
if wasrow
  gamma=gamma.';
end;