%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "Insert title here" 
% Written for GNU Octave
% 
% INSTITUTE OF TELECOMMUNICATIONS  
% University of Stuttgart 
% www.inue-uni-stuttgart.de 
% author: Alexander Knaub  
% date: 05.10.2016
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This example script just prints the values of the given parameters.
% Make sure that the function prototype matches the prototype defined
% in the XML config file. (Parameter names are arbitrary but should be the same
% for better readability)

function [oversampling_used, pdf_link] = ...
    quantization_noise(...
        signal,... 
        clipping, ...
        oversampling,... 
        n_steps,... 
        behav,... 
        round_opt,... 
        method,... 
        tikz_ena,... 
        filename)

if is_octave()
    pkg load signal
end
% define formats
landscape = '-S930,350';
portrait  = '-S640,480';

% select format
output_format = portrait;

% Create an invisible figure.
fh = figure(1);
set(fh, "visible", "off");

% define style


oversampling_used = round(oversampling);
if oversampling_used < 1
    oversampling_used = 1;
elseif oversampling_used > 32
    oversampling_used = 32;
end;

%  n_steps = 8;
n_samples = 30000;
n_visible_samples = 200;
t = linspace(0,2,n_visible_samples)';
t_max = t(end);


if strcmp(signal, 'static')
    %t = linspace(0,2,n_samples);
    x = -sin(t*2*pi)-sin(t*7*pi-pi/6)*0.7;

elseif strcmp(signal, 'inout')
    x = linspace(-2,2, length(t))';

elseif strcmp(signal, 'uni')
    diststr = 'uniform';
    x = rand(n_samples,1)*2;
    if oversampling_used > 1
	x_orig = x;
        x = resample(x, oversampling_used, 1);
    end
    x = x - mean(x);
    x_orig = x_orig - mean(x);
    
    
elseif strcmp(signal, 'gauss') || strcmp(signal, 'compact') || strcmp(signal, 'ulaw')
    diststr = 'normal';
    % create a normal distributed signal with unit power:
    x = randn(n_samples,1);
    % eliminate the values above the clipping threshold:
    x_orig = x;
    if oversampling_used > 1
        x = resample(x, oversampling_used, 1);
    end
    
%  elseif 
%      x = randn(n_samples,1);
end

fprintf('Signal variance: %0.2f\n', var(x));

% do clipping or companding if necessary:
mu = 7;
if strcmp(signal, 'compact')
    x_clip = erf(x/clipping);
    x_range = 2;
elseif strcmp(signal, 'ulaw')
    
    x_clip = sign(x) .* (log(1+mu*abs(x)) ./ log(1+mu));
    x_range = 2*max(abs(x_clip));
%      x_clip(abs(x_clip) > 1) = sign(x(abs(x_clip)>1));
else
    x_clip = x;
    x_clip(abs(x) > clipping) = sign(x(abs(x)>clipping)) * clipping;
    x_range = 2*max(abs(x_clip));
%      x_clip(abs(x_clip) > clipping) = sign(x_clip(abs(x_clip)>clipping)) * clipping;
end;

% prepare the quantization parameters:
quantization_step = x_range/n_steps;

tic

%  if round_opt && mod(n_steps, 2) == 0
%      fprintf('Rescaling for better rounding...\n');
%      % renormalize the signal so that it lies within a smaller range:
%      scaling_factor = quantization_step * (n_steps/2-0.5)*(1-eps) / max(abs(x_clip))
%      x_clip = x_clip * scaling_factor;
%      x = x * scaling_factor;
%  end

if (~strcmp(method, 'round')  && ~strcmp(method, 'dn'))
    fprintf('Unknown method %s! Applying rounding.\n', method)
    method = 'round';
end

