% decrease OVSAMP

% based on "pulse_shaping.m"
% - shortened, tailored
% to add
% - fshift by explicit freq.
% - nomenclature, other signals, phi, cos, exp
% - configurable FILTER (R.C., via roll-off)
% - later: result after demod., with/without filter, roll-off... etc.
% - BW of test signal...
% - animation... how it comes up ...
%
% - t, f, t/f-domain display
% - fshift-animation
%
% - for ÜT1, move to complex signals...
% .. SSB vs DSB; Hilbert filter; notation...

% checks
% - re-check normalizations; t-/f-domain...
% - reduce oversampling...
% - idea: indicate redundant parts in spectrum
% note: with carrier, one symbol only... makes no sense, edge effects in spectrum
% - idea: select what signal to plot (pull down); then t/f-domain
% - idea: compute difference to transmitted BB-signal as MSE in dB



% t-domain:     modulation(10000,2.5e3, 1 ,0 ,1,    2, 0, 0, 0.1, 0, 50, 1    ,8,0, 1, 3 ,0,0,"dummy")
%               modulation(1000,2.5e3, 2 , 0  ,1,    2, 1,   1.5, 0.1, 0, 50, 0    ,8,0,1, 3 ,3 ,0,0,"dummy")
% in lect., AM2 (demod.): modulation(2000,2.5e3, 2 , 0, 0  ,1,    2,    4  ,   1.5, 0.1, 0, 50, 0   ,8,0,1,    3    ,3 ,0,0,"dummy")
% .. modulation(2000,2.5e3, 10e3 , 0, 0, 200,   10e3,    0  ,   1.5, 0.1, 0, 50, 0   ,8,0,1,    0    ,3 ,0,0,"dummy")
% .. modulation(2000,2.5e3, 1, 10e3 , 0, 0 ,   20   ,   10e3,    0  ,   1.5, 0.1, 0, 50, 1   ,8,0,1,    3    ,3 ,0,0,"dummy")
%
% WS2016
% old modulation(1,2.5e3, 2, 5e3 , 2, 0, 200, 5e3, 0, 1.5, 0.1, 0, 20, 1, 8, 0, 1, 0, 3, 0,0,"dummy")
% old modulation(1, 2.5e3, 0, 5e3 , 0, 0, 200, 5e3, 0, 1.5, 0.1, 0, 20, 1, 8, 0, 1, 0, 3, 0,0,"dummy")
% SS2017: audio signal test
% modulation(1, 8e3, 10, 16e3 , 0, 0, 200, 16e3, 0, 1.5, 0.1, 0, 2000, 0, 1, 0,0, 2, 101, 0, 1, 0, 3, 0,0,"dummy")
% modulation(1000, 4e3, 20, 16e3 , 0, 0, 200, 16e3, 0, 1.5, 0.1, 0, 200,1, 0, 0, 1, 0,0, 2, 101, 0, 1, 0, 3, 0,0,"dummy")
% modulation(1, 4e3, 0, 16e3 , 0, 0, 200, 16e3, 0, 1.5, 0.1, 0, 200,2, 0.003, 0.006, 1, 0,0, 100, 2, 101, 0, 1, 0, 3, 0,0,"dummy")

