% WEBDEMO chromatic_dispersion.m-file, Institute of Telecommunications, University of Stuttgart, Germany
% www.inue.uni-stuttgart.de
% author: Michael Bernhard, michibernhard@gmx.de
% date: 17.10.2013
%
% Written for GNU Octave 
%
% EXAMPLE:  chromatic_dispersion(10, 5, 17, 1, 1, 0, 0, "chromatic_dispersion")
% PARAMETERS:
%  fiber_length,    % in km
%  symbol_rate,     % in GHz  
%  chromatic_dispersion, % in ps/km/nm
%  impulse_shaper,  % 1: rectangular, 4: raised cosine  
%  //random_data,     % 1: create random data
%  svg_en,          % if 1, then also output svg-file
%  webdemo_en,      % if 0, then test-mode in Octave/Matlab commandline, no need for webserver running
%  filename         % output file for webdemo





% BEGIN, main function (control)

function chromatic_dispersion(
  fiber_length,     % in km
  symbol_rate,      % in GHz  
  chromatic_dispersion, % in ps/km/nm
  impulse_shaper,   % 1: rectangular, 4: raised cosine  
  svg_en,           % if 1, then also output svg-file
  webdemo_en,       % if 0, then test-mode in Octave/Matlab commandline, no need for 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");
  endif


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

  output_format = landscape; %portrait;     % default


  set(0, "defaultlinelinewidth", 2);

  if (webdemo_en)
    set(0, "defaultaxesfontsize", 8);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
  elseif (svg_en)
    set(0, "defaultaxesfontsize", 6);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 6);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 1;
  else
    set(0, "defaultaxesfontsize", 10);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
  endif


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



  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN check input parameters for consistency
  fiber_length = fiber_length * 1e3; % from km to m
  symbol_rate = symbol_rate * 1e9; % from GHz to Hz
  chromatic_dispersion = chromatic_dispersion * 1e-6; % from ps/km/nm in s/m/m

 
  % END check input parameters for consistency
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  hold off;
  t1 = 1/symbol_rate;
 % if random_data == 1
 %   M = 20;
 %   tx_data = zeros(1,M);
 %   for k=1:M
 %     tx_data(1,k) = round(rand);
 %   endfor
 % else
    tx_data = [0 1 0 1 1 0];
 % endif

  oversampling_factor = 64;

  nb_points = oversampling_factor * size(tx_data,2); % total number of points for plot
  [tx_signal_time tx_signal_data] = create_tx_signal(tx_data, impulse_shaper, t1, nb_points);
 
  % plot intensity of tx signal, time axis in ns
  plot(tx_signal_time*1E9, abs(tx_signal_data).^2, "linewidth",  const_lw, "color", "b");
  hold on;

  % compute rx signal
  [rx_signal_time rx_signal_data] = compute_rx_signal(tx_signal_time, tx_signal_data, chromatic_dispersion, fiber_length, symbol_rate, oversampling_factor);

  plot(rx_signal_time*1e9, abs(rx_signal_data).^2, "linewidth", const_lw, "color", "m");
  xlim([0 rx_signal_time(1,end)*1e9]);

  xlabel('time t/ns');
  ylabel('intensity');
  legend('tx signal', 'rx signal');
  hold off

  % Create an image file. This image is displayed on the webpage. Its size
  % must be 640x480.
  print(filename, "-dpng",output_format);
  
  if (svg_en)
    print(["fig_",filename], "-dsvg",output_format);
  endif


