decom1d.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

        rh=[-0.06453888262876,-0.04068941760920,0.41809227322204,...
            0.78848561640637,0.41809227322204,-0.04068941760920,-0.06453888262876];
        rg=[0.03782845550726,0.02384946501956,-0.11062440441844,-0.37740285561283,...
            0.85269867900889,-0.37740285561283,-0.11062440441844,0.02384946501956,...
            0.03782845550726];
        W=zeros(N,L);                 % N=2^(Method-210)
        x3=C_anablock(X,rh,rg);       % x3 is just used temporary here
        i=Method-211;                 % levels left
        while i
            W(((2^i)+1):(2^(i+1)),:)=reshape(x3(:,2),2^i,L);
            x3=C_anablock(x3(:,1),rh,rg);   
            i=i-1;
        end
        W(1:2,:)=x3';
    elseif (Method==218);             % wavelet
        [h,g,rh,rg] = wfilters(p1);      % building blocks of the filterbank
        W=zeros(N,L);                    % N=2^(Method-210)
        x3=C_anablock(X,h,g);              % x3 is just used temporary here
        i=p2-1;                          % levels left
        while i
            W(((2^i)+1):(2^(i+1)),:)=reshape(x3(:,2),2^i,L);
            x3=C_anablock(x3(:,1),h,g);   
            i=i-1;
        end
        W(1:2,:)=x3';
    elseif sum(Method==[220:229]);    % LOT
        if ((NofArgs>=3) & (p1~=0))
            T=GetLOT(N,p1);
        else
            T=GetLOT(N,0.95);
        end
        T=T';                         % the forward transform  size Nx(2N)   
        W=reshape(X,N,L);             % W is NxL, L=Samples/N  
        W=[W;[W(:,2:L),W(:,1)]];      % W is (2N)xL  
        W=T*W;                        % W is NxL
    elseif sum(Method==[230:239]);    % ELT
        if ((NofArgs==3) & (p1~=0))
            T=GetELT(N,p1,0.95);
        elseif ((NofArgs==4) & ((p1*p2)~=0))
            T=GetELT(N,p1,p2);         % 4NxN
        else
            T=GetELT(N,0.7,0.95);
        end
        T=T';                         % the forward transform  size Nx(4N)   
        W=reshape(X,N,L);             % W is NxL, L=Samples/N  
        W=[W;[W(:,2:L),W(:,1)];[W(:,3:L),W(:,1:2)];[W(:,4:L),W(:,1:3)]];  % W is (4N)xL  
        W=T*W;                        % W is NxL
    elseif Method==255
        % S is Savg stored in the frame or fourth argument (p2) if given
        X=reshape(X,N,L);
        % the variables for the frame is already loaded
        if strcmp(VecSel.Prog1,'VSblock') | strcmp(VecSel.Prog1,'VSolap1')
            % the distribution S(l) is found by VecSel.Prog1, since S is sparseness
            % factor it will initially be evenly distributed, but updated during
            % vector selection. This depends on the arguments B, el, and Bf
            VSalg=VecSel.arg1;
            B=VecSel.arg2;
            el=VecSel.arg3;
            Bf=VecSel.arg4;
            if el<1; el=1; end;
            % we do VSfomp2 first time anyway, (GMP+fomp2)
            if strcmp(VecSel.Prog1,'VSblock')
                W=VSblock(X,Fg,S,'VSfomp2',B,0,Bf);   
                R=X-Fg*W;
            else
                W=VSolap1(X,Fg,S,'VSfomp2',B,0);   
                R=X; for p=1:P; R=R-Fg(:,:,p)*W(:,[(L+2-p):L,1:(L+1-p)]); end;
            end   
            rr2=R(:)'*R(:);
            temp=10*log10(xx2/rr2);
            if Display; 
                disp([Mfile,': after VSfomp2,  SNR=',num2str(temp,4)]);
                for i=1:el
                    if strcmp(VecSel.Prog1,'VSblock')
                        W=VSblock(X,Fg,W,VSalg,B,0,Bf);   
                        R=X-Fg*W;
                    else
                        W=VSolap1(X,Fg,W,VSalg,B,0);   
                        R=X; for p=1:P; R=R-Fg(:,:,p)*W(:,[(L+2-p):L,1:(L+1-p)]); end;
                    end
                    rr2=R(:)'*R(:);
                    temp=10*log10(xx2/rr2);
                    Display; disp([Mfile,': Iteration ',int2str(i),' using ',VSalg,...
                            ',  SNR=',num2str(temp,4)]); 
                end
            else
                if strcmp(VecSel.Prog1,'VSblock')
                    W=VSblock(X,Fg,W,VSalg,B,el-1,Bf);   
                else
                    W=VSolap1(X,Fg,W,VSalg,B,el-1);   
                end
            end
            % S=full(sum(W~=0));
        else
            disp([Mfile,': not ready for vector selection using ',VecSel.Prog1,'.']);
            W=zeros(K,L);
        end
        W=full(reshape(W,K,L));         % W is size KxL
    else
        error([Mfile,': Unknown value of Method, see help.']);
    end
