%
%
%
%define time
Fs = 10000;
dt = 1/Fs;
t = 0:dt:1;

%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

%make gaussian noise plus 60 Hz interference
noiseAmp = 1;
interAmp = 1;
vnoise = noiseAmp*randn(size(t)) + interAmp*sin(60*2*pi*t);

%now add niose to APs and 
vtotal = vnoise + vAP ;


nplot=5

subplot(nplot,1,1)
pwelch(vnoise,[],[],[],Fs)
title('Filtering noise')
legend('input noise: RMS=60Hz')
set(gca,'xlim',[0 2000])

subplot(nplot,1,2)
%Median filter
%vmedian=medfilt1(vnoise,7);
vlowmedian = medfilt1(vnoise,length(oneAP));
vfilt = vnoise - vlowmedian;
vbogus = vfilt;

[b,a] = butter(4, 500/(Fs/2)); %3 KHz, 4-pole lowpass
vfilt = filter(b,a,vfilt);
[b,a] = butter(2, 200/(Fs/2), 'high'); %200 Hz, 2-pole highpass
vfilt = filter(b,a,vfilt);
title('median-flitered/butterworth spectrum')
pwelch(vfilt,[],[],[],Fs)
set(gca,'ylim',[-100 0])
set(gca,'xlim',[0 2000])
legend('median-flitered/bu+tterworth spectrum')

subplot(nplot,1,3)
[b,a] = butter(4, 500/(Fs/2)); %3 KHz, 4-pole lowpass
vfilt = filter(b,a,vnoise);
[b,a] = butter(2, 200/(Fs/2), 'high'); %200 Hz, 2-pole highpass
vfilt = filter(b,a,vfilt);
title('butterworth-flitered spectrum')
pwelch(vfilt,[],[],[],Fs)
set(gca,'ylim',[-100 0])
set(gca,'xlim',[0 2000])
legend('butterworth-flitered spectrum')

subplot(nplot,1,4)
vfilt = conv(vnoise,oneAP)/length(oneAP);
pwelch(vfilt,[],[],[],Fs)
set(gca,'ylim',[-100 0])
set(gca,'xlim',[0 2000])
legend('AP matched-flitered spectrum')

subplot(nplot,1,5)
vfilt = conv(vbogus,oneAP)/length(oneAP);
pwelch(vfilt,[],[],[],Fs)
set(gca,'ylim',[-100 0])
set(gca,'xlim',[0 2000])
legend('AP matched-flitered/median spectrum')


%