% define genetic search parameter
generations = 5;   % number of generations
N = 500;          % population size
M = 10*N;           % number of offspring
pmut = .05;         % probability of mutation
std = .2;            % std of random mutations
delay = 2.5;          % delay of equalized signal path in sample times
D = ceil(delay);

% define baud rate
fbaud = 125e6;
Tbaud = 1/fbaud;

% frequency scale
fscale = 2*pi*70e6;

% initialization stuff
coax_mdl_noplot;  % setup channel model
seed = sum(100*clock);
rand('state',seed);
clear pop offspring e;

% create initial population
for k = 1:N
   % random orthonormal ladder
   b = (rand(1,3)+.1).*sign(rand(1,3)-.5);
   pop(k,:) = b;
end;

% perform evolution
for k = 1:generations

  % show current generation 
  k

  % PZ map of current population
  if (sum(k == [0]))
  figure
  for l = 1:N
    pzmap(ol(fscale*[pop(l,1) 0 pop(l,2) 0 pop(l,3)],fscale*a));
    hold on
  end;
  hold off
  axis([-8 +8 -8 +8]); 
  title(sprintf('Generation %d',k));
  end;

  % calculate fitness function for each individual
  for l = 1:N
    imp = impulse(cable_sys_s*ol(fscale*[pop(l,1) 0 pop(l,2) 0 pop(l,3)],fscale*a),Tbaud*(mod(delay,1):delay+30)); 
    e(l) = sum([imp(1:D-1);imp(D+1:end)].^2)/imp(D)^2; 
  end;
  emin = min(e)

  % create offspring
  numoffspring = round(M*(1./e)/sum(1./e)); 
  m = 0;
  for l = 1:N 
    for x = m+1:m+numoffspring(l);
      offspring(x,:) = pop(l,:);
    end;
    m = m + numoffspring(l);
  end;

  % create next generation
  l = find(e==emin);
  pop(1,:) = pop(l(1),:);
  for l = 2:N
    parent1 = ceil(m*rand(1));
    parent2 = ceil(m*rand(1));
    brk = ceil(3*rand(1));
    pop(l,:) = [offspring(parent1,1:brk) offspring(parent2,brk+1:end)];
  end;

  % mutation
  for l = 2:N
    if (rand(1) < pmut)
      pop(l,:) = pop(l,:) + std*randn(1,3);
    end;
  end;

end;

% calculate fitness function for each individual
for l = 1:N
  imp = impulse(cable_sys_s*ol(fscale*[pop(l,1) 0 pop(l,2) 0 pop(l,3)],fscale*a),Tbaud*(mod(delay,1):delay+30)); 
  e(l) = sum([imp(1:D-1);imp(D+1:end)].^2)/imp(D)^2; 
end;

% choose the best of the survivors
emin = min(e)
k = find(e==emin);  k = k(1);

% plot impulse response
bopt = [pop(k,1) 0 pop(k,2) 0 pop(k,3)]
sys = ol(fscale*[pop(k,1) 0 pop(k,2) 0 pop(k,3)],fscale*a);
chan_imp = impulse(cable_sys_s,Tbaud*(mod(delay,1):delay+30));
eq_imp = impulse(cable_sys_s*sys,Tbaud*(mod(delay,1):delay+30));
hold off; stem(Tbaud*(mod(delay,1):delay+30),chan_imp,'r');
hold on; stem(Tbaud*(mod(delay,1):delay+30),eq_imp,'g'); 
impulse(cable_sys_s*sys,Tbaud*(delay+30),'g');
hold off;