% BEGIN, main function (control)
function modulation(
  Nsym,                         % nb of PAM symbols
  fmax,                         % max freq. of BB signal in Hz
  source_signal_type,           % 0=real, BPSK, raised cos. pulse shaping (just as an example)
                                % 1=complex, QPSK, raised cos. pulse shaping
                                % 2=complex, analytical (im is +Hilbert-trafo of re), USB
                                % 3=complex, analytical (im is -Hilbert-trafo of re), LSB
                                % 10="analog" audio signal, read from file, at fs=16kHz, so fmax MUST BE SET to 8kHz at input (FORCED)
                                % 11="analog" another audio signal, read from file, at fs=16kHz
                                % 20=continuous sinusoid at 440Hz
                                % 21=continuous sinusoid at 1000Hz
                                % 22=sum of 440Hz and 1000Hz sinusoids
  tx_fc,                        % carrier frequency in Hz
  tx_mod,                       % 0=cos, 1=-sin, 2=exp(j..)
  carrier_en,                   % 0=AM without carrier; -1=AM with carrier, auto dc_offset; else carrier_en=dc_offset; -2: carrier only
  SNR_dB,                       % set to 200dB to switch off
  rx_fc,                        % receiver LO freq. in Hz
  rx_demod,                     % 0=cos; 1=-sin, 2=exp(j..), 3=sqr-law, 4=pos, 5=abs, 6=cos, -sin
  rx_filter_fmax_factor,        % 0=off, 1=enabled (in general BW rx_filter_fmax_factor*2*fmax); note: Rs=2fmax here; to show effect of wider LP filters
  rx_filter_alpha,              % 0..1 roll-off
  rx_phi,                       % default 0; DSB-AM has no output for pi/2+k*pi
  Nsym_display,                 % t-domain plot: display from middle section only Nsym_display<Nsym symbols
  tdomain_display_format,       % bit 1 set: [FORCED for single impulse; see below] then t=0 at middle of e.g. a single impulse (makes only sense if signal is a single impulse)
                                % bit 2 set: use display_t_min, display_t_max for time domain display (in sec)
  display_t_min,                % in sec., use if bit 2 set in tdomain_display_format
  display_t_max,                % in sec., use if bit 2 set in tdomain_display_format
  N_t_steps_spectrogram,        % 100; to increase spectral resolution, use small nb
  spectrum_display_format,      % bit 1 (1) set: show neg. freq. axis (should be on by default);
                                % bit 2 (2) set: spectrum in dB
                                % bit 3 (4) set: use display_f_min, display_f_max for spectral display
                                % bit 4=2^3=8 set: display t/f-uncertainty in spectrogram
  display_f_min,                % use if bit 3 of spectrum_display_format is set
  display_f_max,                % use if bit 3 of spectrum_display_format is set
  spectrum_display_bw_factor,   % display "spectrum_display_bw_factor"-times the Nyquist bandwidth R_s in spectrum
  N_nyq_fbins,                  % fix nb of freq. bins in Nyq. bandwidth top show in spectrum plot
  rnd_data_en,                  % 0=not random
  display_rx_filter_en,         % show rx filter in spectrum plot
  display_signal,               % 0=at TX, BB
                                % 1=at TX, after modulation, PB
                                % 2=at RX, after channel, PB
                                % 3=at RX, after demodulation, BB
                                % 4=at RX, after filter, BB
  td_fd_display_options,        % switch, subplots to be displayed;
                                % bit 1=1: t-domain
                                % bit 2=2: f-domain
                                % bit 3=4: spectrogram
                                % e.g. t-/f-domain: set this value to 1+2=3
  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
)

  warning_message_str="";


  t_display_ctr_zero_en = bitand(tdomain_display_format, 1);

  % if single t-domain impulse, force display t-axis to be symmetric to zero (impulse centered at t=0)
  if ((Nsym==1) && (source_signal_type<10))
    t_display_ctr_zero_en=1;
  end

t_display_ctr_zero_en=0; % does not work well yet; need a "symbol" and "absolut" time axis display option

  tdomain_display_tmin_tmax_en=bitand(tdomain_display_format,2);


  show_neg_freq_en=bitand(spectrum_display_format,1);              % show negative frequency axis (should be on by default)
  spectrum_dB_en=bitand(spectrum_display_format,2);                % plot spectrum y-axis in dB
  spectrum_display_fmin_fmax_en=bitand(spectrum_display_format,4);
  spectrogram_tf_uncertainty_en = (bitand(spectrum_display_format,8)==8);


  rx_filter_bw_en=rx_filter_fmax_factor;

  complex_bb_en=1;                      % NOTE: then tx_mod and rx_demod forced to Re exp(j...), exp(-j...)
  if (source_signal_type==0)  complex_bb_en=0;  end


  % special set-up for audio signals
  if ((source_signal_type>=10) && (source_signal_type<=12))
    complex_bb_en=0;
    if (fmax!=8000)
      warning_message_str="fmax forced to 8kHz for audio signals";
      fmax=8000;      % since audio signals are sampled at 16kHz
    end
  end

  if ((source_signal_type>=20) && (source_signal_type<=22)) complex_bb_en=0;  end


  % time-/frequ.base, all based on baseband bandwidth 2*fmax
  % derive dig. QAM specs for generating signals
  Rs=2*fmax;    % symbol rate = 2*fmax, if perfect sinx/x-pulse shaping used
  Tsym=1/Rs;    % symbol duration in seconds (set to 1 for normalization)








  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN check input parameters for consistency
  td_en = (bitand(td_fd_display_options, 1) == 1);
  fd_en = (bitand(td_fd_display_options, 2)  == 2);
  spectrogram_en = (bitand(td_fd_display_options, 4)  == 4);
  % END check input parameters for consistency
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  % based on pulse_shaping.m
  % gradually remove some of those variables
  nb_bits_per_real_comp=1;        % nb of bits per PAM symbol
  impulse_shape=0;                 % R.C., triangular or rect. impulse
  unipolar_en=0;                  % unipolar impulses also possible
%  Tsym=1;                         % symbol duration in seconds (set to 1 for normalization)
  offset_QAM_en=0;                % shift Q-channel by T/2 if enabled (for phase plot)

  % a little too high...
  ovsamp=64;                       % nb of samples in between PAM symbols (oversampling)
  if ((source_signal_type>=10) && (source_signal_type<=12))
    ovsamp=4;
  elseif ((source_signal_type>=20) && (source_signal_type<=22))
    ovsamp=8;
  end

  spectrum_display_bw_factor=min(spectrum_display_bw_factor,ovsamp);




  rx_filter_sym_len=64;                 % length of t-domain raised cosine filter in PAM-symbols
