gfdualnorm.m

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

function [o1,o2]=gfdualnorm(gamma,g,a,M,L);
%GFDUALNORM  Measure of how close a window is to be a dual window. 
%   Usage:  dn=gfdualnorm(g,gamma,a,M);
%           dn=gfdualnorm(g,gamma,a,M,L);
%           [scal,res]=gfdualnorm(g,gamma,a,M);
%           [scal,res]=gfdualnorm(g,gamma,a,M,L);
%
%   Input parameters:
%         gamma  : input window..
%         g      : window function.
%         a      : Length of time shift.
%         M      : Number of modulations.
%         L      : Length of transform to consider
%   Output parameters:
%         dn     : dual norm.
%         scal   : Scaling factor
%         res    : Residual
%
%   GFDUALNORM(g,gamma,a,M) calculates how close gamma is to be a dual
%   window of the Gabor frame with window g and parameters a and M.
%
%   GFDUALNORM(g,gamma,a,M,L) does the same, but considers a transform
%   length of L.
%
%   [scal,res]=GFDUALNORM(g,gamma,a,M) or
%   [scal,res]=GFDUALNORM(g,gamma,a,M,L) will compute two entities:
%   scal determines if the windows are scaled correctly, it must be 1 
%   for the windows to be dual. res is close to zero if the windows
%   (scaled correctly) are dual windows.
%
%   GFDUALNORM can be used to get the maximum relative reconstruction
%   error when using the two specified windows. Consider the following code
%   for some signal f, windows g, gamma, parameters a and M and 
%   transform-length L (See help on DGT for how to obtain L):
% 
%      fr=idgt(dgt(f,g,a,M),gamma,a); 
%      er=norm(f-fr)/norm(f);
%      eest=gfdualnorm(g,gamma,a,M,L);
% 
%   Then  er < eest for all possible input signals f.
%
%   To get a similar estimate for a tight window gt, simply use
%   
%      eest=gfdualnorm(gt,gt,a,M,L);
%
%   SEE ALSO:  DGT, CANTIGHT, CANDUAL.

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

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

assert_squarelat(a,M,1,'GFDUALNORM',0);

if size(g,2)>1
  if size(g,1)>1
    error('g must be a vector');
  else
    % g was a row vector.
    g=g(:);
  end;
end;

if size(gamma,2)>1
  if size(gamma,1)>1
    error('gamma must be a vector');
  else
    % gamma was a row vector.
    gamma=gamma(:);
  end;
end;
 
Ls=length(g);
Lwindow=length(gamma);

if nargin<5
  [b,N,L]=assert_L(Ls,Lwindow,[],a,M,'GFDUALNORM');
else
  [b,N,L]=assert_L(Ls,Lwindow,L,a,M,'GFDUALNORM');
  g=fir2iir(g,L);
  gamma=fir2iir(gamma,L);
end;

% Handle the Riesz basis (dual lattice) case.
if a>M

  % Calculate the right-hand side of the Wexler-Raz equations.
  rhs=dgt(gamma,g,a,M);
  scalconst=1;
  
else
  
  % Calculate the right-hand side of the Wexler-Raz equations.
  rhs=dgt(gamma,g,M,a);
  
  scalconst=a/M;
  
end;

if nargout<2
  % Subtract from the first element to make it zero, if the windows are
  % dual.
  rhs(1)=rhs(1)-scalconst;

  o1=norm(rhs(:),1);
else
  % Scale the first element to make it one, if the windows are dual.
  o1=rhs(1)/scalconst;
  o2=norm(rhs(2:end),1);
end;