% Number of possible frequencies
n = 50;
% Number of samples
m = 15;
% Number of non-zero frequencies
k = 3;

% The fixed frequencies
frequencies = (1:n);

% The sample points
samples = sort(pi*rand(1,m))';

A = sin( frequencies .* samples );

% Setup the true solution
x_plus = zeros(n,1);
I = randperm( n, k );
x_plus(I) = randn(k, 1);

% and the right-hand side
b = A*x_plus;

% Plotting
N = 201;
t = linspace( 0, pi, N );
figure(1); clf;
plot( t, sum(sin( t .* frequencies') .* x_plus ), '-' );
pause
hold on;
plot( samples, sum(sin( samples' .* frequencies') .* x_plus ), 'x' );
pause

% We use 100 steps of the Douglas-Rachford splitting method:
z = randn(n,1);
iters = 100;
gamma = 1;
I = eye(n);
O = zeros(m);

for i = 1:iters
	x = wthresh( z, 's', gamma );

	yl = [I A'; A O] \ [2*x - z; b];
	y = yl(1:n);

	z = z + y - x;
end

% Plot the solution
plot( t, sum(sin( t .* frequencies') .* x), '--' );
