function [signal_rate, num_bits] = assignment_example ()
  // Simple example of Scilab for Assignment
  // Calling Sequence
  // [signal_rate, num_bits] = assignment_example();
  // Authors
  // Steven Gordon; steven.gordon@ieee.org
  
  // Create a string
  name = 'Steve';
  
  // Convert to binary
  name_bin = text2bin(name);
  
  // Count the number of bits
  num_bits = size(name_bin,2);
  
  // Encode using NRZ-L (data rate = 2bps)
  [out_signal,signal_rate] = ds_nrzl_encode(name_bin,2);
  
  // Plot the output signal
  as_plot(out_signal);
  
  // Save the plot to a file (PNG)
  xs2png(0,name+'-example.png');
  clf;
endfunction


function [data] = text2bin (text)
  // Converts ASCII text to an array of bits
  // Calling Sequence
  // binary_data = text2bin('Steve')
  // Parameters
  // text: string of ASCII characters
  // data: array of bits
  // Authors
  // Steven Gordon; steven.gordon@ieee.org

  a = dec2bin(ascii(text),7);
  for i=1:7,
    b=part(a,i);
    for j=1:size(b,2),
      data(1,(j-1)*7+i)=msscanf(b(j),'%d');
    end
  end
endfunction


function [text] = bin2text(data)
  // Converts array of bits to ASCII text
  // Calling Sequence
  // string_St = bin2text([ 1 0 1 0 0 1 1 1 1 1 0 1 0 0])
  // Parameters
  // data: array of bits
  // text: String of ASCII characters
  // Authors
  // Steven Gordon; steven.gordon@ieee.org
  
  bits = size(data,2);
  // Input data should be a multiple of 7 bits
  // If not add 0's at the end
  if floor(bits/7)*7<>bits,
    data = [data zeros(((floor(bits/7)+1)*7)-bits)];
    bits = size(data, 2);
  end
   
  text = ''; 
  for i=1:7:bits,
    x = strcat(string(data(1,i:(i+6))));
    text = strcat([text ascii(bin2dec(x))]);
  end
endfunction


function [output_signal,signal_rate] = ds_nrzl_encode (input_data, data_rate)
  // Generates a digital signal from input bits using Non-Return to Zero Level (NRZ-L)
  // Calling Sequence
  // [os,sr] = ds_nrzl_encode ([1 0 1 0 0 1 1], 1000)
  // Parameters
  // input_data: array of bits for input data, e.g. [0 1 0 0 1]
  // data_rate: rate at which output data is generated in bits per second, e.g. 2000
  // output_signal: matrix. Column 1 time; Col 2 signal value
  // signal_rate: rate at which signal elements are sent
  // Authors
  // Steven Gordon; steven.gordon@ieee.org
  
  level_1 = -1;
  level_2 = 1;
  points = global_points();
  
  t=0; // Initial value of time
  signal_rate = data_rate;
  element_duration = 1/signal_rate;
  point_duration = element_duration/points;
  
  for i=1:size(input_data,2),
    if input_data(1,i)==0, 
      // If bit is 0
      for j=1:points,
        // Signal is high
        output_signal((i-1)*points+j,:)=[t+(j-1)*point_duration level_2]; 
      end
    else
      // If bit is 1
      for j=1:points,  
        // Signal is low
        output_signal((i-1)*points+j,:)=[t+(j-1)*point_duration level_1];
      end
    end
    t = t+element_duration;
  end
endfunction


