clear all
close all
clc

x_orig = zeros(1024,6);
load ex4mwden
clear x
whos
for e = 1:30
    x_orig(:, 5) = load_signal('Piece-Regular', 1024);
    x_orig(:, 6) = load_signal('cusp', 1024);
    kp = 0;
    kk = 6;
    amp_perc = 0.1;
    for m = 1:kk
        sigma(m) = amp_perc*max(x_orig(:, m));
    end
    x(:, 1) = x_orig(:, 1) + sigma(1)*randn(size(x_orig(:, 1)));
    x(:, 2) = x_orig(:, 2) + sigma(2)*randn(size(x_orig(:, 2)));
    x(:, 3) = x_orig(:, 3) + sigma(3)*randn(size(x_orig(:, 3)));
    x(:, 4) = x_orig(:, 4) + sigma(4)*randn(size(x_orig(:, 4)));
    x(:, 5) = x_orig(:, 5) + sigma(5)*randn(size(x_orig(:, 5)));
    x(:, 6) = x_orig(:, 6) + sigma(6)*randn(size(x_orig(:, 6)));
    
    
    if amp_perc == 0.1
        iter = [1 3 1 1 1 150];
    elseif amp_perc == 0.2
        iter = [5 30 10 10 10 300];
    elseif amp_perc == 0.3
        iter = [10 25 10 15 15 300];
    end
    for k = 1:kk
        
        gamma = 0.5;
        g = firpm(6,[0 0 .45 .5]*2,[0 0 1 1]);
        r = round(size(g, 2)/2);
        im1 = x(:,k);
        
        for j = 1:iter(k)
            im1 = padarray(im1, r-1, 'symmetric');
            x_hpf = conv(im1, g, 'same'); %high-passed signal
            x_or = im1; %main signal
            x_hpf_col = [x_hpf; 0];
            
            for i = 1:2
                x_sgn = [sign(x_hpf_col(1:end-1)); -1];
                x_sgn(x_sgn==0) = 1;
                xxp = [x_hpf_col(1:end-1); 0];
                x_in = (xxp'*x_sgn)/numel(xxp);
                x_p = xxp - x_in.*x_sgn;
                x_pp = x_p.*x_sgn;
                x_hpf_col = [x_p(1:end-1); sum(x_pp(1:end-1))];
                
                
                x_sgn_new = sign(x_hpf_col);
                arg = x_sgn_new-x_sgn;
                [ind1 val1] = find(arg == 2);
                [ind2 val2] = find(arg == -2);
                x_hpf_col(ind1) = 0;
                x_hpf_col(ind2) = 0;
                
            end
            sol = perform_l1ball_projection([x_p(1:end-1);0], sum(x_pp(1:end-1)));
            x_lpf = x_or - x_hpf; %low-passed signal
            x_lpf = [x_lpf; 0];
            xl = sol;
            xl = xl + x_lpf;
            
            im1 = xl(r:end-r);
        end
        imm(:,k) = im1;
    end
    
    
    for d = 1:kk
        SNR_in(d) = snr(x_orig(:,d), x(:,d));
        SNR_out(d) = snr(x_orig(:,d), imm(:,d));
    end
    
    SNR_in1(e, :) = SNR_in;
    SNR_out1(e, :) = SNR_out;
end

% Plotting the original, noisy, and denoised signals
kp = 0;
figure(2)
for f = 1:kk
    subplot(kk,3,kp+1), plot(x_orig(:,f)); axis tight;
    title(['Original signal ',num2str(f)])
    subplot(kk,3,kp+2), plot(x(:,f)); axis tight;
    title(['Observed signal ',num2str(f)])
    subplot(kk,3,kp+3), plot(imm(:,f)); axis tight;
    title(['Denoised signal ',num2str(f)])
    kp = kp + 3;
end

% Calculating the averaged SNR
SNR_In_Ave = mean(SNR_in1, 1)
SNR_Out_Ave = mean(SNR_out1, 1)

% The average SNR for all 6 signals
mean(SNR_Out_Ave)