% AC for toff_track; pulse_shaping(128, 1, 0, 0, 1, 0, 4, 8, 0.1, 0, 1, 1, 1, 1, 1, 2, 20, 0, 0, 0, 1, 0, 0,1, 0,[0 0 0 0 0 0 1], 0, 0, "pulse_shaping")
% pulse_shaping(100, 1, 0, 0, 1, 0, 8, 8, 0.1, 0, 1, 1, 1, 1, 1, 2, 20, 0, 0, 0, 1, 0, 0,1, 0,[0 0 0 0 0 0 1], 0, 0, "pulse_shaping")
% -- also allow oversamp as parameter

% AC-tests: pulse_shaping(10000, 1, 0, 0, 1, 0, 4, 8, 0.1, 0, 1, 1, 1, 1, 1, 2, 20, 0, 0, 0, 1, 0, 0,1, 0,[1 0 0 0 0 0 0], 0, 0, "pulse_shaping")
% SFO tests: pulse_shaping(10000, 1, 0, 0, 1, 0, 2, 8, 0.1, 0, 1, 1, 1, 1, 1, 2, 20, 0, 0, 0, 1, 0, 0,1, 0,[0 0 0 0 0 1 0], 0, 0, "pulse_shaping")

% to do
%- pulse shaping rect; 3dB,30dB BW
%- AWGN

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXAMPLE: pulse_shaping(16, 1, 0, 0, 1, 8, 8, 0.1, 0, 0, 1, 0, 1, 1, 2, [1 0 0 0 0 0 0], 0, 0, "pulse_shaping")
% PARAMETERS:
%   Nsym,                         % nb of PAM symbols
%   nb_bits_per_real_comp,        % nb of bits per PAM symbol
%   impulse_shape,                % R.C., triangular or rect. impulse
%   unipolar_en,                  % unpolar impulses also possible
%   Tsym,                         % symbol duration in seconds (set to 1 for normalization)
%   offset_QAM_en,                % shift Q-channel by T/2 if enabled    (for phase plot)
%   ovsamp,                       % nb of samples in between PAM symbols (oversampling)
%   h_rc_sym_len,                 % length of t-domain raised cosine filter in PAM-symbols
%   h_rc_alpha,                   % raised cosine filter roll-off
%   complex_en,                   % if 1 (enabled), generate complex QAM symbols
%   rnd_data_en,                  % enable random data
%   show_indiv_imp_en,            % show individual PAM-impulses in t-plot
%   exclude_rampup_en,            % if 1, then don't include filter rampup in t-domain signal
%   H_magn_sq_en,                 % plot |H(f)| (to better see raised cosine) or |H(f)|^2
%   spectrum_dB_en,               % plot spectrum y-axis in dB
%   spectrum_display_bw_factor,   % display "spectrum_display_bw_factor"-times the Nyquist bandwidth in spectrum
%   SNR_dB,                       % for scatter plot
%   phi_off,                      % for scatter plot
%   foff_over_Rs,                 % for scatter plot, =1/100000;
%   phi_imb,                      % for scatter plot, =pi/180*0; %(-10);   % Q ampl. imbalance
%   a_imb,                        % for scatter plot, =1; %1.2;       % Q phase imbalance
%   dc_imb,                       % for scatter plot, =0; %0.1;
%   t_samp_off,                   % for scatter plot, =0;
%  ToutOverTin,                   % SFO, e.g., 10001.0/10000 means, that output symbol is 1 sample longer every 10000 input samples; default is 1
%   colormap_idx,                 % for scatter plot
%   plot_en,                      % enable plots 0, 1 [t-domain, f-domain, ampl/phase diagram, h_rc_n, eye-diagram, scatter-plot, auto-corr]
%   svg_en,                       % if 1, then also output svg-file
%   webdemo_en,                   % if 0, then test-mode in Octave/Matlab commandline, no webserver running
%   filename                      % output file for webdemo

% examples
% t-domain:    pulse_shaping(16, 1, 0, 0, 1, 0, 8, 8, 0.1, 0, 0, 1, 0, 1, 1, 2, 200, 0, 0, 0, 1, 0, 0, 1, 0,[1 0 0 0 0 0 0], 0, 0, "pulse_shaping")
% f-domain:    pulse_shaping(200, 1, 0, 0, 1, 0, 8, 8, 0.5, 1, 0, 1, 0, 1, 1, 2, 200, 0, 0, 0, 1, 0, 0, 1, 0,[0 1 0 0 0 0 0], 0, 0, "pulse_shaping")
% Offset-QPSK: pulse_shaping(200, 1, 0, 0, 1, 0, 16, 8, 0.5, 1, 0, 1, 1, 1, 1, 2, 200, 0, 0, 0, 1, 0, 0,1, 0,[0 0 1 0 0 0 0], 0, 0, "pulse_shaping")
% single impulse: pulse_shaping(100, 1, 0, 0, 1, 0, 8, 8, 0.1, 0, 0, 1, 0, 1, 1, 2, 200, 0, 0, 0, 1, 0, 0,1, 0,[0 0 0 1 0 0 0], 0, 0, "pulse_shaping")
% eye-diagram:    pulse_shaping(100, 1, 0, 0, 1, 0, 8, 8, 0.1, 0, 1, 1, 1, 1, 1, 2, 200, 0, 0, 0, 1, 0, 0,1, 0, [0 0 0 0 1 0 0], 0, 0, "pulse_shaping")
% scatter-plot:   pulse_shaping(10000, 2, 0, 0, 1, 0, 32, 8, 0.1, 0, 1, 1, 1, 1, 1, 2, 20, 0, 0, 0, 1, 0, 0,1, 0,[0 0 0 0 0 1 0], 0, 0, "pulse_shaping")
% further extensions, or separate demos:
% - spectrum, with Hamming or Hanning window to avoid edge effects (option)
% - eye diagram (with/without channel; to illustrate 2nd Ny criterion); for fun: after AWGN channel with Eb/N0, Es/N0
% - also: plot constellation after simple AWGN channel; compute EVM in dB

