% WEBDEMO chromatic_dispersion_eye.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:  
%
% 2-AM, RZ, rectangular, 5km:
% chromatic_dispersion_eye(5, 10, 17, 6, 'rz_2', 0, 0, 0, 0, "chromatic_dispersion_eye")
%
% 2-AM, RZ, raised-cosine (alpha=1), 30km:
% chromatic_dispersion_eye(30, 10, 17, 5, 'rz_2', 0, 0, 0, 0, "chromatic_dispersion_eye")
%
% 4-AM, NRZ, raises-cosine (alpha=1), 30km:
% chromatic_dispersion_eye(30, 10, 17, 5, 'nrz_4', 0, 0, 0, 0, "chromatic_dispersion_eye")
%
% color grading, 2-AM, NRZ, rectangular, 30km:
% chromatic_dispersion_eye(30, 10, 17, 6, 'nrz_2', 0, 1, 0, 0, "chromatic_dispersion_eye")
%
% duobinary, color grading, 120 km
% chromatic_dispersion_eye(120, 10, 17, 5, 'duobinary', 0, 1, 0, 0, "chromatic_dispersion_eye")
%
% PARAMETERS:
%  fiber_length,    % in km
%  symbol_rate,     % in GHz  
%  chromatic_dispersion, % in ps/km/nm
%  impulse_shaper,  % 6: rectangular; 4: raised cosine, alpha = 0.5; 5: raised cosine, alpha = 1  
%  pulse_mode,      % nrz_2, nrz_4, rz_2, rz_4, duobinary, nrz_2_dpsk, rz_2_dpsk
%  intensity_plot,  % 0: plot sqrt of intensity; 1: plot intensity
%  color_grading,   % 0: color grading off, 1: color grading on
%  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_eye(
  fiber_length,     % in km
  symbol_rate,      % in GHz  
  chromatic_dispersion, % in ps/km/nm
  impulse_shaper,   % 6: rectangular;  4: raised cosine, alpha = 0.5; 5: raised cosine, alpha = 1  
  pulse_mode,       % nrz_2, nrz_4, rz_2, rz_4, duobinary, nrz_2_dpsk, rz_2_psk
  intensity_plot,   % 0: plot sqrt of intensity; 1: plot intensity
  color_grading,    % 0: color grading off, 1: color grading on
  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

  hold off;
  % 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;
  M = 100;
  if color_grading == 1
    M = 500;
  end
  tx_data = zeros(1,M);
  for k=1:M
    if strcmp(pulse_mode, 'nrz_4') || strcmp(pulse_mode, 'rz_4') % 4-AM
      tx_data(1,k) = round(rand*3) / 3;
    else    
      tx_data(1,k) = round(rand); % 2-AM, binary
    end
  endfor

  oversampling_factor = 64;
  if impulse_shaper == 6
    oversampling_factor = 256;
  end
 
  nb_points = oversampling_factor * size(tx_data,2); % total number of points for plot
 
  if strcmp(pulse_mode, 'duobinary')
    N = size(tx_data,2);
    tx_data *= 2;
    tx_data -= 1;

    [tx_signal_time tx_signal_data] = create_tx_signal(tx_data, impulse_shaper, pulse_mode, t1, nb_points);

    tx_data_2 = zeros(1,N);
    tx_data_2(1,2:N) = tx_data(1,1:N-1);  
    [tx_signal_time tx_signal_data_2] = create_tx_signal(tx_data_2, impulse_shaper, pulse_mode, t1, nb_points);
    tx_signal_data += tx_signal_data_2;
    tx_signal_data /= 2;
  else
    if strcmp(pulse_mode, 'nrz_2_dpsk') || strcmp(pulse_mode, 'rz_2_dpsk')
      N = size(tx_data,2);
      tx_data *= 2;
      tx_data -= 1;

      [tx_signal_time tx_signal_data] = create_tx_signal(tx_data, impulse_shaper, pulse_mode, t1, nb_points);
      tx_signal_data = tx_signal_data / sqrt(2);
  else
      [tx_signal_time tx_signal_data] = create_tx_signal(tx_data, impulse_shaper, pulse_mode, t1, nb_points);
  end
end
 

 

% 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);
  % optical - electrical conversion
  rx_electrical_signal_data = [];
  if strcmp(pulse_mode, 'nrz_2_dpsk')  || strcmp(pulse_mode, 'rz_2_dpsk')  % using two photodiodes
    %UET2-Skript, 20140108scriptUET2_kommentiert.pdf, page 95
    N = size(rx_signal_data,2);    
    rx_electrical_signal_data_2 = zeros(1,N);
    T = oversampling_factor;
    a = rx_signal_data;
  
    a1 = zeros(1,N);
    a2 = zeros(1,N);
    psi =0;
    a1(1:N-T) = 1/2*(a(T+1:end) + 1i*a(1:N-T)*1i*exp(1j*psi));
    a1(N-T+1:end) = a1(1:T); % sinnvolle Daten am Ende
    a2(1:N-T) = 1/2*(1i*a(T+1:end) + a(1:N-T)*1i*exp(1j*psi));
    a2(N-T+1:end) = a2(1:T); % sinnvolle Daten am Ende
    i_out = zeros(1,N);
  
    if intensity_plot == 1
      i_out = abs(a2).^2 - abs(a1).^2;
    else
      i_out = abs(a2) - abs(a1);
    end
  
    rx_electrical_signal_data = i_out;
  else
    if intensity_plot == 1
      rx_electrical_signal_data = abs(rx_signal_data).^2;
    else
      rx_electrical_signal_data = abs(rx_signal_data);
    end
 end

  step = size(rx_signal_data,2)/M;
  xlim_a =rx_signal_time(1+step/2)*1e9;
  xlim_b =rx_signal_time(step*2+step/2+1)*1e9;

  if strcmp(pulse_mode, 'nrz_2_dpsk') || strcmp(pulse_mode, 'rz_2_dpsk')
    rectangle("Position",[xlim_a,-1, xlim_a, 2],"FaceColor",[0.7 0.7 0.7],"EdgeColor",[0.7 0.7 0.7]);
    hold on;
    rectangle("Position",[xlim_b-xlim_a,-1, xlim_a, 2],"FaceColor",[0.7 0.7 0.7],"EdgeColor",[0.7 0.7 0.7]);
  else
    rectangle("Position",[xlim_a,0, xlim_a, 1.4],"FaceColor",[0.7 0.7 0.7],"EdgeColor",[0.7 0.7 0.7]); 
    hold on;
    rectangle("Position",[xlim_b-xlim_a,0, xlim_a, 1.4],"FaceColor",[0.7 0.7 0.7],"EdgeColor",[0.7 0.7 0.7]);
  end
 
  if color_grading == 0
    for m=3:M-3
      a = 1+step*(m-1)+step/2;
      b = a+step*2-1;  
      if strcmp(pulse_mode, 'nrz_2') || strcmp(pulse_mode, 'nrz_4') || strcmp(pulse_mode, 'duobinary') || strcmp(pulse_mode, 'nrz_2_dpsk')% NRZ
	a = 1+step*(m-1);
	b = a+step*2-1;
      end
      plot(rx_signal_time(1+step/2:step*2+step/2)*1e9, rx_electrical_signal_data(a:b), "linewidth", const_lw, "color", "m");
      hold on;
    end
 
    if strcmp(pulse_mode, 'nrz_2_dpsk') || strcmp(pulse_mode, 'rz_2_dpsk')
      ylim([-1 1]);
    else
      ylim([0 1.4]);
    end
    xlim([xlim_a xlim_b]);
    xlabel('time t/ns');
    if strcmp(pulse_mode, 'nrz_2_dpsk') || strcmp(pulse_mode, 'rz_2_dpsk')
      if intensity_plot == 1
	ylabel('differential photo current');
      else
	ylabel('differential photo current');
      end
    else
      if intensity_plot == 1
	ylabel('optical intensity');
      else
	ylabel('optical magnitude');
      end
    end  
  else
    hold off
    offset = step/2
    if strcmp(pulse_mode, 'nrz_2') || strcmp(pulse_mode, 'nrz_4') || strcmp(pulse_mode, 'duobinary') || strcmp(pulse_mode, 'nrz_2_dpsk') % NRZ
      offset = 0;
    end  
    res = create_matrix(rx_electrical_signal_data, step*2, 100, offset);
    
    imagesc(res);
    my_color_map = jet; 
    %my_color_map = flipud(my_color_map);
    my_color_map(1, :) = [1 1 1]; % set the 'blue' value (first row) to 'white'
    colormap(my_color_map);
    %colormap(flipud(colormap));
  end  

  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, pulse_mode, t1, nb_points);
  N=size(tx_data,2);
  tx_signal_data = zeros(1,nb_points);
  tx_signal_time = linspace(0,N*t1,nb_points);

  N1 = N;
  k_step = 1;
  if strcmp(pulse_mode, 'nrz_2') || strcmp(pulse_mode, 'nrz_4') || strcmp(pulse_mode, 'duobinary') || strcmp(pulse_mode, 'nrz_2_dpsk') % NRZ
    N1 = N/2;
    k_step = 1;  
  end
  tx_symbol_1 = zeros(1,nb_points/N1);
  tx_symbol_2 = zeros(1,nb_points/N1);
  tx_symbol = zeros(1,nb_points/N1);

  for k=1:k_step:N % symbol index
    for time=0:nb_points/N1-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/N1/2)/(nb_points/N1/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/N1/2)/(nb_points/N1/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;
	case 6
	  t_n = (time - nb_points/N1/2)/(nb_points/N1/2);
	  alpha = 0.0000000000001;
	  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
    if strcmp(pulse_mode, 'nrz_2') || strcmp(pulse_mode, 'nrz_4') || strcmp(pulse_mode, 'duobinary') || strcmp(pulse_mode, 'nrz_2_dpsk')% NRZ
      a = nb_points/N; %nb_points/N1*3/4
   
      tx_signal_data(1,1+(k-1)*nb_points/N:k*nb_points/N) = [tx_symbol_2(1,a+1:end) + tx_symbol(1,1:a)];   
      tx_symbol_1 = tx_symbol_2;
      tx_symbol_2 = tx_symbol;
    else
      tx_signal_data(1,1+(k-1)*nb_points/N:k*nb_points/N) = tx_symbol(1,:);
    end
  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

  if strcmp(pulse_mode, 'nrz_2') || strcmp(pulse_mode, 'nrz_4') || strcmp(pulse_mode, 'duobinary')|| strcmp(pulse_mode, 'nrz_2_dpsk')% NRZ
    tx_signal_data = tx_signal_data;
  end

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




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




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create_matrix
function res = create_matrix(data_signal, x_resolution, y_resolution, offset)
  res = zeros(y_resolution, x_resolution);
  data_signal = round(data_signal*y_resolution/(max(data_signal)+0.1));%*1.4/y_resolution;
  N = size(data_signal,2)/x_resolution;
  for k=1:N-1
    for index = 1:x_resolution  
     res(y_resolution - data_signal(1,index+(k-1)*x_resolution+offset), index) += 1;  
    end 
  end
end

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


