% WEBDEMO .m-file, Institute of Telecommunications, University of Stuttgart, Germany
% Author: Stephan ten Brink, 09.12.2013
% EXAMPLE: tf_sel_fading(60, 128, 0.1, 1e6, 6, 1, 100, 200e6, -1, 1, 0, 0, 0, "tf_sel_fading")
% PARAMETERS:
%   nb_t_samples,        % nb of t-domain samples
%   nb_f_samples,        % nb of f-domain samples
%   time_interval,       % time interval to plot, in sec
%   bandwidth,           % bandwidth to plot, in Hz
%   nb_taps,             % tapped delay line model, nb of taps
%   delay_stretch,       % factor, unitless; max. delay is defined as (1/bandwidth) * delay_stretch
%   speed_v_kmh,         % Doppler shift defined through vehicle speed;  in km/h
%   fc,                  % carrier frequency in Hz
%   t_offset,            % time offset to improve randomness; if -1, then random offset
%   plot_delay_profile_en, % plot delay profile as inset plot
%   plot_idx_en,         % 0=heatmap, 1=surface, 2=surface
%   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
% DAB, Ts=1.246us, v=100km/h, fc=200MHz: tf_sel_fading(60, 128, 0.1, 1e6, 6, 1, 100, 200e6, -1, 0, 0, 0, 0, "tf_sel_fading")
% WLAN 802.11a: Ts=4us, v=0.01km/h, fc=2.45e9Hz, bw=20MHz: tf_sel_fading(60, 128, 0.1, 20e6, 6, 0.6, 0.001, 2.45e9, -1, 0, 1, 0, 0, "tf_sel_fading")



% NOTE
% - only f-selective (fading coeff set to const.; e.g. by very slow speed, v=0.0001)
% - only t-selective (by having no delay spread, or simply setting nb_taps=1, or delay_stretch to close to 0)
% possible add-ons:
% - plot in fig.: TD, Bcoh for f-sel; and Tcoh, Ds for t-sel;
% - plot suggested pilot grid
% - different standard scenarios: DAB; WLAN; LTE 10MHz
% - delay profile: choose delay profile; uniform, exp., random, or fixed



% BEGIN, main function (control)

function tf_sel_fading(
  nb_t_samples,        % nb of t-domain samples
  nb_f_samples,        % nb of f-domain samples
  time_interval,       % time interval to plot, in sec
  bandwidth,           % bandwidth to plot, in Hz
  nb_taps,             % tapped delay line model, nb of taps
  delay_stretch,       % factor, unitless; max. delay is defined as (1/bandwidth) * delay_stretch
  speed_v_kmh,         % Doppler shift defined through vehicle speed;  in km/h
  fc,                  % carrier frequency in Hz
  t_offset,            % time offset to improve randomness; if -1, then random offset
  plot_delay_profile_en, % plot delay profile as inset plot
  plot_idx_en,         % 0=heatmap, 1=surface, 2=surface
  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
)

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % 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 = portrait;     % default


  set(0, "defaultlinelinewidth", 2);

  if (webdemo_en)
