demo5.m

% % % % % % % % % % % % % % % % % % % % % % % % % % %
% Test sparsity closing to m/n=1
% % % % % % % % % % % % % % % % % % % % % % % % % % %

clear
clc

% % % % % % % % % % % % % % % % % % % % % % % % % % %
% Prepare raw data
% % % % % % % % % % % % % % % % % % % % % % % % % % %

RawInpLoad = load('15814m_ltdbECG_1h.mat');
RawInpLoad = RawInpLoad.val;
n_dl = 102;
epochs = floor(length(RawInpLoad) / n_dl);    % 4517
RawInpLoad = RawInpLoad(1:n_dl * epochs);

% % % % % % % % % % % % % % % % % % % % % % % % % % %
% Setting parameters for training
% % % % % % % % % % % % % % % % % % % % % % % % % % %

param.lambda = 0.15;            % sparsity constraint 
param.numThreads = -1; 
param.batchsize = 400;
param.verbose = false;
param.iter = 10; 

% % % % % % % % % % % % % % % % % % % % % % % % % % %
% Prepare training and testing data
% % % % % % % % % % % % % % % % % % % % % % % % % % %

RawInp = RawInpLoad(1:n_dl*epochs);
RawInp = reshape(RawInp , n_dl, epochs);
crossValidFactor = 0.7;
TrainInp = RawInp(:, 1:floor(epochs*crossValidFactor));
TrainInp = TrainInp - repmat(mean(TrainInp),[size(TrainInp,1),1]);
TrainInp = TrainInp ./ repmat(sqrt(sum(TrainInp.^2)),[size(TrainInp,1),1]);

TestInp = RawInp(:, (size(TrainInp,2)+1):epochs);
TestInp = TestInp - repmat(mean(TestInp),[size(TestInp,1),1]);
TestInp = TestInp ./ repmat(sqrt(sum(TestInp.^2)),[size(TestInp,1),1]);

% % % % % % % % % % % % % % % % % % % % % % % % % % %
% Compressive sensing
% % % % % % % % % % % % % % % % % % % % % % % % % % %
rsnr_dl = zeros(1,50);
res_dl = zeros(1,50);
sparsity_dl = zeros(1,50);
m_dl_cnt = 1;
rsnr = 0;
res = 0;
spar = 0;

param.K = 512;
disp('Starting to  train the dictionary');
D = mexTrainDL(TrainInp,param);
% alpha = mexLasso(TrainInp,D,param);
% MSE = mean(0.5*sum((TrainInp-D*alpha).^2));
% fprintf('objective function: %f\n',MSE);

psi_dl = D;

% for m_dl = floor(0.95*n_dl: 0.01*n_dl: n_dl)
for m_dl = floor(n_dl/50: n_dl/50: n_dl)
    phi_dl = randn(m_dl,n_dl);
    A_dl = phi_dl * psi_dl;

    for ep = 1:size(TestInp,2)
        y_dl = phi_dl * TestInp(:,ep);
        x0_dl = pinv(A_dl) * y_dl; 
        xs_dl = l1eq_pd(x0_dl, A_dl, [], y_dl); 
        xhat_dl = psi_dl * xs_dl;
        rsnr = rsnr + 20 * (log10 (norm(TestInp(:,ep),2) / norm(TestInp(:,ep) - xhat_dl,2)));   
        res = res + norm(TestInp(:,ep) - xhat_dl,2);
        spar = spar + length(find(abs(xs_dl)>0.001) );
    end
    rsnr_dl(m_dl_cnt) = rsnr / size(TestInp,2);
    res_dl(m_dl_cnt) = res / size(TestInp,2);
    sparsity_dl(m_dl_cnt) = 1 - spar / size(TestInp,2) / length(xs_dl);
    m_dl_cnt = m_dl_cnt + 1;
    rsnr = 0;
    res = 0;
    spar = 0;
 end


% % % % % % % % % % % % % % % % % % % % % % % % % % %
% Plot results
% % % % % % % % % % % % % % % % % % % % % % % % % % %

figure
plot(floor(n_dl/50: n_dl/50: n_dl)./n_dl,sparsity_dl);
xlabel('m/n, n = 102');
ylabel('L0-norm of coefficients');