%  rnd_data_en=0;%1;                  % enable random data
  show_indiv_imp_en=0;            % show individual PAM-impulses in t-plot
  exclude_rampup_en=0;            % if 1, then don't include filter rampup in t-domain signal
  H_magn_sq_en=0;                 % plot |H(f)| (to better see raised cosine) or |H(f)|^2
  phi_off_deg=0;                      % for scatter plot
  foff_over_Rs=0;                 % for scatter plot, =1/100000;
  phi_imb_deg=0;                      % for scatter plot, =pi/180*0; %(-10);   % Q ampl. imbalance
  a_imb=1;                        % for scatter plot, =1; %1.2;       % Q phase imbalance
  dc_imb=0;                       % for scatter plot, =0; %0.1;
  t_samp_off=0;                   % for scatter plot, =0;
  colormap_idx=0;                 % for scatter plot


rx_filter_len=8; % BW is 1/Tsym plus roll-off alpha; *ov_samp



% for full IQ/QAM modulation!
% OUTPUT WARNING MESSAGE
if ((complex_bb_en==1) && (tx_mod!=2) && (display_signal>0))
%if ((complex_bb_en==1) && ((tx_mod!=2) || (rx_demod!=2)))
  tx_mod=2;        % exp(-j..., re)
%  rx_demod=2;     % exp(-j...)
  warning_message_str="complex BB signal, did set tx to exp(..j2pif_c t), Re()";
end




