function [signal_in, signal_noised, signal_out, koef_in, koef_out,  psnr, psnr_sr] = signal_denoising_impl(delka_signalu, ridkost, baze, wavelet, input_snr, hloubka_decompozice)

disp(['Delka signalu: ' num2str(delka_signalu)]);
disp(['Ridkost: ' num2str(ridkost)]);
disp(strcat('Baze: ', baze));
disp(strcat('Wavelet: ', wavelet));
disp(['Hloubka dekompozice: ' num2str(hloubka_decompozice)]);
disp(['Inout_snr: ' num2str(input_snr)]);



x=0:1:delka_signalu;

N = delka_signalu;
k = ridkost;
J = hloubka_decompozice;

J = 3;
% krty krty krt...
switch baze
    case 'DCT'
    U = dctmtx(N)';
    
    case 'DWT'
    [Lo_D,Hi_D,Lo_R,Hi_R] = w_filters(wavelet); 
    U = dwtmatrix(Lo_D,Hi_D,J,N)';
    
end
    
%% Original_signal

x_signal = zeros(N, 1);
perm = randperm(N); %random positions of nonzeros
%vector of coordinates contains random values from [-1,1] at nonzero positions
indexes = perm(1:k);
x_signal(indexes) = 2*rand(k,1)-1;
% max_freq = max(indexes) %maximum 'frequency'

signal_in = U * x_signal;

%% Coeficient_orig
%koeficienty ridkeho signalu
switch baze
    case 'DCT'
         koef_in = abs(dct(signal_in));
    
    case 'DWT'
         koef_in = x_signal;
    
end


%% Add a noise

% Add Gaussian i.i.d. noise

sigma_noise = 10^(-input_snr/20)*std((signal_in));
noise = (randn(size(signal_in)))*sigma_noise/sqrt(2);
signal_noised = signal_in + noise;


% Tolerance on noise
epsilon = 0.4*sqrt(chi2inv(0.99, 2*numel(signal_in))/2)*sigma_noise;

%% Encoding (measurements)
y = signal_noised;


%% Decoding 

% Solve BPDN problem (without weights) by Douglas-Rachford splitting algorithm

        A_ = @(x_) U*x_; 
        At_ = @(x_) U'*x_;

    % Sparsity operator
        Psit_ = @(x_) x_; 
        Psi_ = Psit_;

    % Proximity operator
        prox_lambda=0.0001;
        Psit_ = @(x_) prox_l1(x_, prox_lambda);
        Psi_ = Psit_;

    % Parameters for BPDN
        param.verbose = 1;      %  0 no log, 1 print main steps, 2 print all steps
        param.tol_l1 = 1e-8;    % is stop criterion for the loop. The algorithm stops if
    %
    %                            (  n(t) - n(t-1) )  / n(t) < tol,
    %      
    %                          where  n(t) = Psi(x)|| is the objective function at iteration t*
    %                          by default, tol=10e-4.
        param.maxit_l1 = 300;    %  max. nb. of iterations (default: 200)
        param.nu_l1 = 1;
        param.tight_l1 = 1;



sparse_coord = solve_bpdn(y, epsilon, A_, At_, Psi_, Psit_, param);

signal_out = U*sparse_coord;               %reconstructed sparse vector of coordinates



% sparse_coord_lin = solver_l1_linprog(signal_noised,U,epsilon);          %l1-minimization by linear programming
% signal_out_lin_prog = U*sparse_coord_lin;                      %reconstructed sparse vector of coordinates

%% Coeficient_recontr.
%koeficienty zrekonstruovaneho ridkeho signalu
switch baze
    case 'DCT'
         koef_out = abs(dct(signal_out));
    
    case 'DWT'
         koef_out = sparse_coord;
    
end


%% PSNR rekonstrukce ku originalu
signal_in    = double(signal_in);
signal_out = double(signal_out);
absD = abs(signal_in-signal_out);
moc    = absD.^2;
mse  = sum(moc(:))/numel(signal_in);
psnr = 10*log10((max(signal_in)-min(signal_in))^2/mse);

%% PSNR_sr      % puvodni signalu k zasumenemu signalu
signal_in    = double(signal_in);
signal_noised = double(signal_noised);
absD = abs(signal_in-signal_noised);
moc    = absD.^2;
mse  = sum(moc(:))/numel(signal_in);
psnr_sr = 10*log10((max(signal_in)-min(signal_in))^2/mse);