%    set(0, "defaultaxesfontsize", 8);
    set(0, "defaultaxesfontsize", 7);
    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, "defaultaxesfontsize", 5);
    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;
  endif

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


  N0 = 16;  % nb of oscillators used to generate complex Rayleigh fading phasor

  Ts = time_interval/nb_t_samples;   % duration of one time-sample
  delta_f = bandwidth/nb_f_samples;  % bandwidth of one sample in freq.-direction

  % max Doppler-shift
  fD = fc*(speed_v_kmh*1000/3600.0)/3e8;

  % note: only Ts*fD relevant; Ts, fD separate only for labeling of plot axes
  TsfD = Ts*fD;



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

  if (t_offset<0)
    t_offset = floor(rand(1)*1e6);   % random time offset
  end

  if (nb_t_samples*nb_f_samples>64000)
    nb_t_samples = 800;
    nb_f_samples = 800;
    warning("setting  '(nb_t_samples*nb_f_samples>64000)'");
  end

  if ((nb_taps<0) || (nb_taps>6))
    nb_taps = 6;
  end

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


  % RGB color definitions, http://www.rapidtables.com/web/color/RGB_Color.htm
  gray =    [[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;


  % for simplicity, define fixed delay profile
  tap_weights = [1, 0.8, 0.2, 0.7, 0.4, 0.1];
  tap_weights = tap_weights(1:nb_taps);

  tap_power = sum(tap_weights.*tap_weights);
  tap_weights = tap_weights ./ sqrt(tap_power);  % normalize tap weights such that overall power it 1, i.e., gain 0dB

  tap_delays = [0, 1, 4, 8, 10, 14]*(1/bandwidth)*delay_stretch;    % normalized to Ny-BW (relative, for simplicity)
  tap_delays = tap_delays(1:nb_taps);

  % generate time-variant time-domain impulse response h(t,tau)
  for tap_idx=1:nb_taps
    h_ray(tap_idx,:) = tap_weights(tap_idx) .* get_time_sequence(TsfD, nb_t_samples, N0, t_offset);
    t_offset = t_offset*tap_idx*2376423;    % pseudo-random offset to have uncorrelated Rayleigh fading channel-taps
  end


  % generate time-variant frequency-domain transfer function H(t,f)
  for t_discrete=1:nb_t_samples
    [H_f, f_axis] = get_frequency_response(t_discrete, h_ray, tap_delays, nb_taps, delta_f, nb_f_samples);
    H_t_f_dB(t_discrete,:) = 10*log10(H_f .* conj(H_f));
  end



  t_axis = (1:nb_t_samples)*Ts;

  nb_points =  length(t_axis)*length(f_axis)


  clf reset
  hold on

  title(["Speed v=" num2str(speed_v_kmh) "km/h; max. Doppler f_D=" num2str(fD) "Hz"]);

  colormap("default");

  if (plot_idx_en==0)
    % probably ok wrt run-time, for up to 60000 points
    if (nb_points < 60000)
      xlabel("time in sec.");
      ylabel("frequency in Hz");
      imagesc(t_axis, f_axis, H_t_f_dB', [-30,10]);
      set(gca,'ydir','normal');
      colorbar("EastOutside");
      title(colorbar,"|H(t,f)|^2 in dB");
      inset_pos = [0.5 0.6 0.2 0.2];
    end
    output_format = portrait;
  elseif (plot_idx_en==1)
    if (nb_points < 10000)
      ylabel("time in sec.");
      xlabel("frequency in Hz");
      zlabel("|H(t,f)|^2 in dB");
      meshz(f_axis, t_axis, H_t_f_dB); %, [-30,10]);
      view(135,45);
      set(gca,'xdir','reverse');
      inset_pos = [0.70 0.6 0.2 0.2];
    end
    output_format = portrait;
  elseif (plot_idx_en==2)
    if (nb_points < 10000)
%      surfc(f_axis, t_axis, H_t_f_dB);
      ylabel("time in sec.");
      xlabel("frequency in Hz");
      surf(f_axis, t_axis, H_t_f_dB);
      view(135,45);
      set(gca,'xdir','reverse');
    end
    output_format = portrait;
  end

%      main_ax = gca;

  % inset plot
  if (plot_delay_profile_en)
    sub_ax = axes('position', inset_pos);
    set(sub_ax, 'color', [1 1 1]);
    stem(tap_delays, tap_weights, 'color', 'b');
    set(sub_ax, 'XTick', tap_delays);
%%          set(sub_ax, 'XLim', [-1.1*A 1.1*A]);
%          set(sub_ax, 'XLabel', 'delay');
%          set(sub_ax, 'YLabel', 'path weights');

    % how can this be done better? see next 6 lines...
    tmp_str = {'0'; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''; ''};
%          tmp_str = cell(1,nb_taps);
%          tmp_str(1) = '0';
%          for tmp=2:nb_taps-1     tmp_str(tmp)=' ';   end
    tmp_str(nb_taps) =  ['tau_{max}=',num2str(tap_delays(end)),'s'];
    tmp_str = tmp_str(1:nb_taps);

    set(sub_ax, 'XTickLabel', tmp_str);
%%          set(sub_ax, 'YLim', [-1.1 1.1]);
    set(sub_ax, 'YTick', [0 0.2 0.4 0.6 0.8 1]);
%%          set(sub_ax, 'XGrid', 'on');
%%          set(sub_ax, 'YGrid', 'on');
    set(sub_ax, 'Title', 'delay profile');

%          rect_pos = [0.65 0.45 0.24 0.3];
%          axes(main_ax);
%          rectangle('Position', [0.05, 600000, 0.01, 100000], 'FaceColor',  gray(7,:), 'EdgeColor', 'none');
%          rectangle('Position', fig2ds_pos(main_ax, rect_pos), 'FaceColor',  gray(7,:), 'EdgeColor', 'none');
  end
  
  
  
  
  
  hold off


% Create an image file. This image is displayed on the webpage. Its size
% must be 640x480.
%print(filename, "-dgif","-S640,480", "-r0");
%print(filename, "-djpg","-S1280,960", "-r0");
  if (webdemo_en)
%    print(filename, "-djpg",output_format);
    print(filename, "-dpng",output_format);
  end
  if (svg_en)
    print(filename, "-dsvg",output_format);
  end


end

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






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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% only for ONE generator
function [is, qs] = get_coeff(N0, wm, t)

  is = 0;
  qs = 0;

  for n=1:N0
    beta_n = pi*n/N0;                   % modified Jakes model ("Dent model")
    tmp = cos(wm*t*cos(pi*(n-0.5)/2/N0 ));   % note:  wm = 2*pi*TsfD

    is = is + sin(beta_n)*tmp;
    qs = qs + cos(beta_n)*tmp;
  end

  is = is*sqrt(2/N0);
  qs = qs*sqrt(2/N0);

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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = get_time_sequence(TsfD, nb_t_samples, N0, t_offset)

  % some constants
  wm = 2*pi*TsfD;

  is_vec = zeros(1,nb_t_samples);
  qs_vec = zeros(1,nb_t_samples);

  t = t_offset;

  for n=1:nb_t_samples
    [is, qs] = get_coeff(N0, wm, t);

    is_vec(n) = is;
    qs_vec(n) = qs;

    t=t+1;    % discrete time index
  end

  % check power normalization
%  mean_is = sum(is_vec .* is_vec)/nb_t_samples;  % checked 2013/11/27, E[is^2]=0.5, ok
%  mean_iq = sum(qs_vec .* qs_vec)/nb_t_samples;  % checked 2013/11/27, E[qs^2]=0.5, ok

  y = is_vec + j*qs_vec;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [y, f_axis] = get_frequency_response(t_discrete, h_ray, tap_delays, nb_taps, delta_f, nb_f_samples)

  y=zeros(1,nb_f_samples);

  for f_discrete=1:nb_f_samples
    sum_H = 0;
    f = f_discrete*delta_f;
    for tap_idx = 1:nb_taps
      sum_H = sum_H + h_ray(tap_idx, t_discrete)*exp(-j*2*pi*f*tap_delays(tap_idx));
    end
    y(f_discrete) = sum_H;
  end

  f_axis = (1:nb_f_samples)*delta_f;

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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NOTE: taken from Mach-Zehnder Mod. webdemo
function ds_pos = fig2ds_pos(ah, pos)
% FIG2DS_POS Convert a position from normalized figure to data space units.
%
% IN:
%   ah  Handle of axes that contains the dataspace
%   pos Four element position vector [left bottom width height] in normalized
%       figure units. (Values between 0 and 1.0).
% OUT:
%  ds_pos Four element position vector in data space units of the given axes.

    % Position that is to be converted
    x = pos(1);
    y = pos(2);
    width = pos(3);
    height = pos(4);
    
    % Limits of data space of given axes
    ax_xlim = get(ah, 'XLim');
    ax_ylim = get(ah, 'YLim');
    dsx_max = ax_xlim(2);
    dsx_min = ax_xlim(1);
    dsy_max = ax_ylim(2);
    dsy_min = ax_ylim(1);
    
    % Position of given axes inside the figure
    ax_pos = get(ah, 'Position');
    fig_x = ax_pos(1);
    fig_y = ax_pos(2);
    fig_width = ax_pos(3);
    fig_height = ax_pos(4);
    
    % Conversion
    ds_pos(1) = (dsx_max - dsx_min)/fig_width * (x - fig_x) + dsx_min;
    ds_pos(2) = (dsy_max - dsy_min)/fig_height * (y - fig_y) + dsy_min;
    ds_pos(3) = (dsx_max - dsx_min)/fig_width * width;
    ds_pos(4) = (dsy_max - dsy_min)/fig_height * height;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



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