%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


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

  if (rem(N_nyq_fbins,2)==0) N_nyq_fbins=N_nyq_fbins+1; end  % make odd

  if ((N_nyq_fbins<51) || (N_nyq_fbins>1001))
    N_nyq_fbins = 101;  % fix nb of freq. bins in Nyq. bandwidth
  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 (rx_filter_alpha>1)  rx_filter_alpha=1;    warning("setting  'rx_filter_alpha = 1'");     end
  if (rx_filter_alpha<0)  rx_filter_alpha=0;    warning("setting  'rx_filter_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 (rx_filter_alpha==0) rx_filter_alpha = 0.00001; end
  if (rx_filter_alpha==1) rx_filter_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_bb_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 (fd_en)
    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=64; 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";

%  landscape = "-S1860,700";
%  portrait  = "-S1280,960";

%  output_format = landscape;     % default
  output_format = portrait;

  % new
  def_lw_thin = 1;
  def_lw = 2;
  def_lw_thick = 4;


  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


  h_rc_alpha=0.1;
  h_rc_sym_len=64;

  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
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% RX filter,  normalized to 0dB gain
rx_filter_t_stretch=rx_filter_bw_en; % control BW here.. times R_s
[h_rx_filter, h_rx_filter_k, h_rx_filter_pow] = get_tdomain_imp_resp_h_raised_cosine(rx_filter_alpha, rx_filter_sym_len, ovsamp, normalize_en, rx_filter_t_stretch);


hilbert_filter_len=16;       % seems to be enough
[h_hilbert_filter, h_hilbert_filter_av_pow] = get_tdomain_imp_resp_h_hilbert(hilbert_filter_len, ovsamp, 1);


% comms link, start


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, generate impulse drain
  %
  % impulse drain: just used as source signal at the moment


  % first check: is source an audio signal?
  if ((source_signal_type>=10) && (source_signal_type<=12))

    if (source_signal_type==10) wavefile="man1_wb.wav"; end
    if (source_signal_type==11) wavefile="woman1_wb.wav"; end
    if (source_signal_type==12) wavefile="piano_beethoven_fmax8kHz.wav"; end

    [w_y,w_fs,w_nbits]=wavread(wavefile);

    if (w_fs != 2*fmax)
      warning_message_str="wavefile: (w_fs != 2*fmax)";
    end
    % sound(w_y, w_fs, w_nbits);  % works
    % at 16kS/s; set fmax to 8kHz
%  size(w_y)
%% %  w_fs
%% %  w_nbits
% cut short
%%  s=s(1:tmp_len);
    s=w_y(:,1)';
    s=s(1:length(s)/2); % just cut short for now, otherwise too long...
    voice_activity_factor=0.3; % FUDGE FACTOR... NOT CLEAR YET; HOW TO PROPERLY SET THE NORMALIZATION FOR AUDIO SIGNAL?

      % normalize
      s_active_len=length(s)*voice_activity_factor;
      pow_sqrt=sqrt(sum(s.^2)/s_active_len);
      s=s./pow_sqrt;

%  size(s)
% w_t=(0:1:length(s)-1).*1/w_fs;
% plot(w_t, s, 'linewidth', def_lw_thick, "m");
%return;
    tmp_ovsamp=(1/w_fs) / Tsamp;  % if fmax=8kHz, then tmp_ovsamp==ovsamp
%  tmp_ovsamp
    s=upsample(s,tmp_ovsamp);
    s=conv(s,h_rc);          % LPF, fmax MUST BE fmax=8000Hz, since w_fs=16kHz

    % how many Nsym would that correspond to
    Nsym=length(s)/ovsamp;
%    Nsym

  elseif ((source_signal_type>=20) && (source_signal_type<=22)) % continuous tone
    % determine Nsym such that signal is 2 second long
    Nsym=2.0/Tsym;
    t_samples=(1:1:Nsym*ovsamp)*Tsamp;

    if     (source_signal_type==20)
      f_tone=440;
      s=sin(2*pi*f_tone .* t_samples); % Hz
    elseif (source_signal_type==21)
      f_tone=1000;
      s=sin(2*pi*f_tone .* t_samples); % Hz
    elseif (source_signal_type==22)
      f_tone=440;
      s=sin(2*pi*f_tone .* t_samples);
      f_tone=1000; % Hz
      s=s+sin(2*pi*f_tone .* t_samples);
    end


    % apply Hamming window to smoothen out edge effects (dc-offset...)   problem... power offset...
    % only at start/end
    s_len=length(s);
    win_len=round(s_len/20);  % 1/20 of entire signal should be used for ramp-up/ramp-down
    if (rem(win_len,2)==1) win_len=win_len+1; end % make even

    for (n=1:win_len/2)
      s(n)=s(n)*(0.54-0.46*cos(2*pi*(n-1)/(win_len-1)));     % ramp-up at start
      s(s_len-(n-1))=s(s_len-(n-1))*(0.54-0.46*cos(2*pi*(n-1)/(win_len-1)));     % ramp-down at end
    end

    s=s./ sqrt(ovsamp); % since re-normalized later again...

  else %%%%%%%%%%%%%%%%%%%%%%%%%%          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);

    if (complex_bb_en==0)
      s=s_re;
    else

      if (source_signal_type==1)  % normal complex signal
        [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);
      % analytical signal, i.e., im is Hilbert-trafo of re (or -im for Lower SB)
      elseif ((source_signal_type==2) || (source_signal_type==3))
        s_im=conv(s_re, h_hilbert_filter);
        if (source_signal_type==3)  s_im=-s_im; end % for lower sideband (LSB)
        s_re=[ zeros(1,hilbert_filter_len*ovsamp/2), s_re, zeros(1,hilbert_filter_len*ovsamp/2-1) ];
      end

      norm=1.0/sqrt(2);  % due to complex constellation
      s = norm*(s_re+j*s_im);
    end
  end



  if (carrier_en==-1)
    dc_offset=abs(min(real(s)));
  elseif (carrier_en==-2)       % carrier only switch
    dc_offset=0;
  else
    dc_offset=carrier_en;
  end

  s = s.* sqrt(ovsamp);  % account for zeroes in between symbols prior to pulse shaping... (see also pulse-shaping demo; adjust...)


% NORMALIZE AGAIN...
%s=normalize_time_signal(s);

  s=s+dc_offset;


  s_tx_bb=s;  % store

  s_len = length(s);

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





  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, modulate

%  tx_fc=tx_fc/Tsym;      % unnormalize, to use multiples of Rs as fc

  if (tx_mod==0)
    s_carrier=cos(   2*pi*tx_fc .* (1:1:s_len)*Tsamp);
  elseif (tx_mod==1)
    s_carrier=-sin(   2*pi*tx_fc .* (1:1:s_len)*Tsamp);
  elseif (tx_mod==2)
    s_carrier=exp(j*2*pi*tx_fc .* (1:1:s_len)*Tsamp);
  end

  s=s .* s_carrier;

  if (tx_mod==2)  s = real(s);  end

  s_tx_pb=s;  % store

  % END, modulate
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



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

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

% only on re?
% also, ovsamp etc... check

  if (SNR_dB<100)
    SNR_lin=10^(SNR_dB/10.0);
    av_power_sqrt=sqrt(compute_av_power(s));
%    sn=av_power_sqrt/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));
    sn=av_power_sqrt/sqrt(SNR_lin) * sqrt(ovsamp);
    n=sn*randn(1,s_len);
    s=s+n;
  end

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


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


  s_rx_pb=s;  % store


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, demodulate

