In [8]:
using PyPlot;
srand(1);
n = 100;
x = zeros(n);
sel = randperm(n);
x[sel[1:4]] = randn(4);
dftmtx = fft(eye(n,n),1);
A = inv(dftmtx);
y_true = A*x;
y = y_true + 0.004*randn(n);
figure();
plot(1:n, y);
title("Measured sound signal");
xlabel("time");
figure();
plot(1:n, x);
xlabel("frequency");
title("Original frequency components");
In [16]:
# fully sampled case
xhat = A\y;
figure();
plot(1:n, x, label="Original signal");
plot(1:n, xhat, label="L2-reconstruction (fully sampled)");
legend();
title("Reconstructed frequency components");
figure();
plot(1:n, y_true, label="Original signal");
plot(1:n, dftmtx*xhat/n, label="L2-reconstruction (fully sampled)");
legend();
title("Reconstructed sound signal");
In [17]:
# underdetermined case (randomly take 30% of samples), Tykhonov reg.
using Convex, ECOS
m2 = round(Int64, 0.3*n);
sel = randperm(n);
y2 = y[sel[1:m2]];
A2 = A[sel[1:m2],:];
xh = Variable(n);
#objective = norm(A2*xh-y2)^2 +1.4e-3*sum(abs(xh));
objective = norm(A2*xh-y2)^2 +1e-3*norm(xh);
problem = minimize(objective);
solve!(problem);
xhat2 = evaluate(xh);
figure();
plot(1:n, x, label="Original signal");
plot(1:n, xhat2, label="L2-reconstruction (30% undersampled)");
legend();
title("Reconstructed frequency components");
figure();
plot(1:n, y_true, label="Original signal");
plot(1:n, dftmtx*xhat2/n, label="L2-reconstruction (30% undersampled)");
legend();
title("Reconstructed sound signal");
In [18]:
# underdetermined case (randomly take 30% of samples), LASSO
using Convex, ECOS
xh = Variable(n);
objective = norm(A2*xh-y2)^2 +1e-3*sum(abs(xh));
problem = minimize(objective);
solve!(problem);
xhat2 = evaluate(xh);
figure();
plot(1:n, x, label="Original signal");
plot(1:n, xhat2, label="L1-reconstruction (30% undersampled)");
legend();
title("Reconstructed frequency components");
figure();
plot(1:n, y_true, label="Original signal");
plot(1:n, dftmtx*xhat/n, label="L1-reconstruction (30% undersampled)");
legend();
title("Reconstructed sound signal");
