%% In this file we use various methods to reconstruct undersampled versions
% of the Scarlatti piano song, 	Sonata in D minor K. 9, as played by 
% A. Benedetti Michelangeli.
clear all; clc;
load pianoBasis.mat

addpath ./GPSR_6.0

% Load in the first few seconds of the Scarlatti song.
scarlatti = wavread('scarlattiOrig.wav');
n = length(scarlatti);

% We can reduce the basis size since the clip of the Scarlatti song is 
% not as long.
Uscar = pianoBasis(:,1:n);
T = length(t);
M = round(n/T);

% Let's listen to the best approximation of our Scarlatti song with the
% basis. This uses FULL signal information:

ww = Uscar'\scarlatti;
sound(Uscar'*ww,fs)
wavwrite(Uscar'*ww,fs,'scarlattiBest.wav');


%% Here we create a measurement matrix that simulates undersampling. Again
% there are two types of undersampling: random measurements and uniform
% measurements. To undersample, we take every 40th sample, which gives us 
% slightly more than 6000 samples total.

% Every iunder-th sample
iunder = 40;
ii=iunder:iunder:n;
R_meas = sparse(1:length(ii),ii,ones(size(ii)),length(ii),n);


%% For the l1 minimization code, we need to fit the data using a matrix
% product of our measurement matrix and our known basis.

R = R_meas*Uscar';

hR = @(x) R*x;
hRt = @(x) R'*x;

y = R_meas*scarlatti;

% l1 regularization parameter
tau = 0.1*max(abs(R'*y));
debias = 1;
stopCri = 4;

[x,x_debias,obj,...
    times,debias_start,mse]= ...
         GPSR_BB(y,hR,tau,...
         'Debias',debias,...
         'AT',hRt,... 
         'Monotone',1,...
         'Initialization',0,...
         'StopCriterion',stopCri,...
         'MaxiterA', 2000, ...
         'ToleranceD',0.001);

% l1 result
sound(Uscar'*x_debias,fs);
wavwrite(Uscar'*x_debias, fs, 32, 'scarlatti_6063unif_l1.wav');


%% l2 result
xl2 = R\y;

% Debias the l2 result
inds = find(xl2>0);
xl2(inds) = R(:,inds)\y;

% Let's hear how it sounds:
sound(Uscar'*xl2,fs);
wavwrite(Uscar'*xl2, fs, 32, 'scarlatti_6063unif_l2.wav')

%% Linear interpolation result

tfull = 1/fs:1/fs:(0.25*M);
tii = tfull(ii);
y_interp = interp1(tii,y,tfull,'linear');
sound(y_interp,fs)

wavwrite(y_interp, fs, 'scarlatti_6063unif_interp.wav')

