function eval_midi_sync(pe, fe)
N = 5000;
t = (0:N-1)/N;

x1 = mod(20*t, 1) >=0.5;
x2 = mod((20+fe)*t+pe, 1) >=0.5;

dx1 = 50/N
dx2 = (50+fe)/N

qi1 = 0;
qi2 = 0;
xi1 = 0;
xi2 = 0;
xx1 = 0;
xx2 = 0;

lead = 0;
lag = 0;
ddx = 0;
for i=1:N,
%[qi1, qi2, xi1, xi2] = pd(x1(i), x2(i), qi1, qi2, xi1, xi2);
[qi1, qi2, xi1, xi2] = pd(xx1 >=0.5, xx2 >=0.5, qi1, qi2, xi1, xi2);
err = (qi1-qi2);
lead = 0.*err;
lag = lag + 0.000001*err;
ddx = (lead+lag);
xx1 = mod(xx1 + dx1, 1);
xx2 = mod(xx2 + dx2 + ddx, 1);


q1(i) = qi1;
q2(i) = qi2;

end;

dx2 = dx2 + ddx
dx1 = dx1

subplot(3, 1, 1)
plot(t, x1, t, x2); grid; axis([0 1 -0.2 1.2])
subplot(3, 1, 2)
plot(t, q1); grid; axis([0 1 -0.2 1.2])
subplot(3, 1, 3)
plot(t, q2); grid; axis([0 1 -0.2 1.2])


function [Q1, Q2, xo1, xo2] = pd(x1, x2, qi1, qi2, xi1, xi2)
xo1 = xi1;
xo2 = xi2;
qo1 = qi1;
qo2 = qi2;
for i=1:length(x1)
   % Q1 
   % Detect edge
   if (x1(i) - xo1) > 0.5
       qo1 = 1;
   end;
   xo1 = x1(i);
   
   % Q2 
   % Detect edge
   if (x2(i) - xo2) > 0.5
       qo2 = 1;
   end;   
   xo2 = x2(i);
   
   % Asynchronous reset
   if (qo1 * qo2) > 0.5
       qo1 = 0;
       qo2 = 0;
   end;

   Q1(i) = qo1;
   Q2(i) = qo2;
   
end;