function Y=SignalExpansion(TestNo,Signal,N,p1,p2)
% SignalExpansion   Here we test many different methods for signal expansion. 
%                   A signal expansion is the representation of a signal as
% a linear combination of some basic synthesis signals, which depend on the
% transform, filter bank or frame used.
% The functions tested in this file are:
%   dct function in Matlab (for comparision) in TestNo 1
%   Ttimes function in TestNo 2, 3, and 4
%   Decom1D function in TestNo 5
%   Decom2D function in TestNo 6
%
% Y=SignalExpansion(TestNo,Signal,N,p1,p2);   
%-----------------------------------------------------------------------------------
% arguments:
%  Y          - the expansion coefficients
%               for 1D signal the size of Y depends on the expansion
%               for 2D signals (and orthogonal expansions) the size of Y is as size of X 
%  TestNo     - which test to do
%                1 - use NxN DCT (Matlab functions only)
%                    Y=SignalExpansion(1,1,8);  % 8x8 DCT, 1D AR(1) signal
%                    Y=SignalExpansion(1,2,8);  % 8x8 DCT (dct function), 2D image
%                    Y=SignalExpansion(1,2,8,'dct2');  % 8x8 DCT (dct2 function), 2D image
%                2 - use NxN KLT (Karhunen Loewe Transform) and Ttimes function (case a)
%                    Y=SignalExpansion(2,1,8);  % 8x8 KLT, 1D AR(1) signal
%                    Y=SignalExpansion(2,2,8);  % 8x8 KLT, two KLTs (rows and columns) 
%                    Y=SignalExpansion(2,2,8,1);  % 8x8 KLT, the same KLT for rows and columns
%                3 - use 4NxN ELT (Extended Lapped Transform) and Ttimes function (case b)
%                    Y=SignalExpansion(3,1,8);  % 32x8 ELT, 1D AR(1) signal
%                    Y=SignalExpansion(3,2,8);  % 32x8 ELT, image
%                4 - use a defined method in Ttimes function (case c), note that argument N 
%                    is not used and argument p1 gives the method used in Ttimes
%                    Y=SignalExpansion(4,1,0,1);  % 2x2 Haar wavelet
%                    Y=SignalExpansion(4,1,0,7);  % 8x8 DCT
%                    Y=SignalExpansion(4,2,0,12);  % 32x16 LOT
%                5 - use the methods in Decom1D, only for 1D signal,
%                    note the third argument, N, now is Method used in Decom1D, 
%                    Y=SignalExpansion(5,1,8);  % Decom1D, 8x8 DCT
%                    Y=SignalExpansion(5,1,203);  % Decom1D, IIR filter bank
%                    Y=SignalExpansion(5,1,213);  % Decom1D, Daubechies 7-9 wavelet
%                    Y=SignalExpansion(5,1,218,'db3',3);  % Decom1D, filter from wfilters
%                    Y=SignalExpansion(5,1,235);  % Decom1D, 64x16 ELT
%                    Y=SignalExpansion(5,1,255,'FrameEx2s20',0.25);  % Decom1D, frame
%                6 - use the methods in Decom2D, only for 2D signal, Decom2D is quite similar to
%                    Ttimes and it uses the Ttimes function to do the work.
%                    Y=SignalExpansion(6,2,8);   % Decom2D, 16x16 DCT
%                    Y=SignalExpansion(6,2,12);  % Decom2D, 32x16 LOT
%                    Y=SignalExpansion(6,2,2,'db4',3);   % Decom2D, a wavelet
%                    Y=SignalExpansion(6,2,2,'db79',4);   % Decom2D, a wavelet
%  Signal     -  1 - a simple one-dimensional signal, AR(1)
%                2 - a two-dimensional signal, the image Lena (or Barbara)
%                X - the signal X, 1D or 2D
%  N          - the transform (filter bank) size parameter N, default is 8
%  p1         - an extra argument used by some tests, default 0
%-----------------------------------------------------------------------------------

%----------------------------------------------------------------------
% Copyright (c) 2002.  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:  dd.mm.yyyy
% Ver. 1.0  28.11.2002  KS: function made 
%----------------------------------------------------------------------

Mfile='SignalExpansion';
Message='ok';
Y=0;

if nargin<3; N=8; end;
if nargin<4; p1=0; end;
if nargin<5; p2=0; end;
if nargin<2;
    Message=[Mfile,': wrong number of arguments, see help.'];
    disp(Message);
    return
end;

[Mi,Ni]=size(Signal);
if ((Mi*Ni)>1)     % the signal is given
    if ((Mi==1) | (Ni==1))
        X=Signal(:);
        Signal=1;     % now Signal tells number of dimensions in signal X
        L=length(X);
    else
        X=Signal;
        Signal=2;     % now Signal tells number of dimensions in signal X
    end
else
    if (Signal==1)
        L=2560;
        randn('state',6502);
        X=filter(1,[1,-0.95],randn(L,1));   % AR(1) signal
        X=X(:);
    elseif (Signal==2)
        if exist('GetIm.m')==2  % GetIm is my own function for loading images
            X=GetIm(1);     % load an image (Lena) 
            % X=GetIm(6,2);     % load an image (Barbara) 
        else
            % here you should put the code that load a grayscale image
            disp([Mfile,', is not able to find an image to load, edit mfile (line 98)']);
            return
        end
        X=double(X);
        [Mi,Ni]=size(X);
    else
        Message=[Mfile,': the signal is not (correctly) given, see help.'];
        disp(Message);
        return
    end
end

% ********** now do test 1 ****************************************************** 1 **

