clc
clear
close all



%% % LOADING AND ANALZYING THE SIGNAL, SHOW SPECTRUM % % %
 song = 'humenecko' % 'abba' % 'speech' % 'guitar' % 
 
switch song %setting the input - sparse audio signal
    case 'abba'
        [y, Fs] = audioread('test_1.wav');

    case 'humenecko'
        [y, Fs] = audioread('test_2.wav');
        
    case 'speech'
        [y, Fs] = audioread('test_3.wav');
        
    case 'guitar'
        [y, Fs] = audioread('test_4.wav');
        
end
% sound(y, Fs);

N = length(y); %length of the signal
k = 5;  %sparsity of the signal in the dictionary

y_show = ones (20, 1);
show_length = Fs * 2;
starting_point_orig = N / 2;
ending_point_orig = starting_point_orig + show_length - 1;
y_show(1:show_length) = y(starting_point_orig:ending_point_orig); %at least two beats - average slow song has 80bpm - we need at least two seconds of the song. There are Fs samples per second - we need to display double the value
fig_spect = figure;
Y = sgram(y_show, Fs, 'dynrange',50);
hold on;
% title('spectrum of the part of the input signal');
title('Gaborovo spektrum části signálu');

bits = '10' % '4' % '8' % '16' % '24'% %'32'% %setting the bit depth of quantization processed after application of fade-out effect
% % % \LOADING AND ANALZYING THE SIGNAL, SHOW SPECTRUM % % %

%% % DEFINITION OF THE ENVELOPE % % %

b = floor(N/5); %beginning of the envelope
l = N-b; %lenght of the envelope
cenv = 1/(l+1); %coefficient of the envelope
env = ones (N, 1); %prepares the vector

for j = b:N
    env(j) = 1-((j-b)*cenv);  %vymyslet efektivnější řešení
end

ienv = ones (N, 1);
ienv(b:N) = 1./env(b:N);
Check_env_ienv = env .* ienv;

% % % \DEFINITION OF THE ENVELOPE % % %


%% % APPLICATION OF THE ENVELOPE, SHOW RESULTS % % %

y_env = y .* env;
fig_orig_env = figure;
A = plot(y);
hold on;
B = plot(y_env);
% title('original vs. affected signal');
title('Původní signál vs. efektovaný signál');

% % % \APPLICATION OF THE ENVELOPE, SHOW RESULTS % % %


%% % SETTING THE SIGNAL BOUNDARIES, QUANTIZE FADED SIGNAL % % %

switch bits
   
    case '4'
        d=2^4;
        
    case '8'
        d=2^8;
        
    case '10'
        d=2^10;
    
    case '16'
        d=2^16;
        
    case '24'
        d=2^24;
        
    case '32'
        d=2^32;
end

        min_y_env = min(min(y_env));    
        max_y_env = max(max(y_env));
        range_env = max_y_env - min_y_env;
        
        quant_step_env = range_env/(d-1);
        dec_bounds_env = (min_y_env + quant_step_env/2) : quant_step_env : (max_y_env - quant_step_env/2);  %decision boundaries lie in the middle
        quant_levels_env = min_y_env : quant_step_env : max_y_env; %quantization levels
        
        lower_dec_bounds_env = y_env - (quant_step_env/2);
        upper_dec_bounds_env = y_env + (quant_step_env/2);
        
        [index_env, y_quant_env] = quantiz(y_env, dec_bounds_env, quant_levels_env);
        y_quant_env = y_quant_env';
       
        min_quant_level_env = quant_levels_env(1);
        max_quant_level_env = quant_levels_env(end);
        upper_dec_bounds_env(upper_dec_bounds_env > max_quant_level_env) = max_quant_level_env;
        lower_dec_bounds_env(lower_dec_bounds_env < min_quant_level_env) = min_quant_level_env;

        % fig_orig_quant = figure;
        % plot(y_env);
        % hold on;
        % plot(y_quant_env);
        % title('affected and quantized signal');        


fig_analog_digital = figure;
plot(y_env-y_quant_env);
% title('difference after artificial quantization');
title('Rozdíly po umělé kvantizaci');

% % %  \SETTING THE SIGNAL BOUNDARIES, QUANTIZE FADED SIGNAL % % %

%% % Multiplication only% % %

y_multiplied = y_quant_env .* ienv;

% % % \Multiplication only % % %
%% % DEFINITON OF THE FRAME FOR THE SPARSE DEQUANTIZATION % % %
% [param.w, param.a] = win_param(64, fs, 4);
length_win = 30;
rel_shift = 4;
        
w = length_win*Fs/1000;
w_temp = (w/rel_shift - floor(w/rel_shift)); 
    if w_temp ~= 0
        w = (round(w/rel_shift))*rel_shift;
    end

a = w/rel_shift;
M = w;
% h_support =  1 + 3*rand([500 1]);
% 
% test = fft(h_support);

F = frame('dgtreal', {'hann',w}, a, M);
F = frametight(F);

% % % /DEFINITON OF THE FRAME FOR THE SPARSE DEQUANTIZATION % % %



