%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% .m script of web demo "Sampling of analog signals" 
% Written for GNU Octave
% 
% INSTITUTE OF TELECOMMUNICATIONS  
% University of Stuttgart 
% www.inue-uni-stuttgart.de 
% author: Alexander Knaub
% date: 10.10.2014

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% BEGIN, main function (control)
function [] = sampling_1(...
    signal_type,... % type of the signal. possible: 1) 4-PAM
    f_signal,... % maximum signal frequency which should be in the generated signal
    fs,... % sampling frequency
    fc,... % Cut-off frequency of the reconstruction low pass filter.
    alpha_lpf,... % Quality of reconstruction filter. Better for smaller values.
    display_options,... % switch of subplots to be displayed
    show_reconstr_spect_en,... % enable ploting the spectrum of reconstructed signal
    show_lpf_tf_en,... % enable ploting the transfer function of reconstruction filter
    webdemo_en,... % if 0, then test-mode in Octave/Matlab commandline, no webserver running
    svg_en,... % if 1, then also output svg-file
    filename... % output file for webdemo
)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % BEGIN check input parameters for consistency
    td_en = bitand(display_options, 1);
    fd_en = bitand(display_options, 2);
    % END check input parameters for consistency
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %% Create signals
    if signal_type == 1
        [signal, n_analog, fs_analog, n_min_disp, n_max_disp, h] = generate_4pam(f_signal);
    end;
    
    
    %% Perform sampling in time domain:
    sampling_factor = fs_analog/fs % NON-INTEGER: Take care when using it!
    % sample time instances needed to place points in the plot correctly
    % in digital domain, there is no relation to the continuous time
    n_sampled = zeros(1, floor(length(n_analog)/sampling_factor));
    for i=1:length(n_sampled)
        n_sampled(i) = n_analog(floor(i*sampling_factor));
    end;
    
    signal_sampled = zeros(1, floor(length(n_analog)/sampling_factor));
    for i=1:length(signal_sampled)
        signal_sampled(i) = signal(floor(i*sampling_factor));
    end;
    
    % we describe sampled signal by inserting zeros to see the aliases
    signal_upsampled = zeros(1, length(signal));
    for i= 1:sampling_factor:length(signal)
        % make signal nonzero on sampling points:
        signal_upsampled(floor(i)) = signal(floor(i));
    end;
    
    f_analog = linspace(-fs_analog/2, fs_analog/2, length(signal));
    analog_spectrum = fftshift(fft(signal, length(signal)));

    sampled_spectrum = fftshift((fft(signal_upsampled, length(signal))))*sampling_factor;
    
    
    % create filter for signal reconstruction. 
    % Use RC filter with bandwidth = symb_rate=2*fc:
    [h_lpf] = get_rc_imp_resp(alpha_lpf, 20, fs_analog, 2*fc);
    % transfer function of the filter in frequency domain:
    tf_h_lpf = fftshift(abs(fft(h_lpf, length(f_analog))));
    tf_h_lpf = tf_h_lpf/max(tf_h_lpf);
    
    % do actual filtering:
    signal_reconstructed = conv(signal_upsampled, h_lpf, 'same')*fc*2/fs;
    reconstructed_spectrum = fftshift(fft(signal_reconstructed, length(signal)));
    
    
    
    %% Display the results
    % Create a figure.
    if (webdemo_en)
        fig = figure(1); set(fig, 'visible', 'off');
    else
        fig = figure(1); set(fig, 'visible', 'on');
    end;
    
    set(0, 'defaultlinelinewidth', 2);


    if (webdemo_en)
        set(0, 'defaultaxesfontsize', 8);
        set(0, 'defaultaxesfontname', 'Helvetica');
        set(0, 'defaulttextfontsize', 8);
        set(0, 'defaulttextfontname', 'Helvetica');
        const_lw = 2;
    elseif (svg_en)
        set(0, 'defaultaxesfontsize', 6);
        set(0, 'defaultaxesfontname', 'Helvetica');
        set(0, 'defaulttextfontsize', 6);
        set(0, 'defaulttextfontname', 'Helvetica');
        const_lw = 1;
    else
        set(0, 'defaultaxesfontsize', 10);
        set(0, 'defaultaxesfontname', 'Helvetica');
        set(0, 'defaulttextfontsize', 8);
        set(0, 'defaulttextfontname', 'Helvetica');
        const_lw = 2;
    end;
    
    if td_en & fd_en
        subplot(2, 1, 1)
    end
    
    if td_en
        plot(n_analog, signal(1:length(n_analog)), 'linewidth', 3);
        title(['Time domain sampling and corresponding spectrum, (f_s/f_{max} = ', num2str(fs/f_signal), ')']);
        xlabel('time t \cdot f_{max}');
        ylabel('x(t)');
        axis_ylimit = max(max(abs(signal_reconstructed)), max(abs(signal)));
        axis([n_min_disp n_max_disp -axis_ylimit axis_ylimit]);
        set(gca,'xtick',[-2:1:10])
        hold all; grid on;
        stem(n_sampled, signal_sampled(1:length(n_sampled)), 'color', 'm', 'markersize', 3);
        plot(n_analog, signal_reconstructed(1:length(n_analog)), 'r')
        legend('original (analog)', 'sampled', 'reconstructed');
        hold off;
    end
    if td_en & fd_en
        subplot(2, 1, 2)
    end
    if fd_en
        n_aliases = ceil(5/fs) % the smaller the frequency, the more aliases are visible
        RGBpalette=jet(n_aliases);
        norm_value = max(abs(analog_spectrum));
        analog_db = 20*log10(abs(analog_spectrum)/norm_value);
        h1 = fill(f_analog, analog_db, 'b');
        axis([-5, 5, -15, 2]);
        set(h1, 'EdgeColor','none');
        set(h1, 'FaceAlpha',0);
        set(gca,'xtick',[-5:1:5])
        xlabel('frequency f/f_{max}');
        ylabel('|X_1(f)|^2 [dB]');
        
        hold all; grid on;
        
        for curr_alias = 1:n_aliases
            RGBcolor = RGBpalette(n_aliases-curr_alias+1,:);
            h1 = fill(f_analog+curr_alias*fs, analog_db, RGBcolor);
            set(h1, 'FaceColor', RGBcolor);
            set(h1, 'linewidth', 0.1);
            set(h1, 'EdgeColor','none');
            set(h1, 'facealpha',0.7);
            % add one point to complete the polygon
            h1 = fill([-5, f_analog-curr_alias*fs], [-15, analog_db], RGBcolor);
            set(h1, 'FaceColor', RGBcolor);
            set(h1, 'linewidth', 0.1);
            set(h1, 'EdgeColor','none');
            set(h1, 'facealpha',0.7);
        end;
        h1 = fill(f_analog, analog_db, 'b');
        set(h1, 'EdgeColor','none');
        set(h1, 'FaceAlpha',1);
        
        if show_reconstr_spect_en == 1
            plot(f_analog, 20*log10(abs(reconstructed_spectrum)/norm_value), 'r','linewidth', 4);
        end;
        if show_lpf_tf_en == 1
            plot(f_analog, 20*log10(tf_h_lpf), 'c--','linewidth', 4);
        end;
    %      plot(f_analog, 20*log10(abs(sampled_spectrum)), 'linewidth', 1, 'r');
        hold off;
    end;
    
    print(filename, '-dpng','-S640,480', '-F:8');
    
    if (svg_en)
        print([filename, '.svg'], '-dsvg','-S640,480');
    end;
    