end

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






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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create TX signal
function [tx_signal_time tx_signal_data] = create_tx_signal(tx_data, impulse_shaper, t1, nb_points);
  N=size(tx_data,2);
  tx_signal_data = zeros(1,nb_points);
  tx_signal_time = linspace(0,N*t1,nb_points);
  tx_symbol = zeros(1,nb_points/N);

  for k=1:N % symbol index
    for time=0:nb_points/N-1 % time index
      switch (impulse_shaper)
	case 1 % rectangular
	  tx_symbol(1,time+1) = tx_data(1,k);
	case 2 % rectangular with filter (see some lines below)
	  tx_symbol(1,time+1) = tx_data(1,k);
	case 3 % rectangular with filter (see some lines below)
	  tx_symbol(1,time+1) = tx_data(1,k);
	case 4
	  t_n = (time - nb_points/N/2)/(nb_points/N/2);
	  alpha = 0.5;
	  if (abs(t_n) <= (1-alpha)/2)
	    tx_symbol(1, time+1) = 1 * tx_data(1,k);
	  else
	    if (((1-alpha)/2) < abs(t_n) && abs(t_n) <= (1+alpha)/2)
	      tx_symbol(1,time+1) = (cos(pi/2/alpha * (abs(t_n)-(1-alpha)/2) ))^2 * tx_data(1,k);
	    else
	      tx_symbol(1,time+1) = 0;
	    endif
	  endif;
	case 5
	  t_n = (time - nb_points/N/2)/(nb_points/N/2);
	  alpha = 1;
	  if (abs(t_n) <= (1-alpha)/2)
	    tx_symbol(1, time+1) = 1 * tx_data(1,k);
	  else
	    if (((1-alpha)/2) < abs(t_n) && abs(t_n) <= (1+alpha)/2)
	      tx_symbol(1,time+1) = (cos(pi/2/alpha * (abs(t_n)-(1-alpha)/2) ))^2 * tx_data(1,k);
	    else
	      tx_symbol(1,time+1) = 0;
	    endif
	  endif;
	otherwise       
      endswitch
    endfor
    tx_signal_data(1,1+(k-1)*nb_points/N:k*nb_points/N) = tx_symbol(1,:);
  endfor
  
  if impulse_shaper == 2
    freq = fft(tx_signal_data);
    M = nb_points/3;
    freq(1,nb_points/2-M:nb_points/2+M) = 0;
    tx_signal_data = ifft(freq);
  endif
  if impulse_shaper == 3
    freq = fft(tx_signal_data);
    M = nb_points/2-10;
    freq(1,nb_points/2-M:nb_points/2+M) = 0;
    tx_signal_data = ifft(freq);
  endif

endfunction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute RX signal
function [rx_signal_time rx_signal_data] = compute_rx_signal(tx_signal_time, tx_signal_data, d_c, fiber_length, symbol_rate, oversampling_factor);
  N=size(tx_signal_data,2);
  rx_signal_data = zeros(1, N);
  rx_signal_time = tx_signal_time;

  c = 299792458; % speed of light in m/s
  lm0 = 1550e-9; % reference wavelength in m
  D = d_c;%17e-6; % chromatic dispersion coefficient in s/m^2 at lm0
  L = fiber_length; %100e3; % fiber length in m
  B = oversampling_factor * symbol_rate; % bandwidth
  f = linspace(-B/2, B/2, N);
  
  H = exp( 1j * pi * D * lm0^2 / c * L * f.^2);
  H = my_fftshift(H);

  freq = fft(tx_signal_data);
  freq = freq .* H;
  rx_signal_data = ifft(freq);

endfunction



function retval = my_fftshift (x, dim)
%http://lists.gnu.org/archive/html/octave-bug-tracker/2011-07/msg00014.html

  retval = 0;

  if (nargin != 1 && nargin != 2)
    print_usage ();
  endif

  if (nargin == 2)
    if (!isscalar (dim))
      error ("fftshift: dimension must be an integer scalar");
    endif
    nd = ndims (x);
    sz = size (x);
    sz2 = ceil (sz(dim) / 2);
    idx = cell ();
    for i = 1:nd
      idx{i} = 1:sz(i);
    endfor
    idx{dim} = [sz2+1:sz(dim), 1:sz2];
    retval = x(idx{:});
  else
    if (isvector (x))
      # Error in Version 3.4.0: x must not be overwritten!
      # x = length (x);  
      # xx = ceil (x/2);
      # retval = x([xx+1:x, 1:xx]);
      # Fix: use xl for the length of x
      xl = length (x);  
      xx = ceil (xl/2);
      retval = x([xx+1:xl, 1:xx]);
    elseif (ismatrix (x))
      nd = ndims (x);
      sz = size (x);
      sz2 = ceil (sz ./ 2);
      idx = cell ();
      for i = 1:nd
        idx{i} = [sz2(i)+1:sz(i), 1:sz2(i)];
      endfor
      retval = x(idx{:});
    else
      error ("fftshift: expecting vector or matrix argument");
    endif
  endif

endfunction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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