%% % SPARSE DEQUATISATION, FADE-OUT RECONSTRUCTION % % %
        param.lower_lim = lower_dec_bounds_env;
        param.upper_lim = upper_dec_bounds_env;
        
        param.lower_lim = param.lower_lim .* ienv;
        param.upper_lim = param.upper_lim .* ienv;
        
        indi_thr = 10e-5; % threshold for identifying the point to lie in the set
        
        f1.eval = @(x) norm(x,1);
        f1.prox = @(x,T) sign(x).*max(abs(x)-T, 0);

        f2.eval = @(x) eps;
        f2.prox = @(x,T) x - frana(F, ((postpad(frsyn(F, x),N) - proj_box(postpad(frsyn(F, x),N),[],param)))); 

        
        %%%%%%%% UNLOCBOX version %%%%%%%%%%%%%
        % setting different parameter for the simulation
        paramsolver.verbose = 5;  % display parameter
        paramsolver.maxit = 300;        % maximum iteration
        paramsolver.tol = 10e-7;        % tolerance to stop iterating
        paramsolver.lambda = 1;   % step for DR algorithm (default 1)
        paramsolver.gamma = 1e-2;        % here threshold for soft thresholding
        
        [sol, info] = douglas_rachford(frana(F, y_quant_env), f1, f2, paramsolver);
        info
        
        sol = f2.prox(sol,[]); %final projection into the constraints
        y_dequant = frsyn(F, sol);
        y_dequant = postpad(y_dequant,N);

% % % \SPARSE DEQUATISATION, FADE-OUT RECONSTRUCTION % % %

fprintf('hotovo')

%% % ANALYZING RESULTS - SNR % % %
    
    NormY = norm(y);
    NormYxYdeq = norm(y-y_dequant);
    NormYxYmult = norm(y-y_multiplied);

    snr_deq = 20 * log10((norm(y))/(norm(y-y_dequant)));
    snr_multi = 20 * log10((norm(y))/(norm(y-y_multiplied)));
    
% % %     snr_dequa=snr(y, y-y_dequant);
% % %     snr_multip=snr(y, y-y_multiplied);
    
% % % \ANALYZING RESULTS - SNR % % %
%% % SHOW RESULTS % % %

%define colors
grey = 0.6* ones(1,3);
lightgrey = 0.8* ones(1,3);
black = [0,0,0];
blue = [0.251,0.584,0.808];
orange = [0.851,0.325,0.098];
green = [0 1 0];

fig_y_Previous = figure;
plot(y);
title('Signál původní - před obálkou');
        
fig_y_Gab = figure;
plot(y_dequant);
title('Signál rekonstruovaný pomocí řídké dekvantizace');

fig_time = figure;
h_orig = plot(y, '.-', 'Color', black);
hold on;
h_orig_env = plot(y_env, '.-', 'Color', blue);
h_quant_env = plot(y_quant_env, '.-', 'Color', green);
h_dequant = plot(y_dequant, '.-', 'Color', orange);
% title('Original, faded, quantized and reconstructed signal');
title('Původní, efektovaný, kvantizovaný a rekonstruovaný signál');

uistack(h_orig_env,'top');
uistack(h_dequant,'top');
% h_legend_time_env = legend([h_orig, h_orig_env, h_quant_env, h_dequant ], 'original', 'original - fade out', 'quantized - fade out', 'reconstructed');
h_legend_time_env = legend([h_orig, h_orig_env, h_quant_env, h_dequant ], 'Původní', 'Původní efektovaný', 'Kvantovaný efektovaný', 'Rekonstruovaný');

axis tight

fig_diff = figure;
h_difference_multiplied = plot(y-y_multiplied, '.-', 'Color', orange);
hold on;
h_difference_sparse_deq = plot(y-y_dequant, '.-', 'Color', blue);
% title('Difference between original and reconstructed signal');
title('Rozdíl mezi původním a rekonstruovaným signálem');
% h_legend_diff = legend([h_difference_sparse_deq, h_difference_multiplied], 'sparse dequantization', 'multiplied only');
h_legend_diff = legend([h_difference_sparse_deq, h_difference_multiplied], 'Řídká dekvantizace', 'Pouhé pronásobená');

% fig_diff_2 = figure;
% h_difference = plot(y_multiplied-y_dequant, '.-', 'Color', blue);
% hold on;
% % title('Difference between multiplied and sparsly dequantized');
% title('Rozdíl mezi pronásobeným a řídce dekvantovaným');
% h_legend_diff = legend([h_difference], 'diff multi deq');

fig_orig_deq = figure;
h_orig = plot(y, '.-', 'Color', black);
hold on;
h_dequant = plot(y_dequant, '.-', 'Color', orange);
% title('Original and sparsly dequantized signal');
title('Původní a řídce dekvantovaný signál');
h_legend_orig_deq = legend([h_orig, h_dequant], 'Původní', 'Řídce dekvantovaný');

fig_orig_multi = figure;
h_orig = plot(y, '.-', 'Color', black);
hold on;
h_multi = plot(y_multiplied, '.-', 'Color', blue);
% title('Original and multiplied signal');
title('Původní a pronásobený signál');
h_legend_orig_multi = legend([h_orig, h_multi], 'Původní', 'Pronásobený');


% % % \SHOW RESULTS % % %