Matlab Code for Week 7 - Adam Oberman Math

%% Code for deniosing a signal using optimization
% The natural way to solve this problem is with simple filter.
% we are doing it this way because it generalizes to more complicated
% signals, e.g. signals with discontinuities, where smoothing is a problem.
%
% You may also run the MATLAB demo
% xpsound
 
%% Pure Note
Fs = 4000;
t = (1:Fs)/Fs;
yt = sin(Fs*t/20)';
%soundsc(yt,Fs);
n = Fs;
delta=40;
 
%% Add noise
A = .2;
noise = randn(size(yt));
yn = yt + A*noise;
%soundsc(yn,Fs);
NormNoise =  norm(noise);
 
%%  Now make the difference matrix M
e = ones(n,1);
Dx = spdiags([-e e], 0:1, n, n);
 
% Correct constraints on the endpoints.
Dx(n,:) = [];
 
%%  Now minimize the "Harshness" norm, 
cvx_begin
    variable x(n);
%    minimize ( sum_square(x - yn) + delta*sum_square(Dx*x) );
   minimize ( sum_square(x - yn) + delta*norm(Dx*x,1) );
cvx_end
 
%% plot just a portion of signal to make it easier to visualize
i = 400:600;
%figure(2), plot(t(i),yn(i),'r', t(i),x(i),'b'), axis tight
figure(3), plot(t(i),yt(i),'r', t(i),x(i),'b'), axis tight

%% Code for deniosing a signal using optimization
% The natural way to solve this problem is with simple filter.
% we are doing it this way because it generalizes to more complicated
% signals, e.g. signals with discontinuities, where smoothing is a problem.
%
% You may also run the MATLAB demo
% xpsound
%
 
%% Pure Note
Fs = 1024;
t = (1:Fs)/Fs;
yt = sin(5*Fs*t)';
soundsc(yt,Fs);
n = Fs;
 
 
%% Add noise
A = .2;
noise = randn(size(yt));
yn = yt + A*noise;
soundsc(yn,Fs);
NormNoise =  norm(noise);
 
%%  Now make the difference matrix M
e = ones(n,1);
Dx = spdiags([-e e], 0:1, n, n);
Dx = Dx/n;
% Correct constraints on the endpoints.
Dx(n,:) = [];
 
%%  Now minimize the "Harshness" norm, 
% subject to the "Amplitude" norm is close to the noisy signal.  
% How close? Based on the amplitude of the  noise.
cvx_begin
    variable x(n);
    minimize( sum_square (Dx*x) );
    subject to
    norm(x - yn) <= .25*NormNoise;
cvx_end
 
%% Note, depending on the value of the constraint parameter,
% the solution is smaller amplitude,
% but the shape of the solution seems to fit well with the original.
% we can simply rescale to fix the amplitude:
x = x/max(abs(x));
 
%% plot just a portion of signal to make it easier to visualize
i = 400:500;
figure(2), plot(t(i),yn(i), t(i),x(i)), axis tight
figure(3), plot(t(i),yt(i), t(i),x(i)), axis tight
soundsc(yn,Fs), soundsc(x,Fs);
soundsc(yt,Fs)

%% quadratic versus L1 smoothing Bi objective optimization.
 
% toggle this parameter to make the signal sharp
makesharp = 1;
% Choose 2 norm or 1 norm regularization
p = 1; %p = 2;
 
%% The smooth signal
n = 400;
t = linspace(0,1,n);
Sexact = 2*sin(2*pi*t).*sin(0.01*n*t);
Sexact = Sexact';
 
%% make signal sharp
if makesharp, 
    Sexact = round(Sexact);
end
%% The corrupted signal
Scorrupt =  Sexact + 0.05*randn(size(Sexact));
figure(1), plot(t, Sexact, t, Scorrupt);
 
%% The smoothing matrix
e = ones(n,1);
D = spdiags([-e e], -1:0, n, n);
 
%% run for several values of parameter or just one
if 1
    nopts = 6;
    lambdas = logspace(-1,1,nopts);   
else  % the best seems to be near lambda = 1
    lambdas = 1;    nopts = 1;
end
 
%% perform the minimization
for i=1:nopts
    disp(['* delta = ' num2str(lambdas(i))]);
    cvx_begin
    variable x(n)
    minimize ( norm(x - Scorrupt) + lambdas(i)*norm(D*x,p) )
    cvx_end
 
    % plot the solution against the signal
    figure(3), plot(t,Sexact,t,x);  
    title(['delta = ', num2str(lambdas(i)), ' p = ', num2str(p)]);
    pause(.5)    
end

%% quadratic versus L1 smoothing Bi objective optimization.
% Choose 2 norm or 1 norm regularization
p = 1; 
 
%% The smooth signal
n = 400;
t = linspace(0,1,n);
Sexact = 2*sin(2*pi*t).*sin(0.01*n*t);
Sexact = Sexact';
 
%% make signal jumpy and with a straight line
    Sexact = round(Sexact);
    t0 = .4;
    ii = t< t0;
    Sexact(ii) = (2/t0)*t(ii);
 
%% Add noise
Scorrupt =  Sexact + 0.05*randn(size(Sexact));
figure(1), plot(t, Sexact, t, Scorrupt);
 
%% The smoothing matrix, compare the two
e = ones(n,1);
Dx = spdiags([-e e], -1:0, n, n);
% Or make it
Dxx = spdiags([e -2*e e], -1:1, n,n);
 
%% Choose delta
delta1 = .5;
delta2 = .5;
 
%% perform the minimization
cvx_begin
variable x(n)
minimize ( norm(x - Scorrupt) + delta1*norm(Dx*x,p) )
cvx_end
 
x1 = x;
 
%% Compare to the other Smoothing matrix
cvx_begin
variable x(n)
minimize ( norm(x - Scorrupt) + delta2*norm(Dxx*x,p) )
cvx_end
 
x2 = x;
 
%% Combine both
cvx_begin
variable x(n)
minimize ( norm(x - Scorrupt) + .5*delta1*norm(Dx*x,p) + .5*delta2*norm(Dxx*x,p) )
cvx_end
x3= x;
 
 
%% plot the solution against the signal
figure(2), 
plot(t, x1, t, x2, t, x3); axis([0 1 -2 2]);