%  rx_fc=rx_fc/Tsym;      % unnormalize
%  rx_phi=pi/2;

  if (rx_demod==0)
    s_carrier=cos(   2*pi*rx_fc .* (1:1:s_len)*Tsamp+rx_phi);
    s=s .* s_carrier;
  elseif (rx_demod==1)
    s_carrier=-sin(   2*pi*rx_fc .* (1:1:s_len)*Tsamp+rx_phi);
    s=s .* s_carrier;
  elseif (rx_demod==2)
    s_carrier=exp(-j*2*pi*rx_fc .* (1:1:s_len)*Tsamp+rx_phi);
    s=s .* s_carrier;
  elseif (rx_demod==3)
    s=s.^2;           % envelope detection: square law detector
  elseif (rx_demod==4)
    for (n=1:s_len)
      if (real(s(n))<0)  s(n)=0;  end   % envelope detection: cut neg. ampl.
    end
  elseif (rx_demod==5)
    s=abs(s);         % envelope detection: abs
  elseif (rx_demod==6)  % same as exp(j...); just for verification
    s_carrier=cos(2*pi*rx_fc .* (1:1:s_len)*Tsamp+rx_phi);
    s_re=s .* s_carrier;
    s_carrier=-sin(2*pi*rx_fc .* (1:1:s_len)*Tsamp+rx_phi);
    s_im=s .* s_carrier;
    s=s_re+j*s_im;
  end

%  s=s .* sqrt(2);

  % END, demodulate
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  s_rx_bb_demod=s;  % store


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, rx filter

  if (rx_filter_bw_en>0)
    s=conv(s, h_rx_filter);
    s_len=length(s);

    % normalize (filter should have 0dB gain in passband)
    %
    s=s./sqrt(ovsamp);

    % 0dB ONLY IN PASSBAND of R_s bandwidth...!
    s=s .* sqrt(rx_filter_t_stretch);

%    tmp_pow_sqrt=sqrt(sum(s.^2));
%    s=s./tmp_pow_sqrt*sqrt(ovsamp*rx_filter_sym_len);  % normalization...? scales with lengths... due to window in fft?
  end

  % END, rx filter
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  s_rx_bb_filt=s;  % store













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


  % select signal to plot
  if     (display_signal==0)  s=s_tx_bb; title_str="TX message (BB)";
  elseif (display_signal==1)  s=s_tx_pb; title_str="TX modulated output (PB)";
  elseif (display_signal==2)  s=s_rx_pb; title_str="RX input (PB)";
  elseif (display_signal==3)  s=s_rx_bb_demod; title_str="RX after demodulation (BB)";
  elseif (display_signal==4)  s=s_rx_bb_filt;  title_str="RX after demod./filter. (BB)";
  else
  end



  clf

  if ((td_en && fd_en) || (td_en && spectrogram_en))
    subplot(2, 1, 1)
  end


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %% BEGIN, plot time domain signal

  if td_en
    hold all;

    grid on;


    s_len = length(s);
    ymin=-2-dc_offset; ymax=2+dc_offset;

    ylim([ymin, ymax]);
    xlabel("time t in sec.");
    ylabel("amplitude x(t)");


    if (tdomain_display_tmin_tmax_en)
      xmin=floor(display_t_min/Tsamp);
      xmax=floor(display_t_max/Tsamp)-1;
      if (xmax>s_len-1) xmax=s_len-1; end
      if (xmin<0) xmin=0; end

      x_axis=(xmin:xmax)*Tsamp;
      xlim([xmin, xmax]*Tsamp);
    else
      % only plot from middle section of t-domain signal to keep "readable"
      xmin=floor(max(0,       round(s_len/2-Nsym_display/2*ovsamp)));
      xmax=floor(min(s_len-1, round(s_len/2+Nsym_display/2*ovsamp)));

      x_axis=(xmin:xmax)*Tsamp;
      xlim([xmin, xmax]*Tsamp);
    end


    xcenter=(xmin+(xmax-xmin)/2)*Tsamp;
    % center display around t=0 (only for display...)
    if ((t_display_ctr_zero_en) && ~(tdomain_display_tmin_tmax_en))
      x_axis = x_axis .- xcenter;
      xlim([xmin, xmax]*Tsamp-xcenter);
      xcenter=0;
    end



    if ((fd_en==0) && (spectrogram_en==0))  title([ title_str, ", t-domain"]);
    elseif (fd_en && (spectrogram_en==0))   title([ title_str, ", t- and f-domain"]);
    elseif ((fd_en==0) && (spectrogram_en)) title([ title_str, ", t-domain, spectrogram"]);
    end