function [output_data,data_rate]=ds_nrzl_decode (input_signal,signal_rate)
  // Generates data from received digital signal using NRZ-L
  // Calling Sequence
  // [od,dr] = ds_nrzl_decode (os,100)
  // Parameters
  // input_signal: matrix. Column 1 time; Col 2 signal value
  // signal_rate: rate at which signal elements are sent
  // output_data: array of bits for data, e.g. [0 1 0 0 1]
  // data_rate: rate at which output data is generated in bits per second, e.g. 2000
  // Authors
  // Steven Gordon; steven.gordon@ieee.org
  
  level_1 = -1;
  level_2 = 1;
  points = global_points();
  data_rate = signal_rate;

  // Repeat for all signal elements
  for i=1:size(input_signal,1)/points,
    // Sample the signal at the middle of the interval
    if (input_signal((i-1)*points+1+floor(points/2),2)) >= 0,
      // If signal is high
      // Bit is 0
      output_data(1,i)=0;
    else
      // If signal is low
      // Bit is 1
      output_data(1,i)=1;
    end
  end
endfunction


function [output_signal,signal_rate] = ds_nrzinvert_encode (input_data, data_rate)
  // Generates a digital signal from input bits using Non-Return to Zero Invert on Ones
  // Calling Sequence
  // [os,sr] = ds_nrzinvert_encode ([1 0 1 0 0 1 1], 200)
  // Parameters
  // input_data: array of bits for input data, e.g. [0 1 0 0 1]
  // data_rate: rate at which output data is generated in bits per second, e.g. 2000
  // output_signal: matrix. Column 1 time; Col 2 signal value
  // signal_rate: rate at which signal elements are sent
  // Description
  // Assumes signal is low initially
  // Authors
  // Steven Gordon; steven.gordon@ieee.org
  
  level_1 = -1;
  level_2 = 1;
  points = global_points();
  
  t=0; // Initial value of time
  signal_rate = data_rate;
  element_duration = 1/signal_rate;
  point_duration = element_duration/points;
  
  current_level = level_1; // Assume signal is initially low
  for i=1:size(input_data,2),
    if input_data(1,i)==0, 
      // If bit is 0
      for j=1:points,
        // Do not change the signal
        output_signal((i-1)*points+j,:)=[t+(j-1)*point_duration current_level]; 
      end
    else
      // If bit is 1
      // Change the level
      if current_level == level_1,
        current_level = level_2;
      else
        current_level = level_1;
      end
      for j=1:points,  
        output_signal((i-1)*points+j,:)=[t+(j-1)*point_duration current_level];
      end
    end
    t = t+element_duration;
  end
endfunction


function [output_data,data_rate]=ds_nrzinvert_decode (input_signal,signal_rate)
  // Generates data from received digital signal using NRZ Invert on ones
  // Calling Sequence
  // [od,dr] = ds_nrzinvert_decode (os,100)
  // Parameters
  // input_signal: matrix. Column 1 time; Col 2 signal value
  // signal_rate: rate at which signal elements are sent
  // output_data: array of bits for data, e.g. [0 1 0 0 1]
  // data_rate: rate at which output data is generated in bits per second, e.g. 2000
  // Description
  // For each signal element, a sample is taken in the middle of the interval
  // Authors
  // Steven Gordon; steven.gordon@ieee.org
  
  level_1 = -1;
  level_2 = 1;
  points = global_points();
  data_rate = signal_rate; 
  current_level = level_1; // Assume signal is initially low
  
  // Repeat for all signal elements
  for i=1:size(input_signal,1)/points,
    // Sample the signal at the middle of the interval
    // E.g. each signal element has p points; take the p/2 value
    if (input_signal((i-1)*points+1+floor(points/2),2)) >= 0,
      // If signal is high
      if current_level==level_2,
        // Bit is 0 (no change)
        output_data(1,i)=0;
      else
        // Bit is 1
        output_data(1,i)=1;
      end
      current_level = level_2;
    else
      // If signal is low
      if current_level==level_1,
        // Bit is 0 (no change)
        output_data(1,i)=0;
      else
        // Bit is 1
        output_data(1,i)=1;
      end
      current_level = level_1;        
    end
  end
