function [signal_in, signal_tex, signal_car, signal_out, psnr, psnr_t, psnr_c, koef1, koef2, koef3, koef4] = ...
    decomposition_1D_impl(delka_signalu, ridkost_k1, ridkost_k2, baze_1, baze_2,wavelet_1, wavelet_2, hloubka_decom_1, hloubka_decom_2)

%global N;
global signal_in_c;
global signal_in_t;

disp(['Delka signalu: ' num2str(delka_signalu{1})]);

disp(['Ridkost k1: ' num2str(ridkost_k1)]);
disp(['Ridkost k2: ' num2str(ridkost_k2)]);


disp(strcat('Baze_1: ', baze_1));
disp(strcat('Wavelet_1: ', wavelet_1));
disp(['Hloubka Dekompozice_1: ' num2str(hloubka_decom_1)]);
disp(strcat('Baze_2: ', baze_2));
disp(strcat('Wavelet_2: ', wavelet_2));
disp(['Hloubka Dekompozice_2: ' num2str(hloubka_decom_2)]);


J_1 = hloubka_decom_1;
J_2 = hloubka_decom_2;

N = str2double(delka_signalu{1});
k1 = ridkost_k1;
k2 = ridkost_k2;

%% Base
 
switch baze_1
    case 'DCT'
    U_1 = dctmtx(N)';                 % baze pro zpracovani 
   
     
    case 'DWT'
    [Lo_D,Hi_D,Lo_R,Hi_R] = w_filters(wavelet_1); 
     U_1 = dwtmatrix(Lo_D,Hi_D,J_1,N)';               % baze pro zpracovani 

end

switch baze_2
    case 'DCT'
    U_2 = dctmtx(N)';         % baze pro zpracovani 
        
    case 'DWT'
    [Lo_D,Hi_D,Lo_R,Hi_R] = w_filters(wavelet_2); 
     U_2 = dwtmatrix(Lo_D,Hi_D,J_2,N)';               % baze pro zpracovani 
     
    
end

%% Original_signal

%1.komponenta
x_signal_1 = 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:k1);
x_signal_1(indexes) = 2*rand(k1,1)-1;

signal_in_t = U_1 * x_signal_1;

%2.komponenta
x_signal_2 = 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:k2);
x_signal_2(indexes) = 2*rand(k2,1)-1;

signal_in_c = U_2 * x_signal_2;

% vstupni signal vytvoren jako kombinace dvou ridkych komponent
signal_in = signal_in_t + signal_in_c;

%% Coeficient_orig
%koeficienty prvni komponenty
switch baze_1
    case 'DCT'
         koef1 = abs(dct(signal_in_t));
    
    case 'DWT'
         koef1 = x_signal_1;
    
end
%koeficienty druhe komponenty
switch baze_2
    case 'DCT'
         koef2 = abs(dct(signal_in_c));
    
    case 'DWT'
         koef2 = x_signal_2;
    
end


%% Decomposition

y = signal_in ;

 %% Initialisation  - UnLockBox
         % Author: Nathanael Perraudin
            % Date: 2 October 2013

%             % Loading toolbox
%             init_unlocbox();

            verbose = 1;    % verbosity level

        %% General parameter for the problem
            mu1=1;
            mu2=1;

        %% Creation of the Problem
           

            % Creation of the operator Psi
            Psi = @(x) U_1*x ;            %1/sqrt(N)*fft(x);
            Psit = @(x) U_1'*x  ;          %sqrt(N)*ifft(x);
            nu_Psi = 1; % Frame bound

            % Creation of th operator Phi
            Phi = @(x) U_2*x ;              %M.*Psi(x);
            Phit = @(x) U_2'*x ;              %Psit(M.*x);
            nu_Phi = 1; % Frame bound


           % Creation of the measurements
            s = y;

        %% Define the UNLocBoX proximity operators and gradient

            %Define mask functions
            Mc = @(x) x(1:N);         % select the variable c           
            Mb = @(x) x(N+1:end);     % select the variable b             
            Cc = @(x,c) [c;Mb(x)];    % recreate all opt variable from c               
            Cb = @(x,b) [Mc(x);b];    % recreate all opt variable from b                
            
        
            % Functional mu1*||c||_1
            param_l1_1.verbose = verbose - 1;
            f1.prox = @(x,T) Cc(x,prox_l1(Mc(x),T*mu1,param_l1_1));
            f1.eval = @(x)  mu1*norm(Mc(x),1);

            % Functional mu2*||c||_1
            param_l1_2.verbose = verbose - 1;
            f2.prox = @(x,T) Cb(x,prox_l1(Mb(x),T*mu2,param_l1_2));
            f2.eval = @(x)  mu2*norm(Mb(x),1);

            % Define the gradient of ||s - Psi c - Phi b||_2^2 
            g.grad = @(x)  [ Psit( Psi(Mc(x)) + Phi(Mb(x)) - s ) ; ...
                 Phit( Psi(Mc(x)) + Phi(Mb(x)) - s ) ];
            g.eval = @(x)    norm((s-(Psi(Mc(x))+Phi(Mb(x))))./sqrt(2))^2;   %/sqrt(2))


        %% Define solver parameter

            paramfb.maxit=1000;
            paramfb.gamma=1/(nu_Phi+nu_Psi);
            paramfb.verbose=verbose;

        % Solve the problem
            x0= eps*ones(2*N,1);
            sol=generalized_forward_backward(x0,{f1,f2},g,paramfb);
          
%         % predkoklad velke ridkost, hodne male hodnoty zanedbame         
%        for i = 1:2*N    
%             if sol(i) <= max(abs(sol))/4 && sol(i) >= -max(abs(sol))/3
%                 sol(i)= 0;
%             else
%                 sol(i)= 2*sol(i);
%             end
%        end
            uv = reshape(sol,N,2);
            
            sol_1 = uv(:, 1);
            sol_2 = uv(:, 2);
   
      
signal_tex = U_1 * sol_1;
signal_car = U_2 * sol_2;



signal_out = signal_tex + signal_car;

%% Coeficient_decom
%koeficienty prvni komponenty
switch baze_1
    case 'DCT'
         koef3 = abs(dct(signal_tex));
    
    case 'DWT'
         koef3 = sol_1;
    
end
%koeficienty druhe komponenty
switch baze_2
    case 'DCT'
         koef4 = abs(dct(signal_car));
    
    case 'DWT'
         koef4 = sol_2;
    
end

 %% PSNR
signal_out = double(signal_out);
signal_in    = double(signal_in);
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_T
%max_out = max(abs(signal_kvant));
signal_in_t    = double(signal_in_t);
signal_tex = double(signal_tex);
absD = abs(signal_in_t-signal_tex);
moc    = absD.^2;
mse  = sum(moc(:))/numel(signal_in_t);
psnr_t = 10*log10((max(signal_in_t)-min(signal_in_t))^2/mse);

%% PSNR_C
%max_out = max(abs(signal_kvant));
signal_in_c    = double(signal_in_c);
signal_car = double(signal_car);
absD = abs(signal_in_c-signal_car);
moc    = absD.^2;
mse  = sum(moc(:))/numel(signal_in_c);
psnr_c = 10*log10((max(signal_in_c)-min(signal_in_c))^2/mse);


