%% 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')