% if ((carrier_en!=0) && (carrier_en!=-2))
%    title(["dc-offset=" num2str(dc_offset)]);
% end


    if (carrier_en==-2) % carrier only
      ylim([-2, 2]);
      plot(x_axis, real(s_carrier(xmin+1:1:xmax+1)), "-", "markersize", 1, 'linewidth',def_lw, "color", "b");
    else
      plot(x_axis, real(s(xmin+1:1:xmax+1)), "-", "markersize", 1, 'linewidth',def_lw, "color", blue(6,:));            % oversampled signal
      if ( (complex_bb_en==1) && ((display_signal==0) || (display_signal>=3)) )
        plot(x_axis, imag(s(xmin+1:1:xmax+1)), "--", "markersize", 1, 'linewidth',def_lw, "color", amber(6,:));            % oversampled signal
      end
    end

    % plot envelope
    if ((tx_fc!=0) && (carrier_en!=-2) && (display_signal>=1) && (display_signal<=3) && (complex_bb_en==0) )
      plot(x_axis, +real(s_tx_bb(xmin+1:1:xmax+1)), "-", "markersize", 1, 'linewidth',def_lw_thick, "color", red(4,:));
      plot(x_axis, -real(s_tx_bb(xmin+1:1:xmax+1)), "-", "markersize", 1, 'linewidth',def_lw_thick, "color", red(7,:));
    end

text(xcenter,0, warning_message_str, 'fontsize', 8,'color', magenta(6,:), 'horizontalalignment', 'center' );

    hold off
  end

  %% END, plot time domain signal
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  if (td_en && fd_en)
    subplot(2, 1, 2)
  end

  if (fd_en && spectrogram_en)
    subplot(2, 1, 1)
  end



  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %% BEGIN, plot frequency domain signal

  if fd_en
    hold all;
    grid on;

    if (carrier_en==-2) s=s_carrier; end

    if (td_en==0)  title([title_str, ", f-domain"]); end

    % keep bins in displayed freq. range constant
    % in BW 2fmax, have N_nyq_fbins bins
    % in BW display_delta_f=xmax-xmin,  N_nyq_fbins_new = N_nyq_fbins * 2fmax/display_delta_f
    if (spectrum_display_fmin_fmax_en)
      display_delta_f = display_f_max-display_f_min;
      N_nyq_fbins_new = max(N_nyq_fbins, round(N_nyq_fbins * (2*fmax)/display_delta_f));   % OVERWRITE WITH NEW VALUE
      if (rem(N_nyq_fbins_new,2)==0) N_nyq_fbins_new=N_nyq_fbins_new+1; end  % make odd
    else
      N_nyq_fbins_new = N_nyq_fbins;
    end


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

    % set up plot output

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


    if (spectrum_display_fmin_fmax_en)  % choose selected portion of spectrum to display
      xmin_tmp=-N_disp_fbins/2*fsamp;
      xmax_tmp= N_disp_fbins/2*fsamp;
      xmin=max(display_f_min, xmin_tmp);
      xmax=min(display_f_max, xmax_tmp);
    elseif (show_neg_freq_en==0)
      xmin=0;
      xmax=N_disp_fbins/2*fsamp;
    else
      xmin=-N_disp_fbins/2*fsamp;
      xmax=+N_disp_fbins/2*fsamp;
    end


    xlim([xmin, xmax]);       % fix to relevant spectral portions


    xcenter=(xmin+(xmax-xmin)/2);


    xlabel("frequency f in Hz");

    if (spectrum_dB_en)
      ylabel("spectrum |X(f)|^2 [dB]");

      ylim([-40, 20]);

      y_plot = 10*log10(H_magn_sq);
    else
      if (H_magn_sq_en)  ylabel("spectrum |X(f)|^2");
                   else  ylabel("spectrum |X(f)|");    end

      ylim([-0.2, 2]);

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

    plot(f_axis*fsamp, y_plot, 'linewidth',def_lw, 'color', green(3,:));


    hold all;


    % also plot theoret. RX-filter curve
    if ((rx_filter_bw_en>0) && (display_signal>=3) && (display_rx_filter_en==1))
      H_ref_magn_sq = get_fdomain_imp_resp_H_raised_cosine(rx_filter_alpha, N_nyq_fbins_new/2, rx_filter_t_stretch) .^ 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', def_lw_thick, "m");
      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', def_lw_thick, "m");
      end
    end

    text(xcenter,0, warning_message_str, 'fontsize', 8,'color', magenta(6,:), 'horizontalalignment', 'center' );

    hold off

  end

  %% END, plot frequency domain signal
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  if ((td_en && spectrogram_en) || (fd_en && spectrogram_en))
    subplot(2, 1, 2)
  end


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %% BEGIN, plot spectrogram
  %% test with   modulation(1, 8e3, 20, 16e3 , 0, 0, 200, 16e3, 0, 1.5, 0.1, 0, 2000,2, 0, 1, 10, 14, 0,1000, 2, 101, 0, 1, 0, 5, 0,0,"dummy")
  %% - print Warning on time/frequency uncertainty... Küpfmüller uncertainty
  %% - what if too many t-steps, such that NFFT too small? (single impulse...)
  %% PROBLEM to fix: neg+pos freq and 440Hz


  if spectrogram_en
    % build t-axis
    % short-term FFT, windowed
    % heatmap plot, colorbar

    hold all;
    grid on;

    if ((td_en==0) && (fd_en==0))  title([title_str, ", spectrogram"]); end


    s_len = length(s);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % BEGIN, t-axis, compare to t-domain signal above
    if (tdomain_display_tmin_tmax_en)
      t_min=floor(display_t_min/Tsamp);
      t_max=floor(display_t_max/Tsamp);
      if (t_max>s_len) t_max=s_len; end
      if (t_min<0) tmin=0; end
      xlim([t_min, t_max]*Tsamp);
      x_axis=(t_min:t_max)*Tsamp;
    else
      % only plot from middle section of t-domain signal to keep "readable"
      t_min=max(0,       round(s_len/2-Nsym_display/2*ovsamp));
      t_max=min(s_len, round(s_len/2+Nsym_display/2*ovsamp));
      xlim([t_min, t_max]*Tsamp);
      x_axis=(t_min:t_max)*Tsamp;
    end

    xcenter=(t_min+(t_max-t_min)/2)*Tsamp;
    % center display around t=0 (only for display...)
    if ((t_display_ctr_zero_en) && ~(tdomain_display_tmin_tmax_en))
      x_axis = x_axis .- xcenter;
      xlim([t_min, t_max]*Tsamp-xcenter);
      xcenter=0;
    end

    % t_min, t_max determined

    % END, t-axis
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    s=s(t_min+1:t_max); % index... t_max without +1!


    % display how many time steps? for now, too many for spectrogram...
    N_t_steps = N_t_steps_spectrogram; %=100;  could be a parameter
    x_len=length(x_axis);
    t_min_tmp=x_axis(1);
    t_max_tmp=x_axis(end);
    t_step=(t_max_tmp-t_min_tmp)/N_t_steps;
    x_axis=t_min_tmp+(0:N_t_steps)*t_step;  % one more than N_t_steps

    x_ax = x_axis;         % time


    % cut input s to respective length


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % BEGIN, f-axis, compare to f-domain signal above

    % keep bins in displayed freq. range constant
    % in BW 2fmax, have N_nyq_fbins bins
    % in BW display_delta_f=xmax-xmin,  N_nyq_fbins_new = N_nyq_fbins * 2fmax/display_delta_f
    if (spectrum_display_fmin_fmax_en)
      display_delta_f = display_f_max-display_f_min;
      N_nyq_fbins_new = round(N_nyq_fbins * (2*fmax)/display_delta_f);   % OVERWRITE WITH NEW VALUE
      if (rem(N_nyq_fbins_new,2)==0) N_nyq_fbins_new=N_nyq_fbins_new+1; end  % make odd
    else
      N_nyq_fbins_new = N_nyq_fbins;
    end

    s_shortFFT_len=floor(length(s)/N_t_steps);
    Nsym_shortFFT = floor(Nsym/N_t_steps);