% BEGIN, main function (control)
function pulse_shaping(
  Nsym,                         % nb of PAM symbols
  nb_bits_per_real_comp,        % nb of bits per PAM symbol
  impulse_shape,                 % R.C., triangular or rect. impulse
  unipolar_en,                  % unipolar impulses also possible
  Tsym,                         % symbol duration in seconds (set to 1 for normalization)
  offset_QAM_en,                % shift Q-channel by T/2 if enabled (for phase plot)
  ovsamp,                       % nb of samples in between PAM symbols (oversampling)
  h_rc_sym_len,                 % length of t-domain raised cosine filter in PAM-symbols
  h_rc_alpha,                   % raised cosine filter roll-off
  complex_en,                   % if 1 (enabled), generate complex QAM symbols
  rnd_data_en,                  % enable random data
  show_indiv_imp_en,            % show individual PAM-impulses in t-plot
  exclude_rampup_en,            % if 1, then don't include filter rampup in t-domain signal
  H_magn_sq_en,                 % plot |H(f)| (to better see raised cosine) or |H(f)|^2
  spectrum_dB_en,               % plot spectrum y-axis in dB
  spectrum_display_bw_factor,   % display "spectrum_display_bw_factor"-times the Nyquist bandwidth in spectrum
  SNR_dB,                       % for scatter plot
  phi_off_deg,                  % for scatter plot
  foff_over_Rs,                 % for scatter plot, =1/100000;
  phi_imb_deg,                  % for scatter plot, =pi/180*0; %(-10);   % Q ampl. imbalance
  a_imb,                        % for scatter plot, =1; %1.2;       % Q phase imbalance
  dc_imb,                       % for scatter plot, =0; %0.1;
  t_samp_off,                   % for scatter plot, =0;
  ToutOverTin,                  % SFO, e.g., 10001.0/10000 means, that output symbol is 1 sample longer every 10000 input samples; default is 1
  colormap_idx,                 % for scatter plot
  plot_en,                      % enable plots 0, 1 [t-domain, f-domain, ampl/phase diagram, h_rc_n, eye-diagram, scatter-plot]
  svg_en,                       % if 1, then also output svg-file
  webdemo_en,                   % if 0, then test-mode in Octave/Matlab commandline, no webserver running
  filename                      % output file for webdemo
)



%new
%SNR_dB=200;
phi_off=phi_off_deg/180*pi;
%foff_over_Rs=1/100000;
phi_imb=phi_imb_deg/180*pi; %(-10);   % Q ampl. imbalance
%a_imb=1; %1.2;       % Q phase imbalance
%dc_imb=0; %0.1;
%t_samp_off=0;

%colormap_idx=0;
colorbar_en=1; % fixed

normalize_PAM_en=1;

if (nb_bits_per_real_comp<0)
  normalize_PAM_en=0;
  nb_bits_per_real_comp=abs(nb_bits_per_real_comp);
end


if ((t_samp_off>0) && plot_en(6)) ovsamp=32 end   % to illustrate sampling frequ. offs.
if ((t_samp_off>0) && plot_en(7)) ovsamp=32 end   % to illustrate sampling frequ. offs.


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN check input parameters for consistency

  % only one plot at a time
  if (sum(plot_en)!=1)
    plot_en = [1 0 0 0 0 0 0];
    warning("setting  'plot_en = [1 0 0 0 0 0 0]'");
  end

  if (plot_en(1) && (Nsym>128))  Nsym = 128; warning("setting  'Nsym = 128'");  end

  if (Nsym>100000)  Nsym=100000; warning("setting  'Nsym = 100000'");  end
  if (Nsym<1)       Nsym=1;      warning("setting  'Nsym = 1'");       end

  if (nb_bits_per_real_comp>8)     nb_bits_per_real_comp=8;    warning("setting  'nb_bits_per_real_comp = 8'");     end
  if (nb_bits_per_real_comp<1)     nb_bits_per_real_comp=1;    warning("setting  'nb_bits_per_real_comp = 1'");     end

  if (ovsamp>32)     ovsamp=32;    warning("setting  'ovsamp = 32'");     end
  if (ovsamp<1)      ovsamp=1;     warning("setting  'ovsamp = 1'");      end

  if (h_rc_sym_len>64)  h_rc_sym_len=64;    warning("setting  'h_rc_sym_len = 64'");     end
  if (h_rc_sym_len<4)   h_rc_sym_len=4;     warning("setting  'h_rc_sym_len = 4'");      end

  if (h_rc_alpha>1)  h_rc_alpha=1;    warning("setting  'h_rc_alpha = 1'");     end
  if (h_rc_alpha<0)  h_rc_alpha=0;    warning("setting  'h_rc_alpha = 0'");     end


  if (spectrum_display_bw_factor<1)  spectrum_display_bw_factor=1; end
  if (spectrum_display_bw_factor>64)  spectrum_display_bw_factor=64; end


  % adjust input parameters
  if (h_rc_alpha==0) h_rc_alpha = 0.00001; end
  if (h_rc_alpha==1) h_rc_alpha = 0.9999;  end


  if (spectrum_dB_en)  H_magn_sq_en=1;  end   % OVERRIDE if in dB, can restrict to magn. sq.

  if (plot_en(5)) exclude_rampup_en=1; end     % eye-diagram; no ramp-up

  if (plot_en(6)) exclude_rampup_en=1; complex_en=1; end % for constellation plot, force to be complex

  if (plot_en(7)) exclude_rampup_en=1; complex_en=1; end % for constellation plot, force to be complex


  % END check input parameters for consistency
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % RGB color definitions, http://www.rapidtables.com/web/color/RGB_Color.htm;   use e.g.  as  magenta(7,:)
  grayc =   [[0,0,0]; [32,32,32]; [64,64,64]; [96,96,96]; [128,128,128]; [160,160,160]; [192,192,192]; [224,224,224]; [255,255,255]]/255;
  violett = [[51,0,25]; [102,0,51]; [153,0,76]; [204,0,102]; [255,0,127]; [255,51,153]; [255,102,178]; [255,153,204]; [255,204,229]]/255;
  magenta = [[51,0,51]; [102,0,102]; [153,0,153]; [204,0,204]; [255,0,255]; [255,51,255]; [255,102,255]; [255,153,255]; [255,204,255]]/255;
  lilac   = [[25,0,51]; [51,0,102]; [76,0,153]; [102,0,204]; [127,0,255]; [153,51,255]; [178,102,255]; [204,153,255]; [229,204,255]]/255;
  blue    = [[0,0,51]; [0,0,102]; [0,0,153]; [0,0,204]; [0,0,255]; [51,51,255]; [102,102,255]; [153,153,255]; [204,204,255]]/255;
  seablue = [[0,25,51]; [0,51,102]; [0,76,153]; [0,102,204]; [0,128,255]; [51,153,255]; [102,178,255]; [153,204,255]; [204,229,255]]/255;
  cyan    = [[0,51,51]; [0,102,102]; [0,153,153]; [0,204,204]; [0,255,255]; [51,255,255]; [102,255,255]; [153,255,255]; [204,255,255]]/255;
  seagreen = [[0,51,25]; [0,102,51]; [0,153,76]; [0,204,102]; [0,255,128]; [51,255,153]; [102,255,178]; [153,255,204]; [204,255,229]]/255;
  green   = [[0,51,0]; [0,102,0]; [0,153,0]; [0,204,0]; [0,255,0]; [51,255,51]; [102,255,102]; [153,255,153]; [204,255,204]]/255;
  lightgreen = [[25,51,0]; [51,102,0]; [76,153,0]; [102,204,0]; [128,255,0]; [153,255,51]; [178,255,102]; [204,255,153]; [229,255,204]]/255;
  yellow  = [[51,51,0]; [102,102,0]; [153,153,0]; [204,204,0]; [255,255,0]; [255,255,51]; [255,255,102]; [255,255,153]; [255,255,204]]/255;
  amber   = [[51,25,0]; [102,51,0]; [153,76,0]; [204,102,0]; [255,128,0]; [255,153,51]; [255,178,102]; [255,204,153]; [255,229,204]]/255;
  red     = [[51,0,0]; [102,0,0]; [153,0,0]; [204,0,0]; [255,0,0]; [255,51,51]; [255,102,102]; [255,153,153]; [255,204,204]]/255;

  % other initializations
  if (plot_en(2))
    normalize_en=1;  % to make f-domain plot look nice (use =1 for t-domain plot)
  else
    normalize_en=0;
  end


  if (impulse_shape==5) ovsamp=32; end       % set to max for rect. imp.

  if (plot_en(2)) ovsamp=32; end   % for maximal resolution at higher frequencies

  % time base
  Tsamp = Tsym/ovsamp; % for labeling x-axis: duration of one discrete time index (in seconds)

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, set up graphical output


  % Create an invisible figure.
  if (webdemo_en)
    fig = figure(1); set(fig, "visible", "off");
  else
    fig = figure(1); set(fig, "visible", "on");
  end


  landscape = "-S930,350";
  portrait  = "-S640,480";

  output_format = landscape;     % default


  set(0, "defaultlinelinewidth", 2);
