% WEBDEMO .m-file, Institute of Telecommunications, University of Stuttgart, Germany
% Author: Stephan ten Brink, 27.11.2013
% EXAMPLE:  rayleigh_fading(2000, 1, 0.01, 0,0,0, -12345, 8, 1, 1, 10,  1, [0 1 0 0 0 0], 0, 0, "rayleigh_fading")
% PARAMETERS:
%   nb_samples,          % nb of t-domain samples (Ts=1) to generate
%   Ts,                  % duration of one discrete-time sample in sec. (set to 1 for normalization)
%   fD,                  % max. Doppler freq. shift in Hz       (note: only Ts*fD relevant)
%   speed_v_kmh,         % alternatively, Doppler shift can be defined through vehicle speed;  in km/h
%   fc,                  % carrier frequency in Hz
%   speed_en,            % 1=use speed input for computing Doppler frequency, not fD
%   t_offset,            % time offset to improve randomness; if -1, then random offset
%   N0,                  % nb of scatterers used to generate complex Rayleigh fading phasor
%   nb_ray_out,          % set to 1 for now; nb of output time series
%   h_magn_sq_dB_en,     % plot t-domain signal power in dB
%   ac_nb_periods,       % how many periods to display in auto-correlation plot
%   spectrum_dB_en,      % plot spectrum y-axis in dB
%   plot_en,             % enable plots 0,1 [t-domain, t-domain magn_sq, t-dom hist, t-dom hist magn_sq, auto-corr, f-dom spectrum]
%   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
% rayleigh_fading(2000, 1e-6, 0, 100,1e9,1, -12345, 8, 1, 1, 10,  1, [0 1 0 0 0 0], 0, 0, "rayleigh_fading")
% DAB, Ts=1.246us, v=100km/h, fc=200MHz: rayleigh_fading(2000, 1.246e-3, 0, 100, 200e6,1, -12345, 8, 1, 1, 10,  1, [0 1 0 0 0 0], 0, 0, "rayleigh_fading")


% a central building block for various indoor/outdoor channel models
% totally deterministic (except for random offset)

% possible extensions:
% - add: spectrum
% [- Walsh/Hadamard to generate multiple independent time series at once]


% BEGIN, main function (control)