if (TestNo==1)  % NxN DCT
    disp(['Test ',int2str(TestNo),' using ',int2str(N),'x',int2str(N),' DCT.']);
    if (Signal==1)
        Y=dct(reshape(X,N,L/N));    % expansion coefficients
        Xr=idct(Y);
        Xr=Xr(:);
        temp=norm(X-Xr);
        disp([Mfile,', Test ',int2str(TestNo),' : Norm of error (X-Xr) is ',num2str(temp)]);
        figure(1);clf;
        subplot(2,1,1);plot(1:L,X);title('Original signal, X');
        subplot(2,1,2);plot(1:L,Xr);title('Reconstructed signal, Xr');
        if exist('geomean.m')==2
            sigma2=std(Y').^2;  % estimate for variance of each of the rows of Y
            temp=mean(sigma2)/geomean(sigma2);
            disp([Mfile,', Test ',int2str(TestNo),' : Estimate for coding gain is ',num2str(temp)]);
        end
    elseif (Signal==2)
        % the matlab function dct2 does not take dct of NxN blocks of the image
        if strcmp(p1,'dct2')  % use dct2  (time is 15.2 s)
            disp([Mfile,', Test ',int2str(TestNo),' : use dct2 function.']);
            Y=zeros(size(X));
            for i=1:N:Mi
                ii=i:(i+N-1);
                for j=1:N:Ni
                    jj=j:(j+N-1);
                    Y(ii,jj)=dct2(X(ii,jj));
                end
            end
            % now do the inverse, note that the inverse from the other clause could be used
            Xr=zeros(size(Y));
            for i=1:N:Mi
                ii=i:(i+N-1);
                for j=1:N:Ni
                    jj=j:(j+N-1);
                    Xr(ii,jj)=idct2(Y(ii,jj));
                end
            end
        else  % use reshape and dct (time is 3.5 s)
            disp([Mfile,', Test ',int2str(TestNo),' : use dct function.']);
            temp=dct(reshape(X,N,(Mi*Ni)/N));       % dct of columns
            temp=reshape(temp,Mi,Ni)';              % reshaped back to image size and transposed  
            temp=dct(reshape(temp,N,(Mi*Ni)/N));    % dct of rows
            Y=reshape(temp,Ni,Mi)';                 % reshaped back to image size and transposed
            % now do the inverse
            temp=reshape(Y',N,(Mi*Ni)/N);
            temp=idct(temp);
            temp=reshape(temp,Ni,Mi)';
            temp=idct(reshape(temp,N,(Mi*Ni)/N));
            Xr=reshape(temp,Mi,Ni);
        end
        %
        temp=norm(X(:)-Xr(:));
        disp([Mfile,', Test ',int2str(TestNo),' : Norm of error (X-Xr) is ',num2str(temp)]);
        figure(1);clf;
        subplot(1,2,1);imagesc(X);title('Original image, X');
        subplot(1,2,2);imagesc(Xr);title('Reconstructed image, Xr');
        if exist('geomean.m')==2
            % each NxN block of Y is made into a vector
            temp=Reorder(Y,[Mi,Ni],[N,N],1);
            sigma2=std(temp').^2;  % estimate for variance of each of the rows of Y
            temp=mean(sigma2)/geomean(sigma2);
            disp([Mfile,', Test ',int2str(TestNo),' : Estimate for coding gain is ',num2str(temp)]);
        end
    else
        Message=[Mfile,': Signal has illegal value in test ',int2str(TestNo)];
        disp(Message);
        return
    end
    Message=[Mfile,': Test ',int2str(TestNo),' finished ok.'];
end

% ********** now do test 2 ****************************************************** 2 **

if (TestNo==2)   % NxN KLT
    disp(['Test ',int2str(TestNo),' using ',int2str(N),'x',int2str(N),' KLT.']);
    if (Signal==1)
        Xc=reshape(X,N,L/N); 
        Rxx=Xc*Xc';
        [U,D]=eig(Rxx);    % the columns of U are the synthesis vectors
        T=U';              % the analysis part
        Y=Ttimes(T,X);     % Y is returned same size as X
        Xr=Ttimes(U,Y); 
        Xr=Xr(:);
        temp=norm(X-Xr);
        disp([Mfile,', Test ',int2str(TestNo),' : Norm of error (X-Xr) is ',num2str(temp)]);
        figure(1);clf;
        subplot(2,1,1);plot(1:L,X);title('Original signal, X');
        subplot(2,1,2);plot(1:L,Xr);title('Reconstructed signal, Xr');
        if exist('geomean.m')==2
            temp=reshape(Y,N,L/N);
            sigma2=std(temp').^2;  % estimate for variance of each of the rows of Y
            temp=mean(sigma2)/geomean(sigma2);
            disp([Mfile,', Test ',int2str(TestNo),' : Estimate for coding gain is ',num2str(temp)]);
        end
    elseif (Signal==2)
        if (p1==1)          % columns and rows used together to get a separable KLT transform
            disp([Mfile,', Test ',int2str(TestNo),' : use 1 KLT tramsform.']);
            Xc=[reshape(X,N,(Mi*Ni)/N),reshape(X',N,(Mi*Ni)/N)]; 
            Rxx=Xc*Xc';
            clear Xc;
            [U,D]=eig(Rxx);    % the columns of U are the synthesis vectors
            T=U';              % the analysis part
            temp=Ttimes(T,X);      % do the colums
            Y=Ttimes(T,temp')';    % and the rows 
            temp=Ttimes(U,Y')'; 
            Xr=Ttimes(U,temp); 
        else         % one KLT for the rows and another for the columns
            disp([Mfile,', Test ',int2str(TestNo),' : use 2 KLT tramsforms (rows + columns).']);
            Xc=reshape(X,N,(Mi*Ni)/N);   % columns
            Rxx=Xc*Xc';