endfunction

  
function [digital_data,min_amp,max_amp] = pcm_encode(analog_data, sampling_rate, num_levels)
  // Convert analog data to digital data using Pulse Code Modulation encoding
  // Parameters
  // analog_data: matrix of analog input data, col 1 time, col 2 magnitude
  // sampling_rate: rate to sample [samples/sec]
  // num_levels: number of code levels
  // digital_data: array of bits
  // Authors
  // Steven Gordon; steven.gordon@ieee.org

  min_amp = min(analog_data(:,2));
  max_amp = max(analog_data(:,2));
  max_time = max(analog_data(:,1));
  code_points = [min_amp:(max_amp-min_amp)/num_levels:max_amp-max_amp/num_levels];
  sample_interval = 1/sampling_rate;
   
  // Assume the analog data starts at time 0
  sample_points = [0:sample_interval:max_time];
  
  s=1;
  for i=sample_points,
    // Find the value for the current sample time
    index = find (analog_data(:,1)==i,1);
    if size(index,1)==0,
      next_index = find(analog_data(:,1)>i,1);
      if next_index > 1,
        closeness = abs(analog_data(next_index,1)-i);
        index = next_index;
        if ((analog_data(next_index-1,1)-i)<closeness),
          index = next_index-1;
        end
      else
        index = next_index;
      end
    end
    current_sample = analog_data(index,2);
    code_level = find(code_points>current_sample,1)-2;
    if code_level==-2,
      code_level=num_levels-1;
    end
    x = dec2bin(code_level,log2(num_levels));
    for j=1:log2(num_levels),
      b=part(x,j);
//      digital_data((s-1)*log2(num_levels)+j,1) = i;
      digital_data((s-1)*log2(num_levels)+j,1)=msscanf(b,'%d');
    end
    s = s+1;
  end
  digital_data = digital_data';
endfunction


function [analog_data] = pcm_decode(digital_data, sampling_rate, num_levels, min_amp, max_amp)
  // Convert digital data to analog data using Pulse Code Modulation decoding
  // Parameters
  // digital_data: array of bits
  // sampling_rate: rate to sample [samples/sec]
  // num_levels: number of code levels
  // analog_data: matrix of analog input data, col 1 time, col 2 magnitude
  // Authors
  // Steven Gordon; steven.gordon@ieee.org
  
  code_points = [min_amp:(max_amp-min_amp)/num_levels:max_amp-max_amp/num_levels];
  sample_interval = 1/sampling_rate;
  max_time = sample_interval*size(digital_data,2)/log2(num_levels)-sample_interval;
  points = global_points();
  sample_points = [0:sample_interval:max_time];
  j=1;
  for i=sample_points,
    x = digital_data(1,((j-1)*log2(num_levels)+1):((j-1)*log2(num_levels)+log2(num_levels)));
    code_level = bin2dec(strcat(string(x)))+1;
    signal_value = code_points(code_level);
    for k=[1:points],
      analog_data((j-1)*points+k,1) = i+(k-1)*(sample_interval/points);
      analog_data((j-1)*points+k,2) = signal_value;
    end
    j=j+1;
  end
    
endfunction


function [analog_signal] = as_create (A, f, p, n, duration)
  // Create an analog signal composed of n sine waves
  // Parameters
  // A: peak amplitude
  // f: frequency
  // p: phase [radians]
  // n: number of sine wave components (must be 1 or greater)
  // duration: time of signal [seconds]
  // analog_signal: matrix of analog signal
  // Authors
  // Steven Gordon; steven.gordon@ieee.org
      
  period = 1/f;
  points=global_points();
  t=[0:period/points:duration];
  s=0;
  if n<1,
    error('n must be greater than 1');
  elseif n==1,
    s = sin(2*%pi*f*t+p);
  else
    for i = 1:2:(2*n-1),
      s = s+(1/i)*sin(2*%pi*(i*f)*t + p);
    end
  end
  s = A * s;
  analog_signal = [t' s'];
endfunction


