#d = [ 1 0 1 0 1 0 0 1 0 1 1 ];
datalen = 30;
d = rand(1,datalen) < 0.5;
n = 0:datalen-1;
Tb = 0.2;
fa = 2000;
freq = [ 50 30 ];
signallen = (datalen*Tb-1/fa)
t = 0:(1/fa):signallen;
f = freq(d(floor(t/Tb)+1)+1);
phi = cumsum(2*pi*f)/fa;
s = cos(phi);

subplot(4,1,1)
stem(n,d)
subplot(4,1,2)
plot(t,s)
axis([0 signallen -1 1])

[b,a] = butter(5,2*mean(freq)/fa);

slow = filter(b,a,s);
st = abs(slow);

[b2,a2] = butter(10,2*(min(freq)*0.8)/fa);
hk = filter(b2,a2,st);

subplot(4,1,3)
cla; hold on;
plot(t,slow,'k')
plot(t,st,'g')
plot(t,hk,'r')
hold off
axis([0 signallen -1 1])

mt = zeros(size(n));
for i = 1:datalen
  ts = (i-1) * Tb;
  te = i * Tb;
  tf = t > ts & t < te;
  mt(i) = mean(hk(tf));
endfor
hold on;
stem((n+0.5)*Tb,mt,'r');
hold off;

thr = (min(mt) + max(mt)) / 2;

ds = mt > thr;

subplot(4,1,4)
err = ds != d;
cla; hold on;
stem(n(~err),ds(~err))
stem(n(err),ds(err),'rx')
hold off;