comp_a.m

linear time－frequency toolbox

function [c,Ls]=comp_a(self,f,dorect,callfun,L)
%COMP_A  Compute Analysis transform
%   Usage: [c,Ls]=comp_a(self,f,dorect,callfun,L)
%
%   Input parameters:
%         tfr     : Time/frequency representation handle.
%         f       : Input data
%         dorect  : Do rectangular coefficient layout.
%         L       : Length of transform to do.
%   Output parameters:
%         c       : Array of coefficients.
%         Ls      : Length of input signal.
%

global TF_STORAGE;

% Verify that TF_STORAGE has been initialized
if prod(size(TF_STORAGE))==0
  error([callfun,' Please run TFR_CREATE to obtain a valid handle.']);
end;

% Change f to correct shape.
[f,Ls,W,wasrow,remembershape]=comp_sigreshape_pre(f,'DGT',0);

% Easy references
a=TF_STORAGE.data{self}.a;
M=TF_STORAGE.data{self}.M;
R=TF_STORAGE.data{self}.R;
dgtinfo=TF_STORAGE.data{self}.dgtinfo;

% Fix window and L
if L==-1
  [winpos,L]=assert_winready(self,-Ls,1,1,0,callfun);
else
  [winpos,L]=assert_winready(self,L,1,1,0,callfun);
end;

wininfo=TF_STORAGE.data{self}.w{winpos};

L_dgt=L*dgtinfo.Lmul;
a_dgt=dgtinfo.a;
M_dgt=dgtinfo.M;
N_dgt=L_dgt/a_dgt;

% Compute Gabor coefficients by either factorization or a FB
% algorithm.

f=postpad(f,L);    

%  ---------- compute the transform - preprocess signal

% Double f for the even boundary condition transform.
if dgtinfo.Lmul==2
  f=[f;flipud(f)];
end;

%  ----------- compute the transform - do the DGT --

% Number of signal times number of windows. These are treated together.
WR=W*R;

% Preprocess to handle c_f different from 0.
if (dgtinfo.c_f~=0)
  halfmod=exp(-2*pi*i*dgtinfo.c_f*(0:L_dgt-1).'/M_dgt);
  f=f.*repmat(halfmod,1,WR);
end;

% Do the actual DGT transform.
if abs(wininfo.Lwin)<L
  % Do the filter bank algorithm
  c=comp_dgt_fb(f,TF_STORAGE.data{self}.w{winpos}.g,a_dgt,M_dgt);

else  
  % Compute the window application and the DFT modulation.
  c=fft(comp_dgt_fw(f,TF_STORAGE.data{self}.w{winpos}.gf,a_dgt,M_dgt))/sqrt(M_dgt);

end;

if TF_STORAGE.data{self}.phaselock
  c=reshape(c,M_dgt,N_dgt,WR);
  c=phaselock(c,a_dgt);
  c=reshape(c,M_dgt,N_dgt*WR);
end;

% Post-process if c_t is different from 0.
if (dgtinfo.c_t~=0)

  % The following is necessary because REPMAT does not work for
  % 3D arrays.
  halfmod=repmat(exp(-2*pi*i*dgtinfo.c_t*((0:M_dgt-1)+dgtinfo.c_f).'/M_dgt),1,N_dgt*WR);
  
  c=c.*halfmod;
end;

% Reshape c to correct size.
c=reshape(c,M_dgt,N_dgt,WR);

% ------------ compute the transform - postprocess coefficients. --

% Postprocess if not doing a dgt.
if ~strcmp(dgtinfo.tt(1:3),'dgt')
  c=feval(['comp_',dgtinfo.tt],c,a);
end;

% ----------- do rectangular layout if needed.

if dorect

  switch dgtinfo.tt
    case {'dwilt','dwiltii'}
      Mout=M*2;
    otherwise
      Mout=M;
  end;
  N=M*L/a/Mout;
  
  if R>1
    c=reshape(c,Mout,N,R,W);
  else
    c=reshape(c,Mout,N,W);
  end;

else

  N=L/a;

  c=reshape(c,M*N*R,W);

end;