function as_plot(signal)
  // Plots an analog signal 
  // Calling Sequence
  // as_plot(analog_signal)
  // Parameters
  // signal: matrix of signal values (col 1: time; col 2 value)
  // Authors
  // Steven Gordon; steven.gordon@ieee.org
  
  plot2d(signal(:,1),signal(:,2));
  // Set the axes limits
  a=gca();
  ymax=1.1*max(signal(:,2));
  ymin=1.1*min(signal(:,2));
  
  a.data_bounds = [0,max(signal(:,1)),ymin,ymax];
  a.x_label.text='time [sec]';
  a.y_label.text='signal';
  a.grid = [1 -1 -1];
endfunction


function [rx,n] = add_uniform_noise (tx,noise_amplitude)
  // Adds random (uniform) noise to a transmitted signal
  // Parameters
  // tx: Signal to be transmitted
  // noise_amplitude: magnitude of the noise
  // Authors
  // Steven Gordon; steven.gordon@ieee.org

  t = max(tx(:,1))+(tx(2,1)-tx(1,1));
  n = as_uniform_noise(noise_amplitude,t,size(tx,1)/t);
  rx = add_noise(tx,n);
endfunction


function [rx,impulse] = add_impulse_noise (tx,impulse_amplitude,impulse_duration)
  // Adds a single impulse of noise
  // Parameters
  // tx: Signal to be transmitted
  // impulse_amplitude: magnitude of the noise
  // impulse_duration: length of impulse in seconds
  // Authors
  // Steven Gordon; steven.gordon@ieee.org

  t = max(tx(:,1))+(tx(2,1)-tx(1,1));
  total_points = size(tx,1);
  // Background noise
  bg_n = as_uniform_noise(0.1,t,total_points/t);
  // Impulse noise
  if grand(1,1,'def')>0.5,
    sign_impulse = 1;
    low = 0.7*impulse_amplitude;
    high = 1.3*impulse_amplitude;
  else
    sign_impulse = -1;
    impulse_amplitude = sign_impulse * impulse_amplitude;
    low = 1.3*impulse_amplitude;
    high = 0.7*impulse_amplitude;
  end
    
  n = grand(impulse_duration*total_points/t,1,'unf',low,high);
  
  // Locate the impulse randomly in the tx signal
  index = grand(1,1,'uin',1,floor(total_points*(t-impulse_duration)/t));
  x = [tx(:,1) zeros(total_points,1)];
  for i=index:(index+size(n,1)-1),
    x(i,2)=n(i-index+1,1);
  end
  in = add_noise(bg_n,x);
  impulse=add_noise(bg_n,in);
  rx = add_noise(tx,in);
endfunction


function [ber,num_errors]=bit_error_rate(tx,rx)  
  // Calculate the bit error rate (and number of errors)
  num_bits = size(tx,2);
  num_errors = num_bits-sum(tx==rx);
  ber = num_errors/num_bits;
endfunction


function rx = add_noise (tx,n)
  // Adds noise to a transmitted signal

  // The times for both signal must be identical
  if size(tx)==size(n),
    rx = [tx(:,1) tx(:,2)+n(:,2)];
  else
    error('Tx and N must use same time scale');
  end
endfunction


function n=as_uniform_noise(amplitude,time,points_per_sec)
  // Creates random noise of the given amplitude
  r=grand(points_per_sec*time,1,'unf',-1*amplitude,1*amplitude);
  t=[0:1/points_per_sec:(time-1/points_per_sec)]';
  n=[t r];
endfunction


function [output_signal] = sine_element(A, f, p, starttime,endtime)
  // Generate a sinusoid for a single signal element
  points = global_points();
  duration = endtime-starttime;
  f = f/duration;
  t = [starttime:duration/points:endtime-(duration/points)];
  s = A*sin(2*%pi*f*t+p);
  output_signal=[t' s'];
endfunction
 
 
function [points]=global_points()
  // Defines the points used for accuracy in plots
  // Calling Sequence
  // global_points()
  // Authors
  // Steven Gordon; steven.gordon@ieee.org

  // You should not use a value less than 20.
  points = 40;
endfunction