if strcmp(method, 'round')
    fprintf('Quantization by rounding\n')
    
    if strcmp(behav, 'riser')
        fprintf('Mid-riser behavior\n')
        
        if mod(n_steps, 2) == 0
            fprintf('Even # of steps.\n')
            x_quant = (floor(x_clip/quantization_step) + 0.5) * quantization_step;
            % make sure that we have the correct number of steps:
            x_quant(x_quant > quantization_step*(n_steps/2.0 - 0.5))... 
                = quantization_step * (n_steps/2.0 - 0.5);
        else
            fprintf('Odd # of steps.\n')
            quantization_step = x_range/(n_steps+1)
            
            if round_opt 
                scaling_factor = quantization_step * ((n_steps-1)/2) / max(abs(x_clip));
                x_clip = x_clip * scaling_factor;
                x = x * scaling_factor;
            end
            
            x_quant = floor((x_clip)/quantization_step)*quantization_step + 0.5*quantization_step;
            
            % make sure that we have the correct number of steps:
            x_quant(x_quant > quantization_step*((n_steps-1)/2 - 0.5))... 
                = quantization_step*((n_steps-1)/2 - 0.5);
        end
    else
        fprintf('Mid-tread behavior\n')
        
        if mod(n_steps, 2) == 0
            fprintf('Even # of steps.\n')
            quantization_step = x_range/(n_steps+1);
            if round_opt 
                fprintf('Rescaling for better rounding...\n');
                % renormalize the signal so that it lies within a smaller range:
                scaling_factor = quantization_step * (n_steps/2-0.5) / max(abs(x_clip));
                x_clip = x_clip * scaling_factor;
                x = x * scaling_factor;
%                  clipping = scaling_factor*clipping;
            end
            x_quant = floor((x_clip)/quantization_step+0.5)*quantization_step;
            % make sure, the range after quantizing is [-N/2 ... N/2-1]:
            x_quant(x_quant > quantization_step*(n_steps/2 - 1))... 
                = quantization_step * (n_steps/2 -1);
        else
            fprintf('Odd # of steps.\n')
            quantization_step = x_range/(n_steps);
            x_quant = floor((x_clip)/quantization_step+0.5)*quantization_step;
            % make sure, the range after quantizing is [-N/2 ... N/2-1]:
            x_quant(x_quant > quantization_step*(n_steps/2 - 0.5))... 
                = quantization_step * (n_steps/2 -0.5);
        end
    end
    
else
    fprintf('Quantization by truncating\n')
    
    if strcmp(behav, 'riser')
        fprintf('Mid-riser behavior\n')
        x_quant = floor((x_clip)/quantization_step)*quantization_step;
        if mod(n_steps, 2) == 0
            x_quant(x_quant > quantization_step*(n_steps/2 - 1))... 
                = quantization_step * (n_steps/2 - 1);
        else
            x_quant(x_quant < -quantization_step*((n_steps-1)/2))... 
                = -quantization_step * ((n_steps-1)/2);
        end
    else
        fprintf('Mid-tread behavior\n');
        x_quant = floor((x_clip)/quantization_step+0.5)*quantization_step - 0.5*quantization_step; 
        if mod(n_steps, 2) == 0
            x_quant(x_quant > quantization_step*(n_steps/2 - 0.5))... 
                = quantization_step * (n_steps/2 - 0.5);
            x_quant(x_quant < -quantization_step*(n_steps/2 - 0.5))... 
                = -quantization_step * (n_steps/2 - 0.5);
        else
            x_quant(x_quant > quantization_step*(n_steps/2 - 1))... 
                = quantization_step * (n_steps/2 - 1);
        end
    end
    
end

% do expanding:
if strcmp(signal, 'compact')
    x_quant = erfinv(x_quant)*clipping;
elseif strcmp(signal, 'ulaw')
    x_quant = sign(x_quant)*(1/mu) .* ((1+mu).^abs(x_quant/quantization_step) - 1);
end
if oversampling_used >1
quantization_error = resample(x_quant,1,oversampling_used)-x_orig;
else
quantization_error = x_quant-x;
end
snr = 10*log10(mean(x_orig.^2) / mean((quantization_error).^2));

%  fprintf('min/max error: %0.2f/%0.2f\n',... 
%          min(quantization_error)/quantization_step,...
%          max(quantization_error)/quantization_step)
%  
%  fprintf('quantized error: %0.2f/%0.2f\n',... 
%          min(quantization_error)/quantization_step,...
%          max(quantization_error)/quantization_step)