for (i=1:N_t_steps)
    % pick respective segment of the signal
    s_shortFFT = s( 1+(i-1)*s_shortFFT_len : 1+(i-1)*s_shortFFT_len + s_shortFFT_len-1);

    % BEGIN, WINDOW
%    if (rem(win_len,2)==1) win_len=win_len+1; end % make even
    for (n=1:s_shortFFT_len)
%      s_shortFFT(n)=s_shortFFT(n)*(0.54-0.46*cos(2*pi*(n-1)/(s_shortFFT_len-1))); % Hamming window
      s_shortFFT(n)=s_shortFFT(n)*(0.5-0.5*cos(2*pi*(n-1)/(s_shortFFT_len-1))); % von Hann window (RC-window)
    end
    % END, WINDOW

    [H_magn_sq, f_axis, N_disp_fbins] = compute_spectrum(s_shortFFT, Nsym_shortFFT, ovsamp, spectrum_display_bw_factor, N_nyq_fbins_new);

    if (i==1)
      z_plot=H_magn_sq';
    else
      z_plot=[z_plot, H_magn_sq'];
    end
end
    % set up plot output

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


    if (spectrum_display_fmin_fmax_en)  % choose selected portion of spectrum to display
      f_min_tmp=-N_disp_fbins/2*fsamp;
      f_max_tmp= N_disp_fbins/2*fsamp;
      f_min=max(display_f_min, f_min_tmp);
      f_max=min(display_f_max, f_max_tmp);
    elseif (show_neg_freq_en==0)
      f_min=0;
      f_max=N_disp_fbins/2*fsamp;
    else
      f_min=-N_disp_fbins/2*fsamp;
      f_max=+N_disp_fbins/2*fsamp;
    end
    % now f_min, f_max determined


    % only spectrum withing f_min, f_max needs to be delivered to imagesc
    % cut down z_plot respectively; then ylim[] could be omitted
    f_min_tmp=f_axis(1)*fsamp;
%    f_max_tmp=f_axis(end)*fsamp;

    f_min_idx=floor(abs(f_min-f_min_tmp)/fsamp)+1;
    f_max_idx=floor((f_max-f_min_tmp)/fsamp)+1;

    % cut out respective portion of spectrum to plot
    z_plot=z_plot(f_min_idx:f_max_idx,:);
    f_axis=f_axis(f_min_idx:f_max_idx);


    % END, f-axis, compare to f-domain signal above
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % BEGIN, set up z-axis
    if (spectrum_dB_en)

      zlabel("spectrum |X(f)|^2 [dB]");

      zlim([-40, 20]);

      z_plot = 10*log10(z_plot);
    else
      if (H_magn_sq_en)  zlabel("spectrum |X(f)|^2");
                   else  zlabel("spectrum |X(f)|");    end

      zlim([-0.2, 2]);

      if (H_magn_sq_en)
        z_plot = z_plot;
      else
        z_plot = z_plot .^ (0.5);
      end
    end
    % END, set up z-axis
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    z_plot=[z_plot, z_plot(1:end,end)];  % add identical column at end; since one more entry in x-axis as well

    y_ax = f_axis*fsamp;   % freq

spectrogram_tf_uncertainty_en=1;
if (spectrogram_tf_uncertainty_en==1)
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    % BEGIN, indicate t/f-resolution in plot   delta f * delta t = 1
    t_step_sec=t_step;
    f_step_Hz = 1.0/t_step;
%    f_step_Hz/fsamp
    z_plot_max=max(max(z_plot));
    z_plot_min=min(min(z_plot));
    for (n=1:floor(f_step_Hz/fsamp)+1)
      if (rem(n,2)==1) z_plot(n+1,1)=z_plot_max;
      else z_plot(n+1,1)=z_plot_min;
      end
    end
    % END
    %%%%%%%%%%%%%%%%%%%%%%%%%%
end


    xlabel("time t in sec.");
    ylabel("frequency f in Hz");
    imagesc(x_ax, y_ax, z_plot);
    set(gca,'ydir','normal');


 % strangely enough, on webserver xlim etc needs to be after imagesc...?
 % ok, since Octave 3.6.4 on server, here 4.0.x
 xlim([t_min, t_max]*Tsamp); % repeat...
 ylim([f_min, f_max]);       % fix to relevant spectral portions



   % BEGIN, COLORMAP, use "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 );
   % END, COLORMAP, use "jet light"

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

    hold off
  end

  %% END, plot spectrogram
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



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




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

  if (webdemo_en==0)
    filename=[filename,num2str(display_signal)];
  end

  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 [y,ampl_vec]=get_random_PAM_vector(nb_bits, unipolar_en, N)
  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
  P_av = 1/M*sum(ampl_vec.^2);            % NORMALIZE to av. power 1 at discrete time samples
  ampl_vec=ampl_vec/sqrt(P_av);
  y=discrete_rnd(ampl_vec, prob_vec, 1, N);

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


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

  if (mod(len,2)==0)
%    t_n = (-floor(len/2):1/ovsamp:floor(len/2)-1/ovsamp);
    t_n = (-floor(len/2)+1/ovsamp*0.5:1/ovsamp:floor(len/2)-1/ovsamp+1/ovsamp*0.5);   % avoid time 0 for filter (inf.)
  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==0)
      h(k)=-10;   % should never happen, see above
    else
      h(k) = 1/(pi.*n);
    end
    k=k+1;
  end

  % normalize
  av_pow=1; % just to initialize
  av_pow = sum(h .^ 2);

  if (normalize_en==1)       % normalize such that 0dB gain overall gain (t-domain)
    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, 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=1; % just to initialize
  av_pow = sum(h .^ 2);

  if (normalize_en==1)       % normalize such that 0dB gain overall gain (t-domain)
    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, t_stretch)

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

  % 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)

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

  % 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] = compute_spectrum(s, Nsym, ovsamp, spectrum_display_bw_factor, N_nyq_fbins)

  s_len = length(s);

% apply Hamming window to smoothen out edge effects (dc-offset...)   problem... power off
%  for (n=1:s_len)
%    s(n)=s(n)*(0.54-0.46*cos(2*pi*(n-1)/(s_len-1)));
%  end


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

% N_nyq_fbins=101; fix

  % 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)

  % NEW in ampl. demo; to normalize to 1...
  y=y./sqrt(ovsamp); % nicer to read...
  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 s_return=normalize_time_signal(s)
  s_len=length(s);
  pow_sqrt=sqrt(sum(s.^2)/s_len);
  s_return=s./pow_sqrt;
end

function av_power_return=compute_av_power(s)
  s_len=length(s);
  av_power_return=sum(s.^2)/s_len;
end



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