else
    disp([t2,' (reconstruction).']);
    % do inverse decomposition
    W=X(:);
    if (rem(length(W),K) > 0); 
        error([Mfile,': K (as given by Method) is not a factor of length of W.']); 
    end
    L=length(W)/K;
    W=reshape(W,K,L);
    if (Method==1);   
        X=W(:);
    elseif sum(Method==[2:128]);       
        X=idct(W); 
        X=X(:);
    elseif sum(Method==[201:207]); 
        if (NofArgs == 4)
            if (p1*p2)==0
                X=C_synfb(W.');            
            else
                X=C_synfb(W.',p1,p2).';     % use parameters
            end
        else
            X=C_synfb(W.');             % or not
        end
    elseif (Method==210);          
        T=GetLOT(16,0.95); 
        X=T*W;
        X=X(1:N,:)+[X((N+1):(2*N),L),X((N+1):(2*N),1:(L-1))];
        X=X(:);
    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"   
        h=[0.03782845550726,-0.02384946501956,-0.11062440441844,0.37740285561283,...
         0.85269867900889,0.37740285561283,-0.11062440441844,-0.02384946501956,...
         0.03782845550726];
        g=[0.06453888262876,-0.04068941760920,-0.41809227322204,0.78848561640637,...
         -0.41809227322204,-0.04068941760920,0.06453888262876];
        W=W.';
        X=W(:,1:2);
        for i=1:(Method-211)
            X=[C_synblock(X,h,g),reshape(W(:,(2^i+1):(2^(i+1)))',L*2^i,1)];
        end
        X=C_synblock(X,h,g);
        W=W.';
    elseif (Method==218);             % wavelet
        if NofArgs<4       % use a default
            p1='db6';p2=4;
        end
        [h,g,rh,rg] = wfilters(p1);      % building blocks of the filterbank
        W=W.';
        X=W(:,1:2);
        for i=1:(p2-1)
            X=[C_synblock(X,rh,rg),reshape(W(:,(2^i+1):(2^(i+1)))',L*2^i,1)];
        end
        X=C_synblock(X,rh,rg);
        W=W.';
    elseif sum(Method==[220:229]);    % LOT
        if ((NofArgs>=3) & (p1~=0))
            T=GetLOT(N,p1);
        else
            T=GetLOT(N,0.95);
        end
        X=T*W;
        X=X(1:N,:)+[X((N+1):(2*N),L),X((N+1):(2*N),1:(L-1))];
        X=X(:);
    elseif sum(Method==[230:239]);    % ELT
        if ((NofArgs==3) & (p1~=0))
            T=GetELT(N,p1,0.95);
        elseif ((NofArgs==4) & ((p1*p2)~=0))
            T=GetELT(N,p1,p2);         % 4NxN
        else
            T=GetELT(N,0.7,0.95);
        end
        X=T*W;                        % 4NxL
        X=X(1:N,:)+[X((N+1):(2*N),L),X((N+1):(2*N),1:(L-1))]+...
            [X((2*N+1):(3*N),(L-1):L),X((2*N+1):(3*N),1:(L-2))]+...
            [X((3*N+1):(4*N),(L-2):L),X((3*N+1):(4*N),1:(L-3))];     % NxL
        X=X(:);
    elseif Method==255
        if P==1
            X=full(Fg*W);
        else
            temp=W;X=zeros(N,L);
            for p=1:P
                X=X+full(Fg(:,:,p)*temp);   
                temp=[temp(:,L),temp(:,1:(L-1))];
            end
        end
        X=X(:);
    else
        error([Mfile,': Unknown value of Method, see help.']);
    end
    W=X(:);   % output argument name is W
end

return