%webdemo_en=1;
  if (webdemo_en)
    set(0, "defaultaxesfontsize", 8);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
    const_msize = 4;
  elseif (svg_en)
    set(0, "defaultaxesfontsize", 6);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 6);
    set(0, "defaulttextfontname", "Helvetica");
  set(0, "defaultlinelinewidth", 1);
    const_lw = 1;
    const_msize = 3;
  else
    set(0, "defaultaxesfontsize", 10);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
    const_msize = 2;
  end

  % END, set up graphical output
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




  % reset random variable generator
  if (rnd_data_en==0)
    rnd_val_gen_seed = (1:200); % arbitrary vector of length <=625 to init. rnd val generator to always same sequence
    rand("state", rnd_val_gen_seed);
  end



  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, get pulse used for pulse shaping

  t_stretch=1; % default, no time-stretch (i.e., no ISI)

  % select impulses
  if (impulse_shape==0)        % R.C., no stretch
    t_stretch=1;
  elseif (impulse_shape==1)   % R.C., with stretch to illustrate ISI
    t_stretch=1.4;
  elseif (impulse_shape==2)   % triangle, duration T, no ISI
    t_stretch=2;
  elseif (impulse_shape==3)   % triangle, duration 2T, no ISI
    t_stretch=1;
  elseif (impulse_shape==4)   % triangle, duration 3T, ISI
    t_stretch=2/3.0;
  elseif (impulse_shape==5)   % rectangular, duration T, no ISI
    t_stretch=1;
  end


  if (impulse_shape<2)
    [h_rc0, t_rc_n0, rc_av_pow0] = get_tdomain_imp_resp_h_raised_cosine(0, h_rc_sym_len, ovsamp, normalize_en, t_stretch);     % alpha=0, for power norm.
    [h_rc, t_rc_n, rc_av_pow] = get_tdomain_imp_resp_h_raised_cosine(h_rc_alpha, h_rc_sym_len, ovsamp, normalize_en, t_stretch);

    % NOTE: power of raised cosine imp. only normalized to 1 for alpha=0 (then spectrum rect.)
    % alpha>0 has greater power... need this factor to correctly plot analyt. ref. curve in spectrum
    rc_pow_corr_factor = rc_av_pow0/rc_av_pow;

  elseif (impulse_shape<5)  % triangular
    [h_rc, t_rc_n, rc_av_pow] = get_tdomain_imp_resp_h_triangle_rectangular(h_rc_sym_len, ovsamp, 1, normalize_en, t_stretch);
    rc_pow_corr_factor = 1;
  else % rect
    [h_rc, t_rc_n, rc_av_pow] = get_tdomain_imp_resp_h_triangle_rectangular(h_rc_sym_len, ovsamp, 0, normalize_en, t_stretch);
    rc_pow_corr_factor = 1;
  end

  h_rc_len = length(h_rc);

  % END, get pulse used for pulse shaping
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, generate impulse drain

  [s_re, rampup_len, rampdown_len, u_re] = generate_t_domain_pulse_drain(Nsym, nb_bits_per_real_comp, unipolar_en, ovsamp, h_rc, rnd_data_en, exclude_rampup_en, normalize_PAM_en);

  if (complex_en==0)
    s=s_re;
  else
    [s_im, rampup_len, rampdown_len, u_im] = generate_t_domain_pulse_drain(Nsym, nb_bits_per_real_comp, unipolar_en, ovsamp, h_rc, rnd_data_en, exclude_rampup_en, normalize_PAM_en);
    norm=1.0/sqrt(2);  % due to complex constellation
    s = norm*(s_re+j*s_im);
  end

  s_re=real(s);
  s_im=imag(s);

  s_len = length(s);

  % END, generate impulse drain
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%












  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, channel

  s_re_orig = s_re;
  s_im_orig = s_im;


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, add noise

