ttimes.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function [Y,varargout]=Ttimes(T,X,a3)
% Ttimes    Do the operation Y=T*X, where T represent a filterbank or a transform. 
%           For examples (and tests) of this function, see SignalExpansion.m
% using tests 2, 3 and 4.
%
% Y=Ttimes(T,X);             
% [Y,T]=Ttimes(T,X,a3);             
% ---------------------------------------------------------------------------
% arguments:
%  X     A NMx1 column vector or a NMxL matrix. 
%        If X is a matrix, the operation is applied to each column.
%  Y     the result, the size of Y will be KMxL (when K==N, same as X)
%  T   a If T is of size KxN, this is the transform y=Tx, where x and y
%        are parts of X and Y, with size Nx1 and Kx1
%      b If T is of size KxNxP, with P>1, T represent a filterbank which is
%        applied on each column of X to give a column of Y.
%        Delay: Note that if you give T as a NxNxP matrix (and K=N), and then use the
%        inverse of T (if it exist, U) again of size NxNxP (like LOT and ELT)
%        the resulting signal will be delayed N*(P-1) samples.
%        That is: Y=Ttimes(T,X); Xr=Ttimes(U,Y); then XR(n)==X(n-N*(P-1))
%        This delay will automatically be compensated if 
%        a3 is given as 1 (true or non-zero). If P==1 delay is no problem.
%      c If T is a single integer (size 1x1), we select transform or filterbank  
%        using the following table. A Negative number is used for the inverse (synthesis).
%          1    : 2x2 Haar wavelet, LP component is split until it has less than
%                 32 elements. (32 elements in LP is split again, 31 is not)
%                 this transform also handle any length in columns of X. 
%                 Use like: Y=Ttimes(1,X);Xr=Ttimes(-1,Y,0.5); 
%                 it is an integer transform, X is integer => Y is integer
%          2- 4 : may be used for other special transforms
%          5- 8 : NxN DCT transform, N=K=[2,4,8,16], P=1
%          9-12 : 2N*N LOT,          N=K=[2,4,8,16], P=2
%         13-16 : 4N*N ELT,          N=K=[2,4,8,16], P=4
%  a3    a third argument may be given. It may be used to indicate
%        - compensation for delay, when T is given as a matrix
%        - multiply result by a3, only for abs(T)==1, 2*2 Haar wavelet
%          (may be used as 1/sqrt(2) to keep transform unitary)
% ---------------------------------------------------------------------------

% ---------------------------------------------------------------------------
% Copyright (c) 1999.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  [e-mail]   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:
% Ver. 1.0  11.08.2000 Karl Skretting, function made 
% Ver. 1.1  27.11.2002  KS: moved from ..\Frames to ..\FrameTools
% ---------------------------------------------------------------------------

Mfile='Ttimes';
Display=0;

% check input and output arguments, and assign values to arguments
if (nargin < 2); 
   error([Mfile,': function must have two input arguments, see help.']); 
end
if (nargout < 1); 
   error([Mfile,': function must have one output arguments, see help.']); 
end
if (nargin < 3); a3=0; end;
if nargout==2; varargout(1)={T}; end;

[N,K,P]=size(T);
if ((K==1) & (P==1))
   Ta=abs(T);
   Ts=sign(T);
   % first handle the spesial case when Ta==1;
   if Ta==1
      T=[1,1;1,-1];
      [M,L]=size(X);M0=M; 
      if Ts==1
         if rem(M,2); X=[X;X(M,:)]; M=M+1; end;   % add one row
         Y=reshape(T*reshape(X,2,(M*L/2)),M,L);
         Y=[Y(1:2:M,:);Y(2:2:M,:)];
         if (M>31); Y(1:(M/2),:)=Ttimes(1,Y(1:(M/2),:),a3); end;
         if (M0<M); Y=Y(1:M0,:); end;
      else
         if rem(M,2); X=[X;zeros(1,L)]; M=M+1; end;   % add one row
         if (M>31); X(1:(M/2),:)=Ttimes(-1,X(1:(M/2),:),a3); end;
         Y=zeros(size(X));Y(1:2:M,:)=X(1:(M/2),:);Y(2:2:M,:)=X((M/2+1):M,:);
         Y=reshape(T*reshape(Y,2,(M*L/2)),M,L);
         if (M0<M); Y=Y(1:M0,:); end;
      end
      if a3; Y=Y*a3; end;   
      return;
   end
   if Ta<5
      if Display; disp([Mfile,': this transform is not ready yet. We set Y=X.']); end;
      Y=X; return;
   end
   % now we should find T matrix, first it will be the synthesis vectors (PNxK) 
   N=2^(mod(Ta-1,4)+1);
   if sum(Ta==[5:8]); T=idct(eye(N)); P=1; K=N;          
   elseif sum(Ta==[9:12]); T=GetLOT(N,0.95); P=2; K=N;
   elseif sum(Ta==[13:16]); T=GetELT(N,0.7,0.95); P=4; K=N; 
   else
      if Display; disp([Mfile,': this transform is not ready yet. We set Y=X.']); end;
      Y=X; return;
   end;
   % since unitary FBs, analysis = (synthesis)' and reshape
   T=reshape(T',K,N,P);  
   % transpose once again if synthesis
   if (Ts==(-1)); for p=1:P; T(:,:,p)=T(:,:,p)'; end; end;   
else
   % T matrix is given
   if a3; Ts=(-1); else Ts=1; end;
end

% here we have: [K,N,P]=size(T);   (assume analysis filter bank)
[MN,L]=size(X);
if rem(MN,N)
   if Display; disp([Mfile,': size of X is not comparable to size of T. We set Y=X.']); end;
   Y=X;  return;
end
M=MN/N;
if Display
   disp([Mfile,': Operator (',int2str(K),'x',int2str(N),'x',int2str(P),...
      ') operates on input arg. (',int2str(N*M),'x',int2str(L),...
      ') to give output arg. (',int2str(M*K),'x',int2str(L),').']);
end

X=reshape(X,N,M*L);
Y=T(:,:,1)*X;
for p=2:P
   X=reshape(X,N*M,L);
   if (Ts==1)
      X=X([(N+1):(N*M),1:N],:);   % shift up (left), first part last
   else
      X=X([((M-1)*N+1):(N*M),1:((M-1)*N)],:);   % shift down (right), last part first
   end
   X=reshape(X,N,M*L);
   Y=Y+T(:,:,p)*X;
end
Y=reshape(Y,K*M,L);

if nargout==2; varargout(1)={T}; end;
return;