function rayleigh_fading(
  nb_samples,          % nb of t-domain samples (Ts=1) to generate
  Ts,                  % duration of one discrete-time sample in sec. (set to 1 for normalization)
  fD,                  % max. Doppler freq. shift in Hz       (note: only Ts*fD relevant)
  speed_v_kmh,         % alternatively, Doppler shift can be defined through vehicle speed;  in km/h
  fc,                  % carrier frequency in Hz
  speed_en,            % 1=use speed input for computing Doppler frequency, not fD
  t_offset,            % time offset to improve randomness; if -1, then random offset
  N0,                   % nb of scatterers used to generate complex Rayleigh fading phasor
  nb_ray_out,          % set to 1 for now; nb of output time series
  h_magn_sq_dB_en,     % plot t-domain signal power in dB
  ac_nb_periods,       % how many periods to display in auto-correlation plot
  spectrum_dB_en,      % plot spectrum y-axis in dB
  plot_en,             % enable plots 0,1 [t-domain, t-domain magn_sq, t-dom hist, t-dom hist magn_sq, auto-corr, f-dom spectrum]
  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);
  set(0, "defaulttextfontname", "Helvetica");
  set(0, "defaulttextfontsize", 8);

  if (webdemo_en)
    set(0, "defaultaxesfontsize", 8);
    set(0, "defaultaxesfontname", "Helvetica");
  elseif (svg_en)
    set(0, "defaultaxesfontsize", 10);
    set(0, "defaultaxesfontname", "Helvetica");
  else
    set(0, "defaultaxesfontsize", 10);
    set(0, "defaultaxesfontname", "Helvetica");
  end


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


  if (speed_en)
    fD = fc*(speed_v_kmh*1000/3600.0)/3e8;
  else % take fD from commandline
    speed_v_kmh = fD/fc/1000*3600*3e8;
  end

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



  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % 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];
    warning("setting  'plot_en = [1 0 0 0 0 0]'");
  end

  if (nb_ray_out!=1)
    nb_ray_out=1;
    warning("setting 'nb_ray_out=1'");
  end

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

  if ((N0<1) || (N0>32))
    N0=8;            % reset to default
    warning("setting 'N0=8'");
  end

  if ((TsfD<0.00001) || (TsfD>1))
    TsfD=0.01;       % reset to default
    Ts=1; fD=TsfD;
    warning("setting 'TsfD=0.01'");
  end
  
  if ((ac_nb_periods<1) || (ac_nb_periods>50))
    ac_nb_periods=10;       % reset to default
    warning("setting 'ac_nb_periods=10'");
  end

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


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

  h_ray = get_time_sequence(TsfD, nb_samples, N0, t_offset);



  clf % other resets/clear commands needed?
  hold on

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

  if (plot_en(1))
    t = (0:1:nb_samples-1)*Ts;
    xlabel("time in sec.");
    ylabel("re, im-part of complex phasor");
    plot(t, real(h_ray), "linewidth", 2, "color", "r");
    plot(t, imag(h_ray), "linewidth", 2, "color", "b");
    output_format = landscape;
  elseif (plot_en(2))
    t = (0:nb_samples-1)*Ts;
    xlabel("time in sec.");
    magn_sq = real(h_ray .* conj(h_ray));
    if (h_magn_sq_dB_en)
      ylim([-40, 10]);
      magn_sq_dB = 10*log10(magn_sq);
      ylabel("squared magnitude of phasor in dB");
      plot(t, magn_sq_dB, "linewidth", 2, "color", "b");
    else
      ylabel("squared magnitude of phasor");
      plot(t, magn_sq, "linewidth", 2, "color", "b");
    end
    output_format = landscape;
  elseif (plot_en(3))
    xlim([-3, 3]);
    xlabel("real and imag. part of complex baseband phasor");
    ylabel("distribution");
    hist(real(h_ray), 50, 1, "facecolor", "r", "edgecolor", "k");
    hist(imag(h_ray), 50, 1, "facecolor", "b", "edgecolor", "k");
  elseif (plot_en(4))
    xlabel("magnitude, and squared magnitude (power) of time-domain signal");
    ylabel("distribution");
    magn = sqrt(real(h_ray .* conj(h_ray)));
    magn_sq = real(h_ray .* conj(h_ray));
    x_step = 0.1;
    x = (0:x_step:6)+x_step/2;
    text(2.5, 0.5, "squared magnitude", "fontsize", 12, "color", "b");
    text(2.5, 0.8, "magnitude", "fontsize", 12, "color", "r");

    hist(magn_sq, x, 1*1/x_step, "facecolor", "b", "edgecolor", "k");
    hist(magn, x, 1*1/x_step, "facecolor", "r", "edgecolor", "k");

    % theoret. pdf for magn_sq: exp
    p = 1;
    pdf_magn_sq = 1/p*exp(-x/p);
    plot(x, pdf_magn_sq, "linewidth", 2, "color", "m");


    % theoret. pdf for magn: Rayleigh
    sigma_sq_R = 0.5; % note: E[rayleigh_ampl^2] = 2*sigma_sq_R;
    pdf_magn = x/sigma_sq_R .* exp(-x .^ 2/(2*sigma_sq_R));
    plot(x, pdf_magn, "linewidth", 2, "color", "m");
  elseif (plot_en(5))
    nb_samples_display = min(nb_samples, floor(ac_nb_periods / TsfD));
    xlim([0, nb_samples_display*Ts]);
    ylim([-0.5, 1.2]);
    xlabel("time lag in sec.");
    ylabel("auto-correlation");
    t = (0:nb_samples_display-1)*Ts;
    ac = compute_auto_correlation(h_ray, nb_samples_display);
    plot(t, ac, "linewidth", 2, "color", "b");
    plot(t, besselj(0,2*pi*fD*t), "linewidth", 2, "color", "m");
    legend("auto-corr. of t-domain signal","reference auto-corr. (Bessel)");
    output_format = landscape;
  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, "-dpng",output_format);
%print(filename, "-djpg","-S1280,960", "-r0");
  if (svg_en)
    print(filename, "-dsvg",output_format);
  end


end

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






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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y_vec = compute_auto_correlation(x, nb_samples_display)
  len = length(x);
  norm = len;

  t_lag_max = nb_samples_display;

  y_vec = zeros(1,t_lag_max);
  for t_lag = (0:t_lag_max-1)
    x1 = [ x, zeros(1,t_lag) ];
    x2 = [ zeros(1,t_lag), x ];

    y_vec(1+t_lag)=sum(x1 .* conj(x2))/norm;
  end

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




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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_samples, N0, t_offset)

  % some constants
  wm = 2*pi*TsfD;

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

  t = t_offset;

  for n=1:nb_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_samples;  % checked 2013/11/27, E[is^2]=0.5, ok
%  mean_iq = sum(qs_vec .* qs_vec)/nb_samples;  % checked 2013/11/27, E[qs^2]=0.5, ok

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



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