%  SNR_dB=20;

  if (SNR_dB<100)
    SNR_lin=10^(SNR_dB/10.0);
    sn=1/sqrt(2*SNR_lin);       % noise per real component; assume signal power normalized to 1; complex constellations
    n=sn*(randn(1,s_len)+j*randn(1,s_len));
    s=s+n;
  end

  % END, add noise
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, SFO, resample
  % note: ovsamp=2 is sufficient as input; 1 is not sufficient, as alpha could be 1.0, i.e., BW = 2Rs then...
  % ToutOverTin=10001.0/10000 means, that output symbol duration is a bit larger, i.e., takes 1 symbol sample more every 10000 input symbol samples

%  ToutOverTin=1;
%  ToutOverTin=10001.0/10000;  % SFO; fclk_ppm=40ppm = 40/1e6=4e-5; thus, would be ToutOverTin=1+fclk_ppm; very small... (CHECK!)

  if (ToutOverTin!=1)

    len_s=length(s)

    nb_adj_imp=5;              % consider in total 2*nb_adj_imp +1 impulses for resampling, 2*3+1=7 should be fine, 2*5+1=11 is too much...

    s_resampled=zeros(1,len_s);

    for k_out=(0:1:len_s-1)
     amp=0;
     k_in=round(k_out*ToutOverTin); % determine discrete time center-index for relevant impulses to overlap
     for kk=(-nb_adj_imp:1:nb_adj_imp)
        if ((k_in+kk>=0) && (k_in+kk<len_s))
          amp=amp+s(k_in+kk+1).*sin(pi*(k_out*ToutOverTin-(k_in+kk))+1e-20)./(pi*(k_out*ToutOverTin-(k_in+kk))+1e-20);
        end
      end
      s_resampled(k_out+1)=amp;
    end

    s=s_resampled;

  end

  % END, SFO, resample
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



  % IQ-phase offset
  if (phi_off!=0)
    s = s .* exp(j*phi_off);
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, frequency offset

%  foff_over_Rs=0.0001;

  if (foff_over_Rs!=0)
    DeltaPhi=2*pi*foff_over_Rs;
    phi_k = (0:1:s_len-1)*DeltaPhi/ovsamp;
    s = s .* exp(j*phi_k);
  end



  % END, frequency offset
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, IQ imbalance

%  phi_imb=pi/180*0;   % Q ampl. imbalance
%  a_imb=1.0;       % Q phase imbalance
%  dc_imb=0;

  if ((phi_imb!=0) || (a_imb!=1) || (dc_imb!=0))
    x_re=real(s);
    x_im=imag(s);

    s_re = x_re - a_imb*sin(phi_imb)*x_im;
    s_im = 0    + a_imb*cos(phi_imb)*x_im + dc_imb;

    s = s_re+j*s_im;
  end

  % END, IQ imbalance
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, sampling time offset

 % t_samp_off=0;

  if (t_samp_off>0)
    s = [zeros(1,t_samp_off), s(1:end-t_samp_off)];
  end

  % END, sampling time offset
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  s_re=real(s);
  s_im=imag(s);

  % END, channel
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%








  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, generate plots
  % note: always generate plots from main function

  clf
  hold on

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot time domain signal (real part only)
  if plot_en(1)
    s_len = length(s_re);
    xlim([0, (s_len-1)*Tsamp]);
    ymin=-2.5;
    ymax=2.5;

    if (Nsym<=16)
      set(gca,'XTick',(-100:1:100)*Tsym);
      ymin=-1.8; ymax=1.8;
      grid on;
    end;

    if (normalize_PAM_en==0)
      P_av = get_PAM_average_power(nb_bits_per_real_comp, unipolar_en);
      ymin=ymin*sqrt(P_av);
      ymax=ymax*sqrt(P_av);
    end

    ylim([ymin, ymax]);
    xlabel("time t in sec.");
    ylabel("discrete-time sample-value");

