function F=GetLOT(N,rxx);
% GetLOT    Lapped Orthogonal Transform of size 2NxN
%           This function returns the synthesis matrix for the lapped orthogonal 
% transform of size 2NxN, where we should have N even. The algorithm is based on
% the article: Malvar, 'The LOT: Transform Coding without Blocking Effects'
% in IEEE Trans on ASSP. vol 37 No 4 April 1989.
%
% F=GetLOT(N);  
% F=GetLOT(N,rxx);  
% ---------------------------------------------------------------------------% arguments:
%  F        The synthesis matrix for the lapped orthogonal transform, size is 2NxN
%  N        The size of F, should be even 
%  rxx      Optional argument for the signal autocorrelation function
%           this could be of size 2Nx1, where rxx(1) is for zero lag (rxx0)
%           or it could be of size 1x1, then it is assumed to be rho for an AR-1 process
%           if omitted the F matrix is not multiplied by unitary matriz Z.
% ---------------------------------------------------------------------------
%----------------------------------------------------------------------
% Copyright (c) 2000.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@tn.his.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:
% Ver. 1.0  21.07.2000  KS: m-file made
% Ver. 1.1  28.11.2002  KS: moved from ..\Frames to ..\FrameTools
%----------------------------------------------------------------------

Mfile='GetLOT';
SortFrequency=1;

% check input and output arguments, and assign values to arguments
if (nargin < 1); 
   error([Mfile,': function must have one input argument, see help.']); 
end
if (nargout ~= 1); 
   error([Mfile,': function must have one output arguments, see help.']); 
end
if mod(N,2)
   error([Mfile,': N is not even, see help.']); 
end
if (nargin >= 2); 
   [n,k]=size(rxx);
   if (n==1) & (k==1)
      rho=rxx;
      rxx=rho.^((0:(2*N-1))');
   end
   [n,k]=size(rxx);
   if ((n~=(2*N)) | (k~=1))
      error([Mfile,': illegal size of rxx, see help.']); 
   end
end

% make some needed matrices
D=idct(eye(N));        % the DCT functions
De=D(:,1:2:N);
Do=D(:,2:2:N);
D=De-Do;
J=flipud(eye(N));
F=0.5*[D,D;J*D,-J*D];

if (nargin >= 2); 
   Rxx=toeplitz(rxx);
   [Z,D]=eig(F'*Rxx*F);
   F=F*Z;
   if ((N==16) & (rxx(2)==0.95))
      % makes this the same as the results using GetLOT3216
      P=zeros(N);    % an orthogonal permutation matrix
      P(1,1)=1;P(3,2)=-1;P(2,3)=1;P(5,4)=-1;
      P(4,5)=-1;P(7,6)=-1;P(6,7)=-1;P(12,8)=1;
      P(8,9)=-1;P(13,10)=-1;P(9,11)=-1;P(16,12)=-1;
      P(10,13)=-1;P(14,14)=-1;P(11,15)=-1;P(15,16)=1;
      F=F*P;
   end
end

if SortFrequency
   % sort vectors in F by frequency
   Ff=fft(F);
   Ff=Ff(1:(N+1),:);
   Ff=abs(Ff.*Ff);
   temp=sum(Ff.*((1:(N+1))'*ones(1,N)));
   [t1,i1]=sort(temp);
   P=zeros(N); % an orthogonal permutation matrix
   for k=1:N;P(i1(k),k)=1;end;
   F=F*P;
end

return

% ---------------------------------------------------------------------------% check the results, is F orthogonal?
N=8;rxx=0.95;
F=GetLOT(N,rxx);
figure(1);clf;
PlotF(F);
disp(norm(F'*F-eye(N)));
W=zeros(2*N);
for k=(N+1):(2*N);
   W(k-N,k)=1;
end
disp(norm(F'*W*F));
disp('All the displayed numbers should be small.');

% display the synthesis vectors
G=reshape(1:(2*N*N),2*N,N);
G=G';
G=reshape(G,N,N,2);
for k=1:2; G(:,:,k)=G(:,:,k)'; end;
f=F(:);
figure(1);clf;
PlotF(BuildFg(f,G));   % this also illustrate the overlap