err_hist_bin_width = 0.05;

% if companding and expanding has been applied, quantization step is 
% not constant and normalization has to be modified
if strcmp(signal, 'compact')
    normalize_value = quantization_step*clipping;
else
    normalize_value = quantization_step;
end

first_err_tap = max(-4, min(-1,... 
                    (round((min(quantization_error)/normalize_value)/err_hist_bin_width)-2)*err_hist_bin_width));

last_err_tap = min(4, max(1,... 
                    (round((max(quantization_error)/normalize_value)/err_hist_bin_width)+2)*err_hist_bin_width));
                    
err_hist_taps = [first_err_tap:err_hist_bin_width:last_err_tap];

[err_nn, err_xx] = hist(quantization_error/normalize_value, err_hist_taps);

err_nn(length(err_nn)-1) = sum((quantization_error/normalize_value) > (last_err_tap-(3*err_hist_bin_width/2)));
err_nn(end) = 0;
err_nn = err_nn/sum(err_nn);


hist_taps = [min(x_quant)-2*normalize_value ...
            : normalize_value/4 ...
            : max(x_quant) + 2*normalize_value];

[x_nn, x_xx] = hist(x, hist_taps);
x_nn(length(x_nn)-1) = sum(x > (hist_taps(end-1)-normalize_value/8));
x_nn(length(x_nn)) = 0;
x_nn = x_nn/sum(x_nn);


[quant_nn, quant_xx] = hist(x_quant, hist_taps);
quant_nn = quant_nn/sum(quant_nn);

fprintf('Calculation time: %0.1fms\n', toc*1000);