%    [dummy_y, y_ampl_values] = get_random_PAM_vector(nb_bits_per_real_comp, unipolar_en, 1);
%    set(gca,'YTick',[y_ampl_values,-1,1]);


    if (exclude_rampup_en==0)
      rectangle("Position",[0,ymin,(rampup_len-1)*Tsamp,ymax-ymin],"FaceColor", magenta(7,:),"EdgeColor", magenta(7,:));              % mark filter ramp-up
      rectangle("Position",[(s_len-rampdown_len)*Tsamp, ymin,(s_len-1)*Tsamp,ymax-ymin],"FaceColor",magenta(7,:),"EdgeColor",magenta(7,:));  % mark filter ramp-down
    end


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%% for illustration, plot individual impulses
    if (show_indiv_imp_en)

      if (exclude_rampup_en)
        t_offset = -floor(length(h_rc)/2);
      else
        t_offset = 0;
      end

      RGBpalette=jet(Nsym);    % returns Nsym x 3 matrix
      for t_n = (1:Nsym)
        RGBcolor=RGBpalette(t_n,:);          % get t_n's row, i.e., returns 1x3 vector
        plot(  (((t_n-1)*ovsamp:h_rc_len-1+(t_n-1)*ovsamp)+t_offset)*Tsamp , u_re(t_n)*h_rc, "o-", "markersize", 2, 'linewidth', 1, "color", RGBcolor);
      end
    end

    plot((0:s_len-1)*Tsamp, s_re, "o-", "markersize", 1, 'linewidth',3, "color", "b");            % oversampled signal
    stem((rampup_len:ovsamp:rampup_len+Nsym*ovsamp-1)*Tsamp, u_re, 'linewidth',2, "markerfacecolor", "r", "color", "r");


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot spectrum
  elseif plot_en(2)

    [H_magn_sq, f_axis, N_disp_fbins, N_nyq_fbins] = compute_spectrum(s, Nsym, ovsamp, spectrum_display_bw_factor);

    fsamp = 1/(N_nyq_fbins*Tsamp*ovsamp);   % sampling frequency for x-axis labeling

    % set up plot output
    xlabel("frequency f in Hz");
    xlim([-N_disp_fbins/2, N_disp_fbins/2]*fsamp);       % fix to relevant spectral portions

    if (spectrum_dB_en)
      % note, then H_magn_sq_en=1 (forced)

      ylabel("squared spectral magnitude |G(f)|^2 in dB");

      rectangle("Position",[-N_disp_fbins/2*fsamp,-60, (N_disp_fbins-N_nyq_fbins)/2*fsamp,80],"FaceColor",magenta(7,:),"EdgeColor",magenta(7,:));
      rectangle("Position",[+N_nyq_fbins/2*fsamp,-60,  (N_disp_fbins-N_nyq_fbins)/2*fsamp,80],"FaceColor",magenta(7,:),"EdgeColor",magenta(7,:));

      if (complex_en==0)
        rectangle("Position",[-N_nyq_fbins/2*fsamp,-60, N_nyq_fbins/2*fsamp,80],"FaceColor",magenta(9,:),"EdgeColor",magenta(9,:));
      end

      ylim([-60, 20]);

      y_plot = 10*log10(H_magn_sq);

      plot(f_axis*fsamp, y_plot, 'linewidth',2, "b");
    else
      if (H_magn_sq_en)  ylabel("squared spectral magnitude |G(f)|^2");
                   else  ylabel("spectral magnitude |G(f)|");    end

      rectangle("Position",[-N_disp_fbins/2*fsamp,-0.5, (N_disp_fbins-N_nyq_fbins)/2*fsamp,2.5],"FaceColor",magenta(7,:),"EdgeColor",magenta(7,:));
      rectangle("Position",[+N_nyq_fbins/2*fsamp,-0.5, (N_disp_fbins-N_nyq_fbins)/2*fsamp,2.5],"FaceColor",magenta(7,:),"EdgeColor",magenta(7,:));

      if (complex_en==0)
        rectangle("Position",[-N_nyq_fbins/2*fsamp,-0.5, N_nyq_fbins/2*fsamp,2.5],"FaceColor",magenta(9,:),"EdgeColor",magenta(9,:));
      end

      ylim([-0.5, 2]);

      if (H_magn_sq_en)
        y_plot = H_magn_sq;
      else
        y_plot = H_magn_sq .^ (0.5);
      end

      plot(f_axis*fsamp, y_plot, 'linewidth',2, "b");
    end

    % also plot theoret. curve
    H_ref_magn_sq = get_fdomain_imp_resp_H_raised_cosine(h_rc_alpha, N_nyq_fbins/2) .^ 2;
    H_ref_magn_sq = H_ref_magn_sq * rc_pow_corr_factor;
    H_ref_len = length(H_ref_magn_sq);
    if (spectrum_dB_en)
      y_plot = 10*log10(H_ref_magn_sq);
      plot((-floor(H_ref_len/2):floor(H_ref_len/2))*fsamp, y_plot, 'linewidth', 4, "r");
    else
      if (H_magn_sq_en)   y_plot = H_ref_magn_sq;  else   y_plot = H_ref_magn_sq .^ (0.5);   end
      plot((-floor(H_ref_len/2):floor(H_ref_len/2))*fsamp, y_plot, 'linewidth', 4, "r");
    end

    output_format = portrait;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%  plot phase diagram ("Ortskurve") for illustrating offset QAM (mainly OQPSK)
  elseif plot_en(3)

    xlim([-1.4,1.4]);
    ylim([-1.4,1.4]);


    axis ("square");

    xlabel("real, inphase component");
    ylabel("imaginary, quadrature component");
%    plot(s_re, s_im, "color", "r");

    if (offset_QAM_en)
      t_offset = ovsamp/2;
    else
      t_offset = 0;
    end

    plot(s_re(1:end-t_offset), s_im(1+t_offset:end), "color", "b"); % to illustrate offset-QPSK
    output_format = portrait;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot individual pulse
  elseif plot_en(4)
    xlim([t_rc_n(1), t_rc_n(end)]);
    set(gca,'XTick',(-100:1:100));

    grid on
    xlabel("normalized time t/T");
    ylabel("discrete-time sample-value");
    plot(t_rc_n, h_rc, 'linewidth',2, "o-", "markersize", 2, "color", "b");   % for debugging

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot eye diagram (real part only)
  elseif plot_en(5)
    s_len = length(s_re);
    xlim([-ovsamp, ovsamp]*Tsamp/Tsym);
    ylim([-2.5, 2.5]);
    xlabel("normalized time t/T");
    ylabel("amplitude");

    set(gca,'XTick',(-10:0.5:10));
    set(gca,'YTick',(-4:0.5:4));
%    [dummy_y, y_ampl_values] = get_random_PAM_vector(nb_bits_per_real_comp, unipolar_en, 1);
%    set(gca,'YTick',[y_ampl_values,-1,1]);

    grid on

    for (k=1:Nsym-2)
      TwoT_idx = ((k-1)*ovsamp+1 : (k-1)*ovsamp+1+2*ovsamp);
      plot( ((0:2*ovsamp)*Tsamp-Tsym)/Tsym, s_re(TwoT_idx), "o-", "markersize", 1, 'linewidth',3, "color", "b");
    end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot scatter-plot (for a complex constellation)
  elseif plot_en(6)
    s_len = length(s_re);

    xmax=1.5;
    ymax=1.5;

    xlim([-xmax, xmax]);
    ylim([-ymax, ymax]);
    xlabel("I-channel (re)");
    ylabel("Q-channel (im)");

%    set(gca,'XTick',(-10:0.5:10));
%    set(gca,'YTick',(-4:0.5:4));
%    [dummy_y, y_ampl_values] = get_random_PAM_vector(nb_bits_per_real_comp, unipolar_en, 1);
%    set(gca,'YTick',[y_ampl_values,-1,1]);

    grid on


    % sampled
    x = s_re((0:1:Nsym-1)*ovsamp+1);
    y = s_im((0:1:Nsym-1)*ovsamp+1);