end
% END, main program (control)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, function definitions
%% Helper function to generate the underlying symbols
function [signal, n_signal, fs_analog, n_min_disp, n_max_disp, h] = generate_4pam(f_max)
    
    % in baseband, bandwidth is only half the symbol rate
    symb_rate = 2*f_max; % actualy *1/(1+alpha), here alpha nearly 0
    
    % we always want to have same number of (analog) samples, independent of 
    % symbol rate => oversampling factor is assigned dynamically:
    % higher symbo rate => less time (samples) between symbols.
    ovsamp = 52;
    
    
    % allow RC filter with narrow roll-off:
    alpha = .0001;
    % define RC length 
    len = 6; % has to be even to avoid jumps
    normalize_en = 0;
    complex_en = 0;
    pam_4_sequence = [1, 3, -1, 1, 1, -1, -3, -3, -1, -3, -1, 3, 3, -3, 1, 1,... 
                      -3, -3, -3, 1, -1, 1, -3, -3, -1, -1, 1, 3];
    pam_4_upsampled = upsample(pam_4_sequence, ovsamp/2);
    % remove trailing zeros
    pam_4_upsampled = pam_4_upsampled(1:end-(ovsamp/2-1));
    
    % create the time axis for the signal:
    i_seq_min = -4;
    i_seq_max = i_seq_min + length(pam_4_sequence) - 1;
    n_seq = [i_seq_min:i_seq_max] / symb_rate;
    n_signal = linspace(n_seq(1), n_seq(end), length(pam_4_upsampled));
    
    fs_analog = 1 / (n_signal(2) - n_signal(1));
    
    [h] = get_rc_imp_resp(alpha, len, fs_analog, symb_rate);
    
    
    % 'same' to remove ramp up
    signal = conv(pam_4_upsampled, h, 'same');
    
    n_min_disp = 0;
    n_max_disp = min(n_signal(end), 10);
    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Pulse shaping filter that is used to obtain a band limited signal
% note: len is nb of taps without considering oversampling
function [h] = get_rc_imp_resp(alpha, len, fs_h, symb_rate)
    % create the time axis for the filter:
    t_min = -len / (2*symb_rate);
    t_max = len / (2*symb_rate);
    
    n_h = [t_min:1/fs_h:t_max];
    
    h = (sinc(n_h * symb_rate) .* cos(pi*alpha*n_h*symb_rate)) ./ (1 - (2*alpha*n_h*symb_rate).^2);
    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% insert zeroes
function y = upsample(x, usf)
  len = length(x);
  y = zeros(1,len*usf);
  for k=0:len-1
    y(1+k*usf) = x(1+k);
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% END, function definitions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%