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

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

%define AP 2
AP2amp = .7;
AP2dur1 = .0015;
AP2dur2 = .003;
AP2prob = .002;
%make single AP
%by concantenating two 1/2 sin cycles
oneT = 0:dt:(AP2dur1+AP2dur2);
oneAP2 = AP2amp*sin(2*pi/(2*AP2dur1)*oneT).*(oneT<=AP2dur1) + ...
    AP2amp*sin(2*pi/(2*AP2dur2)*(oneT+AP2dur1)).*(oneT>AP2dur1) ;

% figure(1)
% plot(oneAP2)
% return
%One cycle
% AP1dur1 = .004;
% AP1dur2 = .003;
% oneT = 0:dt:(AP1dur1);
% oneAP = AP1amp*sin(2*pi/(AP1dur1)*oneT) ;

%make a bunch of APs
locationAP1 = rand(size(t))<AP1prob; %poisson
%locationAP1 = mod([1:length(t)],457)==0; %uniform
vAP = conv(locationAP1,oneAP1);
vAP = vAP(1:length(t)); %truncate
locationAP2 = rand(size(t))<AP1prob; %poisson
vAP2 = conv(locationAP2, oneAP2*AP2amp);
vAP = vAP + vAP2(1:length(t)); %truncate

%make gaussian noise plus 60 Hz interference
noiseAmp = .01;
interAmp = .05;
phase = -1.0 ;
vnoise = noiseAmp*randn(size(t)) + interAmp*sin(60*2*pi*t+phase);
vref = sin(60*2*pi*t) ;
%now add noise to APs to simulate a real signal 
vAPnoise =  vAP + vnoise;


nplot=2
figure(10);clf;

subplot(nplot,1,1)
plot(t,vAP)
title('poisson APs')
set(gca,'ylim',[-1 1.5])

subplot(nplot,1,2)
plot(t,vAPnoise)
hold on
plot(t, vref/4, 'r')
title('poisson APs with noise and 60 Hz ref in red')
set(gca,'ylim',[-1 1.5])



%