%%%%    plot(x,y, "o", "markersize", 1, 'linewidth',1, "color", "b");

    NbBinsHalf=50;
    NbBins=2*NbBinsHalf+1;         % nb of bins per dimension

    x_ax=(-NbBinsHalf:1:NbBinsHalf)/NbBinsHalf*xmax;
    y_ax=x_ax;

    h2d=zeros(NbBins,NbBins);

%    xidx=round(x/xmax*NbBinsHalf+NbBinsHalf);
%    yidx=round(y/xmax*NbBinsHalf+NbBinsHalf);
    xidx=ceil(x/xmax*NbBinsHalf+NbBinsHalf);
    yidx=ceil(y/xmax*NbBinsHalf+NbBinsHalf);

    xidx=min(xidx,NbBins);
    xidx=max(xidx,1);

    yidx=min(yidx,NbBins);
    yidx=max(yidx,1);

    % fill histogram
    for (k=1:Nsym)
      xi=xidx(k);
      yi=yidx(k);
      h2d(xi,yi)=h2d(xi,yi)+1;
%      h2d(xi,yi)=k;          % for plot showing latest symbol
    end


    if (colormap_idx==0)       % jet, light
      % still in testing mode: new color map, good for printing (white background), colorful like jet
      nb_cols=256;
      newmap= jet(nb_cols);
      % generate smooth transition of jet() from "black to white" brightness
      for (k=1:nb_cols)
        newmap(k,1:3) = (newmap(k,1:3)*(k-1) + [nb_cols-k, nb_cols-k, nb_cols-k])/(nb_cols-1);
      end
      jet_light_inue=newmap;
      colormap( jet_light_inue );
    elseif (colormap_idx==1) % HELIX colorbar
      nb_cols=256;
      nlev=nb_cols;
      start=1.0; rots=-1.5; hue=2; gamma_factor=1.0;
      newmap=CubeHelix(nlev,start,rots,hue,gamma_factor);
      newmap=flipud(newmap);
      % generate smooth transition from "gray to white" brightness
      gray_level=0.2;
      for (k=1:nb_cols)
        newmap(k,1:3) = newmap(k,1:3)*(1-gray_level)+[gray_level,gray_level,gray_level];
      end
      colormap( newmap );
    elseif (colormap_idx==2)           % gray
      colormap( flipud(gray(256)));
    elseif (colormap_idx==3) % jet, white
      jet_white = jet(256);
      jet_white(1:24,1) = linspace(1,0,24); % red
      jet_white(1:24,2) = linspace(1,0,24); % green
      jet_white(1:24,3) = linspace(1,1,24); % blue
      colormap(jet_white);
    else  % (colormap_idx==4)       % jet, normal
      colormap( jet(256));
    end

    imagesc(x_ax, y_ax, h2d'/Nsym);
    set(gca,'ydir','normal');
    axis("square");

    if (colorbar_en)
      colorbar("EastOutside");
      title(colorbar,"prob. of occurrence");
    end

%plot( x_ax(xidx(1:Nsym/50:Nsym)), y_ax(yidx(1:Nsym/50:Nsym)), "o", "markersize", 1, 'linewidth',3, "color", "b"); % plot every 50th point; not so insightful...


    % compute EVM, sampled
    x_orig = s_re_orig((0:1:Nsym-1)*ovsamp+1);
    y_orig = s_im_orig((0:1:Nsym-1)*ovsamp+1);
    MSE_dB = 10*log10(sum((x-x_orig).^2 + (y-y_orig).^2)/Nsym+1e-10)

    title(["EVM=" num2str(MSE_dB) "dB"]);

    output_format = portrait;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot auto-corr
  elseif plot_en(7)
  max_time_lag=16;
    xlim([-max_time_lag,max_time_lag]);
%    ylim([-1,1]);
%    set(gca,'XTick',(-100:1:100));

    grid on
    xlabel("time lag l");
    ylabel("auto-correlation");


    acorr=compute_auto_correlation_normalized(s, max_time_lag, Nsym, ovsamp);
%    plot((0:1:max_time_lag), acorr, "o-", "markersize", 1, 'linewidth',3, "color", "m");            % oversampled signal
%    stem((-max_time_lag:1:max_time_lag), [ abs(acorr(max_time_lag+1:-1:2)) abs(acorr) ], 'linewidth',2, "markerfacecolor", magenta(8,:), "markersize", 3, "color", magenta(3,:));
    stem((-max_time_lag:1:max_time_lag), [ abs(acorr) ], 'linewidth',2, "markerfacecolor", magenta(8,:), "markersize", 3, "color", magenta(3,:));
    output_format = portrait;
  end

  hold off

  % END, generate plots
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, output to file

  print(filename, "-dpng",output_format);

  if (svg_en)
    print(["fig_",filename], "-dsvg",output_format);
  end

  % END, output to file
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


end

% END, main program (control)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%














%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, function definitions


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [P_av]=get_PAM_average_power(nb_bits, unipolar_en)
  M = 2^nb_bits;                          % nb of ampl. levels

  if (unipolar_en)
    ampl_vec = (0:1:M-1);                 % M ampl. levels (unipolar)
  else
    ampl_vec = (-M+1:2:M-1);              % M ampl. levels (bipolar)
  end

  P_av = 1/M*sum(ampl_vec.^2);            % NORMALIZE to av. power 1 at discrete time samples
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y,ampl_vec]=get_random_PAM_vector(nb_bits, unipolar_en, N, normalize_PAM_en)
  M = 2^nb_bits;                          % nb of ampl. levels

  if (unipolar_en)
    ampl_vec = (0:1:M-1);                 % M ampl. levels (unipolar)
  else
    ampl_vec = (-M+1:2:M-1);              % M ampl. levels (bipolar)
  end

  prob_vec = ones(1,M)*1/M;               % all equally probable

  if (normalize_PAM_en==1)
    P_av = get_PAM_average_power(nb_bits, unipolar_en);   % NORMALIZE to av. power 1 at discrete time samples
    ampl_vec=ampl_vec/sqrt(P_av);
  end

  y=discrete_rnd(ampl_vec, prob_vec, 1, N);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% note: len is nb of taps without considering oversampling
