% Napredne metode digitalne obradbe signala - Projekt
% Compressed sensing - Potpuna rekonstrukcija uz smanjeni broj uzoraka
% Davor Bonaci, davor.bonaci@fer.hr

%% 1. Initialization
clear all;
close all hidden;
pack;

fs = 1e4;   % Sampling frequency, in Hz
N = 32;     % Number of samples
M = 10;     % Number of measurements
f1 = 100;   % Sine frequency one, in Hz
f2 = 400;   % Sine frequency two, in Hz

%% 2. Input signal and its spectrum
t = [0 : 1/fs : (N-1)/fs];
x = sin(2*pi*f1*t) + cos(2*pi*f2*t);
figure, plot(t, x), grid, title('Real "analog" signal'), 
xlabel('Time, seconds'), ylabel('Amplitude');

f = fs * (0 : N/2) / N;
X = fft(x);
figure,
subplot(2, 1, 1), plot(f, abs(X(1 : N/2+1))), grid, 
title('Amplitude spectrum of the real "analog" signal'), 
xlabel('Frequency, Hz'), ylabel('Amplitude'),
subplot(2, 1, 2), plot(f, rad2deg(angle(X(1 : N/2+1)))), grid,
title('Phase spectrum of the real "analog" signal'), 
xlabel('Frequency, Hz'), ylabel('Phase, deegres');

%% 3. Computing DFT matrix and representation of the input signal
%% in the psi-domain
%psi_inv = dftmtx(N);
%psi = conj(psi_inv) / N;
psi_inv = dctmtx(N)';
psi = dctmtx(N);

s = (psi_inv * x.').';

%% 4. Computing measurement matrix and performing measurements
randn('state', 0);
fi = 1/N .* randn(M, N);
y = fi * x.';

%% 5. Reconstructing initial signal
theta = fi * psi;

cvx_begin
    variable s_est(N, 1);
    minimize( norm(s_est, 1) );
    subject to
        theta * s_est == y;
cvx_end

x_est = (psi * s_est).';

%% 6. Verification
max_error_s = max(abs(s - s_est.'))
max_error_x = max(abs(x - x_est))
max_imag_x = max(imag(x_est))

figure, hold on, grid, title('Input signal and reconstruction in the psi-domain'),
plot(s, 'b'), plot(s_est, 'g'), plot(s_est.' - s, 'r'),
legend('Input', 'Reconstruction', 'Error');

figure, grid, hold on, title('Input signal and reconstructed signal'), 
xlabel('Time, seconds'), ylabel('Amplitude'),
plot(t, x, 'b'), plot(t, x_est, 'g'), plot(t, x-x_est, 'r'),
legend('Input', 'Reconstruction', 'Error');
