% Demo ukazuje komprimacni vzorkovani na jednoduche simulaci, v ruznych
% moznych variantach. Baze pro signal je DCT.

clc
clear all
close all

%% Parameters
N = 256         %signal length
k = 8            %sparsity
m = 256         %number of measurements
% sampling = 'random'
sampling = 'equidistant'
grafy_matic = false;

%% Baze
U = dctmtx(N)'; %basis for the signal
if grafy_matic
    figure
    imagesc(U)
    axis image
    colorbar
    title('Baze')
end

%% Original signal
x = zeros(N, 1);
perm = randperm(N); %random positions of nonzeros
%vector of coordinates contains random values from [-1,1] at nonzero positions
indexes = perm(1:k);
x(indexes) = 2*rand(k,1)-1;
max_freq = max(indexes) %maximum 'frequency'

signal = U * x; %forming signal

figure
subplot(2, 1, 1);
plot(signal);
axis tight
hold('on');
title('Originalni signal v casove oblasti');
subplot(2, 1, 2);
stem(abs(dct(signal)));
axis tight
hold('on');
title('Originalni signal v kmitoctove oblasti (modul kosinove transformace)');

%% Measurement matrix
A = zeros(m, N);
% a/ random sampling:
if strcmp(sampling,'random')
    samples = randperm(N);
    samples = sort(samples(1:m));
end
% b/ equispaced sampling:
if strcmp(sampling,'equidistant')
    samples = 1:(N/m):N;
    samples = round(samples(1:m));
end

for ii = 1:m
    A(ii, samples(ii)) = 1;
end

if grafy_matic
    figure
    imagesc(A)
    axis image
    colorbar
    title('Merici matice')
end

%% Encoding (measurements)
y = A * signal;

%% Decoding (signal recovery from measurements)
f = ones(2*N, 1);
lb = zeros(2*N, 1);
v = linprog(f, [], [], [A*U -A*U], y, lb); %l1-minimization by linear programming
signal_ = U*[eye(N) -eye(N)]*v; %reconstructed sparse vector of coordinates

%% Results in graphs
figure
subplot(2, 1, 1);
plot(signal);
axis([0 N-1 min([signal; signal_]) max([signal; signal_])]);
hold('on');
stem(samples, signal(samples), 'LineStyle', 'none', 'MarkerEdgeColor', 'r', 'MarkerFaceColor', 'r', 'MarkerSize', 5);
hold('off');
title('Originalni signal');
subplot(2, 1, 2);
plot(signal_);
axis([0 N-1 min([signal; signal_]) max([signal; signal_])]);
hold('on');
stem(samples, signal(samples), 'LineStyle', 'none', 'MarkerEdgeColor', 'r', 'MarkerFaceColor', 'r', 'MarkerSize', 5);
hold('off');
title(['Rekonstrukce pomoci l1, l2-chyba: ' num2str(100*round(norm(signal-signal_)*100)/100)]);
% title(['Rekonstrukce s l2-chybou ' num2str(100*round(norm(signal-signal_)/norm(signal)*100)/100) ' %']);