decom1d.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function W=Decom1D(Method, X, p1, p2)
% Decom1D   Decompose a 1D-signal into expansion coefficients, 
%           or the inverse function if Method is negative.
% Most decompositions methods suppurted by this function are orthogonal.
% Many of these methods also use other functions to do the real work.
% A special non-orthogonal case, frames, is also supported, the frame should
% then be stored in a mat-file with variables as specified in FrameFile.m
%
% Examples:
% W = Decom1D(Method, X);            % decomposition by Method
% X = Decom1D(-Method, W);           % inverse decomposition by Method
% W = Decom1D(Method, X, p1, p2);    % also using parameters p1, p2
% X = Decom1D(-Method, W, p1, p2);   % inverse, using parameters p1, p2
% ---------------------------------------------------------------------------
% arguments:
%  X        the signal, a column vector of real data, size Samplesx1
%  Method   An integer for the decomposition method 
%           Negative sign is used for the inverse function
%           1        - No decomposition 
%           2-128    - 2x2 to 128x128 Discrete Cosine Transform, DCT
%           201-207  - Tree structure IIR filter bank, 2^(Method-200) subbands
%                      Two additional arguments may be given, they are optional
%             p1       a0 used in the IIR filterbank.
%             p2       a1 used in the IIR filterbank.
%                      NOTE: energy of W will be 2^(Method-200) times energy of X
%           210      - 32x16 Lapped Orthogonal Transform, LOT. 
%                      (this is the same as Method 225)
%           211-217  - Dyadic filter bank, Daubechies 7-9 biorthogonal filters
%                      1-7 (Method-210) levels, 2-8 (Method-209) subbands
%           218      - wavelet decomposition (The wavelet toolbox is needed!)
%                      Two additional arguments should be given (default 'db6', 4 levels)
%                      db6 is the Daubechies wavlets, filterlength is 12.
%             p1       the wavelet name, as 'wname' used in 'wfilters' command
%             p2       The number of levels, 1-7.
%           220-229  - 2NxN Lapped Orthogonal Transform, LOT from GetLOT
%                      where we use N=[4,6,8,10,12,16,20,24,32,64].
%                      A third argument may be given, it is optional
%             p1       autocorrelation of signal, as 'rxx' used in 'GetLOT.m'
%           230-239  - 4NxN Extended Lapped Transform, ELT from GetELT
%                      where we use N=[4,6,8,10,12,16,20,24,32,64].
%             p1       free variable p, as 'p' used in 'GetELT.m'
%             p2       autocorrelation of signal, as 'rxx' used in 'GetELT.m'
%           255      - Use the frame stored in FrameFile
%             p1       the name of the mat-file for the frame, FrameFile
%             p2       the sparseness factor to use (override Savg in FrameFile)
%  W        the expansion coefficients, this matrix will be of size KxL
%           where K depends on the decomposition method, and L also 
%           depends on the number of samples, often we will have Samples=K*L
% ---------------------------------------------------------------------------

% This function needs the following functions to be available:
% C_anafb, C_synfb, (which again use C_ana and C_syn), C_anablock, C_synblock
% GetLOT, GetELT, and wfilters in Matlab Wavelet Toolbox.
% For frames: BuildFg, VSblock, GlobalMP  ??

%----------------------------------------------------------------------
% 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:  dd.mm.yyyy
% Ver. 1.0  17.07.2000  Karl Skretting, function made in the Frame context
% Ver. 1.1  20.12.2000  KS: some changes 
% Ver. 1.2  22.01.2001  KS: Decomposition using frames added
% Ver. 1.3  13.11.2001  KS: VSblock and VSolap1 possible for frames, Display added
% Ver. 1.4  02.12.2002  KS: moved from ..\Frames to ..\FrameTools
%----------------------------------------------------------------------

Mfile='Decom1D';
Display=1;

% check input and output arguments, and assign values to arguments
if (nargin < 4); p2=[]; end;
if (nargin < 3); p1=[]; end;
if (nargin < 2); 
   error([Mfile,': function must have two input arguments, see help.']); 
end
NofArgs=4;
if length(p2)==0; NofArgs=3; end;
if length(p1)==0; NofArgs=2; end;
if (nargout ~= 1); 
   error([Mfile,': function must have one output arguments, see help.']); 
end
if (Method>0); Decomp=1; else Decomp=0; end;
Method=abs(floor(Method));
if (Method==218) & (NofArgs<4)       % use a default
   p1='db6';p2=4;      
end

