signalexpansion.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

            [U1,D]=eig(Rxx);    % the columns of U are the synthesis vectors
            T1=U1';              % the analysis part
            Xc=reshape(X',N,(Mi*Ni)/N);  % rows
            Rxx=Xc*Xc';
            [U2,D]=eig(Rxx);    % the columns of U are the synthesis vectors
            T2=U2';              % the analysis part
            clear Xc;
            temp=Ttimes(T1,X);      % do the colums
            Y=Ttimes(T2,temp')';    % and the rows 
            temp=Ttimes(U2,Y')'; 
            Xr=Ttimes(U1,temp); 
        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 3 ****************************************************** 3 **

if (TestNo==3)   % 4NxN ELT
    disp(['Test ',int2str(TestNo),' using ',int2str(4*N),'x',int2str(N),' ELT.']);
    F=GetELT(N,0.7,0.95);
    T=zeros(N,N,4);Ti=zeros(N,N,4);
    for i=1:4
        T(:,:,i)=F((1:N)+(i-1)*N,:)';     % anlysis part
        Ti(:,:,i)=F((1:N)+(i-1)*N,:);     % synthesis part
    end
    if (Signal==1)
        Y=Ttimes(T,X);
        Xr=Ttimes(Ti,Y,1);
        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)
        temp=Ttimes(T,X);
        Y=Ttimes(T,temp')';
        temp=Ttimes(Ti,Y',1)';
        Xr=Ttimes(Ti,temp,1);
        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
            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 4 ****************************************************** 4 **

if (TestNo==4)   % methods included in Ttimes
    disp(['Test ',int2str(TestNo),' using methods included in Ttimes.']);
    if p1==0; p1=N; end;  % makes 16x16 DCT an often used default
    if p1==0
        Message=[Mfile,': an extra argument (p1) must be given, used as T in Ttimes.'];
        disp(Message);
        return
    end
    disp([Mfile,', Test ',int2str(TestNo),' : use method ',int2str(p1),' in Ttimes.']);
    % argument N is ignored for this case, N is calulated from p1
    N=2^(mod(p1-1,4)+1);
    if (Signal==1)
        Y=Ttimes(p1,X);
        if (p1==1)
            Xr=Ttimes(-p1,Y,0.5);
        else
            Xr=Ttimes(-p1,Y);
        end
        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 ((p1>1) & (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)
        temp=Ttimes(p1,X);
        Y=Ttimes(p1,temp')';
        if (p1==1)
            temp=Ttimes(-p1,Y',0.5)';
            Xr=Ttimes(-p1,temp,0.5);
        else
            temp=Ttimes(-p1,Y')';
            Xr=Ttimes(-p1,temp);
        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 ((p1>1) & (exist('geomean.m')==2))
            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 5 ****************************************************** 5 **

if (TestNo==5)   % methods included in Decom1D
    Method=N;
    disp([Mfile,', Test ',int2str(TestNo),' : use method ',int2str(Method),' in Decom1D.']);
    if (Signal==1)
        Y=Decom1D(Method,X,p1,p2);
        Xr=Decom1D(-Method,Y,p1,p2);
        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 ((Method<255) & (exist('geomean.m')==2))
            N=size(Y,1);
            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)
        Message=[Mfile,': 2D signal can not be used in test ',int2str(TestNo)];
        disp(Message);
    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 6 ****************************************************** 6 **

if (TestNo==6)   % methods included in Decom2D
    Method=N;
    if (Method>4)
        N=2^(mod(Method-1,4)+1);
    elseif Method==2
        N=2^p2;
    end
    disp([Mfile,', Test ',int2str(TestNo),' : use method ',int2str(Method),' in Decom2D. N=',int2str(N),'.']);
    if (Signal==1)
        Message=[Mfile,': 1D signal can not be used in test ',int2str(TestNo)];
        disp(Message);
    elseif (Signal==2)
        Y=Decom2D(Method,X,p1,p2);
        Xr=Decom2D(-Method,Y,p1,p2);
        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 ((Method>1) & (exist('geomean.m')==2))
            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 999 **************************************************** 999 **

if (TestNo==999)    % this test is not ready yet
    if (Signal==1)
        disp([Mfile,', Test ',int2str(TestNo),' : 1D signal is not ready.']);
    elseif (Signal==2)
        disp([Mfile,', Test ',int2str(TestNo),' : 2D signal is not ready.']);
    else
        Message=[Mfile,': Signal has illegal value in test ',int2str(TestNo)];
        disp(Message);
        return
    end
    Message=[Mfile,': Test ',int2str(TestNo),' finished ok.'];
end

disp(Message);
return