if tikz_ena
    tic
    % read the template into the string variable 'tikz_code':
    basename = strrep(filename, '.png', '');
    tikz_template = fopen('tikz_template.tex');
    nlines = fskipl (tikz_template, Inf);
    frewind(tikz_template)
    tikz_code = '';
    for line=0:nlines
        cur_line = fgets(tikz_template);
        tikz_code = sprintf('%s%s', tikz_code, cur_line);
    end;
    fclose(tikz_template);
    
    % export first subplot: signals and error vs. time:
    time_plot_data = [t(t<(2))/oversampling_used, ...
                        x(t<(2)), ... 
                        x_quant(t<(2)), ... 
                        quantization_error(t<(2))];
    time_plot_file = [basename '_time_plot.dat'];
    save('-ascii', time_plot_file, 'time_plot_data');

    % display the clipping thresholds:
    tikz_code = strrep(tikz_code, '#clipping#',... 
                    sprintf('%0.2f', clipping));
    
    % replace some text placeholders with actual values:
    tikz_code = strrep(tikz_code, '#sigdist#', diststr);
    
    tikz_code = strrep(tikz_code, '#nsteps#',... 
                    sprintf('%d', n_steps));
    tikz_code = strrep(tikz_code, '#nbits#',... 
                    sprintf('%0.1f', log2(n_steps)));
    tikz_code = strrep(tikz_code, '#stepsize#',... 
                    sprintf('%0.1f', quantization_step));
    
    % axis scaling
    tikz_code = strrep(tikz_code, '#tmax#',... 
                    sprintf('%0.2f', t_max/oversampling_used));

    % export second subplot: histogram data of original and quantized signal:
    signal_hist_data = [hist_taps'/normalize_value, x_nn', quant_nn'];
    signal_hist_file = [basename '_signal_hist.dat'];
    save('-ascii', signal_hist_file, 'signal_hist_data');
    
    tikz_code = strrep(tikz_code, '#sighistmax#',... 
                    sprintf('%0.2f', max(quant_nn)*1.4));

    % export third subplot: histogram data of the quantization error:
    error_hist_data = [err_hist_taps', err_nn'];
    error_hist_file = [basename '_error_hist.dat'];
    save('-ascii', error_hist_file, 'error_hist_data');
    
    tikz_code = strrep(tikz_code, '#errhistmax#',... 
                    sprintf('%0.2f', max(err_nn)*1.3));
    tikz_code = strrep(tikz_code, '#errbound#',... 
                    sprintf('%0.2f', max(abs(err_hist_taps))*1.3));
    
%      tikz_code = strrep(tikz_code, '#errhist1stbin#',... 
%                      sprintf('%0.2f', first_err_tap));
%      tikz_code = strrep(tikz_code, '#errhistlastbin#',... 
%                      sprintf('%0.2f', last_err_tap));

    % make TIKZ read correct files:
    tikz_code = strrep(tikz_code, '#time_data_file#', time_plot_file);
    tikz_code = strrep(tikz_code, '#err_hist_file#', error_hist_file);
    tikz_code = strrep(tikz_code, '#sig_hist_file#', signal_hist_file);

    
    % set some variables:

    
    tikz_code = strrep(tikz_code, '#snr#', sprintf('%0.2f', snr));
    tikz_code = strrep(tikz_code, '#mse#', sprintf('%0.2f', mean(quantization_error.^2)));
    
    

    tikz_file = fopen(sprintf('%s.tex', basename), 'w');
    fwrite(tikz_file, tikz_code);
    fclose(tikz_file);

    cmd = 'pdflatex -output-directory=./tmp %s.tex; pdf2ps %s.pdf %s.ps; ps2eps -q %s.ps; convert -density 200 -resize 640x480 %s.eps %s.tikz.png ';
    cmd = strrep(cmd, '%s', basename); 

    fprintf('about to execute the command\n%s', cmd);
    [status, output] = system(cmd);
    copyfile(sprintf('%s.tikz.png', basename), filename);
    
    pdf_link = ['<a href=./',sprintf('%s.pdf', basename),'> Hi-Res PDF image</a>'];
    
    toc
else
    
    pdf_link = '';
    % Plot into the figure
    tic
    set(0, 'defaultlinelinewidth', 2);
    set(0, 'defaultaxesfontsize', 6);
    set(0, 'defaultaxesfontname', 'Helvetica');
    set(0, 'defaulttextfontsize', 8);
    set(0, 'defaulttextfontname', 'Helvetica');
    %set(fh, 'papersize', [640/150, 480/150])
    set(fh, 'paperposition', [0,0, 6.4, 4.8])
    subplot(2,2,[1 2], 'align');
    plot(t(t<1.5), x(t<1.5)); hold all;
    plot(t(t<1.5), x_quant(t<1.5));
    plot(t(t<1.5), quantization_error(t<1.5));
    grid on;
    xlim([0,1.5])
    xlabel('Time')
    l = legend('Original', 'Quantized', 'Quant. error', 'location', 'northeastoutside')
    %  set(l, 'LineSpacing', 0)
    title('Time domain signal')

    subplot(2,2,3, 'align');
    bar(x_xx, quant_nn, 'r');
    hold all
    bar(x_xx, x_nn, 'b');
    legend('Quantized', 'Original', 'location', 'northeastoutside')
    title('Signal distribution')
   

    subplot(2,2,4, 'align');

    %bar(err_xx, err_nn, 'b')
    plot(20*log10(abs(fftshift(fft(quantization_error)))));
    hold all;
    plot(20*log10(abs(fftshift(fft(x)))));
    %  x_pos = min(err_xx) + 0.25*(max(err_xx) - min(err_xx));
    legend(sprintf('MMSE: %0.2f', mean(quantization_error.^2)), 'location', 'northeastoutside')
    %  plot(linspace(-0.5, 0.5, length(t))*1/(t(2)-t(1)), fftshift(abs(fft(quantization_error))))
    % xlim([-0.5, 0.5])
    %  xlabel('Normalized frequency')
    title(sprintf('SNR=%0.2fdB', snr))
    toc


    % Create an image file. This png image is displayed on the webpage. Its size 
    % must be 640x480.

    print(filename, '-dpng', '-r100');
    %  print(filename, '-dpng');
    %  print([filename '.pdf'], '-dpdf', '-color');
end

end

%% subfunction that checks if we are in octave
function r = is_octave()
  persistent x;
  if (isempty (x))
    x = exist ('OCTAVE_VERSION', 'builtin');
  end
  r = x;
end
