% test 60 Hz canceller
clear all
figure(1); clf;

Fs = 48000;
dt = 1/Fs;
t = 0:dt:0.5;

%reference noise
ref = sin(2*pi*60*t);
refshift = sin(2*pi*60*t + pi/2);
ph = 90;
%signal
% noisefree = sin(2*pi*10*t);
% sig = 1.0* sin(2*pi*60*t+ph) + noisefree ;


%define APs
AP1amp = 1;
AP1dur1 = .001;
AP1dur2 = .002;
AP1prob = .001;
AP2amp = AP1amp*.5;

%make single AP
%by concantenating two 1/2 sin cycles
oneT = 0:dt:(AP1dur1+AP1dur2);
oneAP = AP1amp*sin(2*pi/(2*AP1dur1)*oneT).*(oneT<=AP1dur1) + ...
    AP1amp*0.5*sin(2*pi/(2*AP1dur2)*(oneT+AP1dur1)).*(oneT>AP1dur1) ;

%make a bunch of APs
locationAP1 = rand(size(t))<AP1prob; %poisson
%locationAP1 = mod([1:length(t)],457)==0; %uniform
vAP = conv(locationAP1,oneAP);
vAP = vAP(1:length(t)); %truncate
locationAP2 = rand(size(t))<AP1prob; %poisson
vAP2 = conv(locationAP2,oneAP*AP2amp);
vAP = vAP + vAP2(1:length(t)); %truncate
noisefree = vAP+ 0.00*randn(size(t));
sig = noisefree +1* sin(2*pi*60*t+1)  ;

%%%%%%%%%%%%%%%%%%%%%%%%%
w1 = 0;
w2 = 0;
mu2e =2^-11 ;

for i=1:length(t)
    out(i) = sig(i) -  (w1*ref (i)+ w2*refshift(i)) ;
    w1 = w1 + mu2e*out(i)*sign(ref(i)) ;
    w2 = w2 + mu2e*out(i)*sign(refshift(i)) ;
end

np=2;

subplot(np,1,1)
plot(t,sig);

subplot(np,1,2)
plot(t,out)
hold on
plot(t,noisefree,'r')