function [h, t_n, av_pow] = get_tdomain_imp_resp_h_raised_cosine(alpha, len, ovsamp, normalize_en, t_stretch)

  if (mod(len,2)==0)
    t_n = (-floor(len/2):1/ovsamp:floor(len/2)-1/ovsamp);
  else
    t_n = (-floor(len/2)-1.0/2:1/ovsamp:ceil(len/2)-1.0/2-1/ovsamp);
  end

  h = zeros(1,length(t_n));

  k=1;

  for n=t_n * t_stretch;
    numerator   = sinc(n)*cos(pi*alpha*n);
    denominator = 1-(2*alpha*n)^2;

%    if (denominator==0) denominator=1e-5; end
    if (denominator==0)
      numerator = sinc(1/2/alpha)*pi/4;    % limit
      denominator=1;
    end


    h(k) = numerator/denominator;
    k=k+1;
%%    h = sinc(t_n).*cos(pi.*alpha.*t_n)./(1-(2.*alpha.*t_n).^2);
%%    h( t_n == inf ) = sinc(1./(2.*alpha)).*pi./4; % for t_n = 1 / (2*alpha) formula for h_t_n yields 0/0 --> bring to limit value
  end

% partial response impulses; test
%h = h + [zeros(1,ovsamp), h(1:end-ovsamp)];  % PR class 1, duobinary
%h = h - [zeros(1,2*ovsamp), h(1:end-2*ovsamp)];  % PR class 4
%h = h - [zeros(1,4*ovsamp), h(1:end-4*ovsamp)];  % PR class 6

  % normalize
  av_pow = sum(h .^ 2);

  if (normalize_en)

    h = h/sqrt(av_pow);
  end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% note: len is nb of taps without considering oversampling
function [h, t_n, av_pow] = get_tdomain_imp_resp_h_triangle_rectangular(len, ovsamp, rect_trian_en, normalize_en, t_stretch)

  if (mod(len,2)==0)
    t_n = (-floor(len/2):1/ovsamp:floor(len/2)-1/ovsamp);
  else
    t_n = (-floor(len/2)-1.0/2:1/ovsamp:ceil(len/2)-1.0/2-1/ovsamp);
  end

  h = zeros(1,length(t_n));

  k=1;

  denominator=1;

  for n=t_n * t_stretch;
    numerator=0;

    if (rect_trian_en==0)
      % rectangular
      if (n==-0.5) numerator=0.5
      elseif ((n>-0.5) && (n<0.5)) numerator=1;
      elseif (n==+0.5) numerator=0.5
      end
    else
      % triangle
      if ((n>-1) && (n<0))  numerator=(n+1);
      elseif ((n>=0) && (n<1)) numerator=(-n+1);
      end
    end

    h(k) = numerator/denominator;
    k=k+1;
  end


  % normalize
  if (normalize_en)
    av_pow = sum(h .^ 2);
    h = h/sqrt(av_pow);
  end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% note: len is nb of taps without considering oversampling
%function [h, t_n, av_pow] = get_tdomain_imp_resp_h_raised_cosine(alpha, len, ovsamp, normalize_en)
function [h, t_n, av_pow] = get_tdomain_xxxxx_imp_resp_h_raised_cosine(alpha, len, ovsamp, normalize_en)

  if (mod(len,2)==0)
    t_n = (-floor(len/2):1/ovsamp:floor(len/2)-1/ovsamp);
  else
    t_n = (-floor(len/2)-1.0/2:1/ovsamp:ceil(len/2)-1.0/2-1/ovsamp);
  end

  h = zeros(1,length(t_n));

  k=1;
  for n=t_n
    if ((n>=-1/2) & (n<1/2))
      h(k) = 1;
    else
      h(k) = 0;
    end

    k=k+1;
  end


  % normalize
  if (normalize_en)
    av_pow = sum(h .^ 2);
    h = h/sqrt(av_pow);
  end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% normalized?
function H_f_n = get_fdomain_imp_resp_H_raised_cosine(alpha, ovsamp)

  f_n = (0:1/ovsamp:3)/2/pi;

  % compute H(f_n)
  % first positive frequency axis (0 <= f_ n<= 3)
  H_f_n = (cos(pi./4./alpha.*(2*pi*f_n-1+alpha))).^2;
  H_f_n (f_n <=(1-alpha)/2/pi) = 1;
  H_f_n (f_n >= (1+alpha)/2/pi) = 0;
  % add negative axis
  H_f_n = [H_f_n(end:-1:2), H_f_n];
  f_n = [-f_n(end:-1:2), f_n];
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y, rampup_len, rampdown_len, u] = generate_t_domain_pulse_drain(Nsym, nb_bits_per_real_comp, unipolar_en, ovsamp, h_rc, rnd_data_en, exclude_rampup_en, normalize_PAM_en)

  % generate random BPSK data
  [u, dummy_ampl] = get_random_PAM_vector(nb_bits_per_real_comp, unipolar_en, Nsym, normalize_PAM_en);

  % for small nb of PAM symbols, impose constant sign pattern for better illustration
  if ((Nsym<6) && (rnd_data_en==0))
    u = [1,1,-1,1,-1,-1](1:Nsym) .* abs(u);
%    u = [1,1,0,0,0,0](1:Nsym) .* abs(u);  % PR1
%    u = [1,0,-1,0,0,0](1:Nsym) .* abs(u);  % PR4
%    u = [1,0,0,-1,0,0](1:Nsym) .* abs(u);  % PR6
  end

  u_up = upsample(u,ovsamp);

  % remove trailing zeros
  u_up = u_up(1:end-(ovsamp-1));

  rampup_len   = floor(length(h_rc)/2);
  rampdown_len = rampup_len-1;

%printf("rampup_len = %hd\n",rampup_len);

  y = conv(u_up, h_rc);                % convolution in t-domain

  if (exclude_rampup_en)
    y = y(1+rampup_len:end-rampdown_len);
    rampup_len=0;
    rampdown_len=0;
  end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y_magn_sq, f_axis, N_disp_fbins, N_nyq_fbins] = compute_spectrum(s, Nsym, ovsamp, spectrum_display_bw_factor)

  s_len = length(s);

