getelt.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function F=GetELT(N,p,rxx);
% GetELT    Extendend Lapped (orthogonal) Transform, ELT, of size 4NxN
%           This function returns the synthesis matrix for the extendend lapped 
% (orthogonal) transform of size 4NxN, where we should have N even. 
% The algorithm is based on the article: Henrique S. Malvar:
% 'Extendend Lapped Transforms: Fast algorithms and applications'
%
% F=GetELT(N);  
% F=GetELT(N,p);  
% F=GetELT(N,p,rxx);  
% ---------------------------------------------------------------------------
% arguments:
%  F        The synthesis matrix for the ELT, size is 4NxN
%  N        The size of F, should be even 
%  p        The free parameter in the design, 0<=p<=1, default p=0.7
%  rxx      Optional argument for the signal autocorrelation function
%           this could be of size 4Nx1, 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:  [e-mail]   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  24.07.2000  KS: m-file made
% Ver. 1.1  28.11.2002  KS: moved from ..\Frames to ..\FrameTools
%----------------------------------------------------------------------

Mfile='GetELT';
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 (nargin < 2); p=0.7; 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 >= 3); 
   [n,k]=size(rxx);
   if (n==1) & (k==1)
      rho=rxx;
      rxx=rho.^((0:(4*N-1))');
   end
   [n,k]=size(rxx);
   if ((n~=(4*N)) | (k~=1))
      error([Mfile,': illegal size of rxx, see help.']); 
   end
end

% make some needed values
i=0:(N-1);
thi=((1-p)/(2*N)*(2*i+1)+p).*((2*i+1)*pi/(8*N));
si=sin(thi);
ci=cos(thi);
h=zeros(2*N,1);
for i=0:(N/2-1);
   h(N/2-i)=-si(i+1)*si(N-i);
   h(N/2+i+1)=si(i+1)*ci(N-i);
   h(3*N/2-i)=ci(i+1)*si(N-i);
   h(3*N/2+i+1)=ci(i+1)*ci(N-i);
end
h=[h;flipud(h)];
 
F=zeros(4*N,N);
for k=1:N
   for n=1:(4*N)
      F(n,k)=h(n)*sqrt(2/N)*cos(pi/N*(k-0.5)*(n-1+(N+1)/2));
   end
end

if (nargin >= 3); 
   Rxx=toeplitz(rxx);
   [Z,D]=eig(F'*Rxx*F);
   F=F*Z;
end

if SortFrequency
   % sort vectors in F by frequency
   Ff=fft(F);
   Ff=Ff(1:(2*N+1),:);
   Ff=abs(Ff.*Ff);
   temp=sum(Ff.*((1:(2*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;p=0.7;rxx=0.95;
F=GetELT(N,p,rxx);
figure(1);clf; PlotF(F);
disp(norm(F'*F-eye(N)));
W=zeros(4*N);
for k=(N+1):(4*N);
   W(k-N,k)=1;
end
disp(norm(F'*W*F));
disp(norm(F'*W*W*F));
disp(norm(F'*W*W*W*F));
disp('All the displayed numbers should be small.');