% find N and K and t2 (a text describing the method) from the Method
t1=[Mfile,': Method ',int2str(Method),' '];
if sum(Method==[1:128]);       
    N=Method;K=N;
    t2=[t1,int2str(N),'x',int2str(N),' Discrete Cosine Transform, DCT.'];
elseif sum(Method==[201:207]); 
    N=2^(Method-200);K=N;
    t2=[t1,'IIR tree structured filter bank with ',int2str(N),' subbands.'];
elseif (Method==210);          
    N=16;K=N;
    t2=[t1,'16x32 Lapped Orthogonal Transform.'];
elseif sum(Method==[211:217]); 
    N=2^(Method-210);K=N;
    t2=[t1,'Dyadic Daubechies 7-9 filter bank, ',int2str(Method-210),' levels.'];
elseif Method==218; 
    N=2^p2;K=N;
    t2=[t1,p1,'-wavelet  of ',int2str(p2),' levels.'];
elseif sum(Method==[220:229]); 
    if Method==220; N=4; end;
    if Method==221; N=6; end;
    if Method==222; N=8; end;
    if Method==223; N=10; end;
    if Method==224; N=12; end;
    if Method==225; N=16; end;
    if Method==226; N=20; end;
    if Method==227; N=24; end;
    if Method==228; N=32; end;
    if Method==229; N=64; end;
    K=N;
    t2=[t1,int2str(2*N),'x',int2str(N),' Lapped Orthogonal Transform.'];
elseif sum(Method==[230:239]); 
    if Method==230; N=4; end;
    if Method==231; N=6; end;
    if Method==232; N=8; end;
    if Method==233; N=10; end;
    if Method==234; N=12; end;
    if Method==235; N=16; end;
    if Method==236; N=20; end;
    if Method==237; N=24; end;
    if Method==238; N=32; end;
    if Method==239; N=64; end;
    K=N;
    t2=[t1,int2str(4*N),'x',int2str(N),' Extended Lapped Orthogonal Transform.'];
elseif (Method==255)
    FrameFile=p1;N=0;Mdim=0;
    if exist([FrameFile,'.mat'])
        try
            % load the following variables from FrameFile: Class, Type, Mdim, F, SizeF
            % G, Ctab, Dtab, Fbest, Savg, Mdat, PreProc, VecSel, InitialF, History, SNRtot
            load(FrameFile);
            K=SizeF(Mdim+1);
            N=prod(SizeF(1:Mdim));
            P=prod(SizeF((Mdim+2):(2*Mdim+1)));
            S=Savg;
            if (length(p2)); if (p2>0); S=p2; end; end;
            t2=[t1,'using FrameFile ',FrameFile,'.mat. Class=',...
               Class,'  Type=',Type','  Mdim=',int2str(Mdim),'  N=',...
               int2str(N),',  K=',int2str(K),',  P=',int2str(P),'.'];
        catch
            % just continue
        end
    end
    if (N==0) | (Mdim~=1)
        error([Mfile,': Not a frame for 1D-signal in ',FrameFile,'.mat']); 
    end
    if Type=='g'
        Fg=BuildFg(F,G);
    else
        Fg=F;
    end
else
    t2=[Mfile,': An undefined method given.'];
    W=[];
    disp(t2);
    return;
end

if Decomp
    % do decomposition
    disp([t2,' (decomposition).']);
    X=X(:);
    Samples=length(X);
    if (rem(Samples,N) > 0); 
        error([Mfile,': N (as given by Method) is not a factor of length of X.']); 
    end
    L=length(X)/N;
    xx2=X'*X;          
    %
    if (Method==1);   
        W=X;
    elseif sum(Method==[2:128]);       
        W=dct(reshape(X,N,L)); 
    elseif sum(Method==[201:207]); 
        if (NofArgs == 4)
            if (p1*p2)==0
                W=C_anafb(X,N).';    
            else
                W=C_anafb(X,N,p1,p2).';     % use parameters
            end
        else
            W=C_anafb(X,N).';    
        end
    elseif (Method==210);          
        T=GetLOT(16,0.95); 
        T=T';                         % the forward transform  size 16x32   
        W=reshape(X,N,L);             % W is 16xL, L=Samples/16 and N=16 
        W=[W;[W(:,2:L),W(:,1)]];      % W is 32xL  
        W=T*W;                        % W is 16xL
    elseif sum(Method==[211:217]); 
        % The Daubechies 7-9 biorthogonal wavelet is like in table 7.2 page 119
        % in C.S.Burrows et.al "Introduction to Wavelets and Wavelet Transforms"