function[signal_in, signal_noised, signal_out,  psnr, psnr_sr] = signal_denoising_2D_impl(image, ridkost, rozmer,original_picture,baze, wavelet, hloubka_decom, input_snr)

global M;
global N;
disp(strcat('image: ', image));
disp(['Ridkost: ' num2str(ridkost)]);
disp(['Rozmer: ' num2str(rozmer)]);

disp(strcat('Real: ', original_picture));
disp(strcat('Baze: ', baze));
disp(strcat('Wavelet: ', wavelet));
disp(['Hloubka Dekompozice: ' num2str(hloubka_decom)]);

disp(['Input_snr: ' num2str(input_snr)]);


J = hloubka_decom;

L = rozmer;
k = ridkost;

%% devise of picture
velikost_deleni = 16;
switch  L < velikost_deleni 
    case true
        velikost_deleni = L;
        m_ = 1;
        n_ = 1;
    case false
        m_ = M/velikost_deleni;
        n_ = N/velikost_deleni;
      
end

%% Base
 
switch baze
    case 'DCT'
    U = dctmtx(velikost_deleni)';                 % baze pro zpracovani 
    U_or = dctmtx(L)';                            % baze pro vygenerovani ridkeho obrazu
     
    case 'DWT'
    [Lo_D,Hi_D,Lo_R,Hi_R] = w_filters(wavelet); 
     U = dwtmatrix(Lo_D,Hi_D,J,velikost_deleni)';               % baze pro zpracovani 
     U_or = dwtmatrix(Lo_D,Hi_D,J,L)';                          % baze pro vygenerovani ridkeho obrazu
   
    
end


%% Original_signal
switch image
    case 'nahodny'
        N = L;
        M = L;
        S = N*N;
        
        x_signal = zeros(S, 1);
        perm = randperm(S); %random positions of nonzeros
        indexes = perm(1:k);
        x_signal(indexes) = randi(255,k,1);    %2*rand(k1,1)-1; %vector of coordinates contains random values from [-1,1] at nonzero positions
        
        koef = reshape(x_signal,M,N);

        original_picture_n = U_or' * koef * U_or; 
                 
        signal_in = original_picture_n(:);
    
    case 'real'
    
            switch size(size(imread(original_picture)),2)
    
                case 3
     
                original_picture_n = im2double(rgb2gray(imread(original_picture)));

                case 2
                 original_picture_n = im2double(imread(original_picture));
            end

            original_picture_n = imresize(original_picture_n,[64,64]);

            M = size(original_picture_n,1);   % size of picture M x N 
            N = size(original_picture_n,2);
 
            signal_in = original_picture_n(:);
end

%% Dictionary
 A = kron(U',U');

 %% 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 = signal_in + noise;
signal_noised = reshape(signal,M,N);

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



%% Denoising

signal_out = zeros(M,N);
for i = 1:m_
    for j = 1:n_
        
        m_dilci = ((i-1)*velikost_deleni + 1):i*velikost_deleni;
        n_dilci = ((j-1)*velikost_deleni + 1):j*velikost_deleni;
        
        part_of_picture = signal_noised(m_dilci,n_dilci);
        y = part_of_picture(:);
        
        % Solve BPDN problem (without weights) by Douglas-Rachford splitting algorithm

        A_ = @(x_) A*x_; 
        At_ = @(x_) A'*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 = 200;    %  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_v = A*sparse_coord;                      %reconstructed sparse vector of coordinates        

signal_out(m_dilci,n_dilci) = reshape(signal_out_v,velikost_deleni,velikost_deleni);


    end
end

disp('konec');
%% PSNR
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 = double(signal);
absD = abs(signal_in-signal);
moc    = absD.^2;
mse  = sum(moc(:))/numel(signal_in);
psnr_sr = 10*log10((max(signal_in)-min(signal_in))^2/mse);




