

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

%=== settable parameters =======================
lowpass_cutoff = 25; %Hz
highpass_cutoff = 1; %Hz
noise_amp = .2 ; % thermal noise relative to QRS
interference_amp = 100; % 60 Hz nose amp relative to QRS
mu2e =.002 ; % LMS convergence parameter 
baseline_amp = 1 ; %baseline drift amplitude
Fs = 360*8; % sample rate
%90 degrees at 60 Hz is 12 samples
%60 degress at 60 Hz is 8 samples
%30 degress at 60 Hz is 4 samples
%=== end settable parameters ===================

dt = 1/Fs;
t = 0:dt:15 ;

%reference noise
f = 60.5 ;
ref = sin(2*pi*f*t) +0.1*sin(2*pi*3*f*t);

%define QRS
AP1amp = 1;
AP1dur1 = .050;

%make single QRS
oneT = 0:dt:(AP1dur1);
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 QRS
%locationAP1 = double(rand(size(t))<AP1prob); %poisson
locationAP1 =  double(0==mod(t,0.95));
%locationAP1 = mod([1:length(t)],457)==0; %uniform
vAP = conv(locationAP1,oneAP);
vAP = vAP(1:length(t)); %truncate
noisefree = vAP+ noise_amp*randn(size(t));
sig = noisefree + baseline_amp * sin(2*pi*0.2*t);
sig =sig + interference_amp * sin(2*pi*f*t+1)+ interference_amp*0.1*sin(2*pi*3*f*t+3) ;

[b,a] = butter(4, lowpass_cutoff/(Fs/2));
sig1 = filter(b, a, sig );

[b,a] = butter(2, highpass_cutoff/(Fs/2),'high');
sig2 = filter(b, a, sig1 );

%%%%%%%%%%%%%%%%%%%%%%%%%
n = 20 ;
w = zeros(1,n+1); % the weight vector
offsets = 0:-1:-n ; % offsets into the reference signal samples

for i=n+1:length(t)
    history = ref(i+offsets);
    out(i) = sig2(i) - dot(w,history) ;
    w = w + mu2e * out(i) * history ;
end

np=2;

subplot(np,1,1)
plot(t,sig);
hold on
plot(t,sig2, 'r');
legend('unfiltered', 'filtered(before LMS)')

subplot(np,1,2)
plot(t,out)
hold on
plot(t,vAP-2,'r')
set(gca, 'ylim', [-2 2])
legend('with LMS', 'original')