%sum(s .* conj(s))  % power ok

  N_nyq_fbins = 101;  % fix nb of freq. bins in Nyq. bandwidth


  % Note:
  % - each frequency bin has bandwidth (1/(Ts/ovsamp))/s_len; if Ts is the PAM symbol duration, then the Nyquist bw is 1/Ts
  % - fix to  (1/Ts)/N_nyq_fbins *1/2  (*1/2, since a total bandwidth of 2 Ny BW should be displayed)
  % - for this: 1. append zeros prior to fft such that s_len/ovsamp is integer mult. of 2*N_nyq_fbins = N_disp_fbins

  N_disp_fbins = spectrum_display_bw_factor*N_nyq_fbins;   % display a spectrum of "display_bw_factor"-times the Nyq-BW

  N_append_zeros = rem(N_disp_fbins-rem(s_len/ovsamp, N_disp_fbins),N_disp_fbins) * ovsamp;

  s = [s, zeros(1, N_append_zeros)];


  s_len = length(s);

  y = fft(s)/sqrt(ovsamp*Nsym);    % wrong: sqrt(s_len); correct: sqrt(ovsamp*Nsym)
  y_magn_sq = real(y .* conj(y));

%sum(y .* conj(y))   % power ok

  % re-sort neg.freq. for plotting
  y_magn_sq = [y_magn_sq(1+floor(s_len/2):s_len), y_magn_sq(1:floor(s_len/2))];

  %             2. after fft, sum up, each, energy of (s_len/ovsamp)/N_nyq_fbins freq. bins

  N_sum_bins = (s_len/ovsamp)/N_nyq_fbins;
  N_tot_fbins = s_len/N_sum_bins;         % total nb of freq. bins after summation

  % sum over freq. bins
  for new_bin = 0:N_tot_fbins-1-1
%    y_magn_sq(1+new_bin) = sum(y_magn_sq( 1+new_bin*N_sum_bins : 1+new_bin*N_sum_bins+N_sum_bins-1 ))/N_sum_bins;
    % NOTE: done this way (average of neighboring bin-sums over 2) to get sym. spectrum for s real
    y_magn_sq(1+new_bin) = (sum(y_magn_sq( 1+new_bin*N_sum_bins : 1+new_bin*N_sum_bins+N_sum_bins-1 ))/N_sum_bins
                           +sum(y_magn_sq( 1+1+new_bin*N_sum_bins : 1+1+new_bin*N_sum_bins+N_sum_bins-1 ))/N_sum_bins ) /2;
  end

  y_magn_sq = y_magn_sq(1:N_tot_fbins);  % set new length

%N_sum_bins
%y_magn_sq = y_magn_sq/16;

%sum(y_magn_sq)    % power ok

  f_axis = ((-floor(N_tot_fbins/2):ceil(N_tot_fbins/2)-1) + 0.5); % +0.5 due to averaging of neighboring bins

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Usage:
%         colormap(CubeHelix(...))
%
%==========================================================
% Calculates a "cube helix" colour map for MATLAB. The
% colours are a tapered helix around the diagonal of the
% RGB colour cube, from black [0,0,0] to white [1,1,1].
% Deviations away from the diagonal vary quadratically,
% increasing from zero at black to a maximum and then
% decreasing to zero at white, all the time rotating in
% colour.
%
% The input values are:
%   nlev  = number of colour steps
%   start = colour to begin at (1=red, 2=green, 3=red;
%           e.g. 0.5=purple)
%   rots  = number of rotations
%   hue   = hue intensity scaling, 0=B&W
%   gamma = intensity correction
%
% The routine returns an nlev-by-3 matrix that can be used
% as a colourmap for MATLAB (function 'colormap').
%
% Use (256,0.5,-1.5,1.2,1.0) as defaults.
%
% See arXiv:1108.5083 for more details.
%
%----------------------------------------------------------
% Original written in Fortran77 by Dave Green, 2011 Jan 10
% Transcribed to MATLAB by Philip Graff, 2011 Sept 1
%
% https://www.mrao.cam.ac.uk/~dag/CUBEHELIX/
% https://www.mrao.cam.ac.uk/~dag/CUBEHELIX/cubetry.html
%
%==========================================================
function map = CubeHelix(nlev,start,rots,hue,gamma)
  map=zeros(nlev,3);
  A=[-0.14861,1.78277;-0.29227,-0.90649;1.97294,0];
  for i=1:nlev
      fract=(i-1)/(nlev-1);
      angle=2*pi*(start/3+1+rots*fract);
      fract=fract^gamma;
      amp=hue*fract*(1-fract)/2;
      map(i,:)=fract+amp*(A*[cos(angle);sin(angle)])';
      for j=1:3
          if map(i,j)<0
              map(i,j)=0;
          elseif map(i,j)>1
              map(i,j)=1;
          end
      end
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y_vec = compute_auto_correlation_normalized(x, max_lag, Nym, ovsamp)
  len = length(x);
  norm=sum(x.*conj(x));    % so that we get a normalized auto-corr.

  y_vec = zeros(1,2*max_lag+1);
  for t_lag = (-max_lag:1:max_lag)
    if (t_lag>=0)
      x1 = [ x, zeros(1,t_lag) ];
      x2 = [ zeros(1,t_lag), x ];
    else
      x1 = [ zeros(1,-t_lag), x ];
      x2 = [ x, zeros(1,-t_lag) ];
    end
% MUSS CORR PER Ts machen, und toff muss man sehen; keine full AC... konfiguierbar machen! in Zug
% sampled
%x1 = x1((0:1:Nsym-1)*ovsamp+1);
%x2 = x2((0:1:Nsym-1)*ovsamp+1);

    y_vec(1+max_lag+t_lag)=sum(x1 .* conj(x2))/norm;
%y_vec(1+t_lag)=10*log10(abs(y_vec(1+t_lag)));

  end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%function y_vec = compute_diff_metric(x, max_lag, ovsamp)
%%  len = length(x);
%%  norm=sum(x.*conj(x));    % so that we get a normalized auto-corr.
%%
%%  y_vec = zeros(1,max_lag);
%%  for t_lag = (0:max_lag-1)
%%    for k=(1:len)
%%
%%    end
%%
%%    x1 = [ x, zeros(1,t_lag) ];
%%    x2 = [ zeros(1,t_lag), x ];
%%
%%    y_vec(1+t_lag)=sum(x1 .* conj(x2))/norm;
%%
%%
%%  end
%%
%%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



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