% Demonstration of the detection of peak overlap by derivatives.
% To notice a change in any derivative caused by an overlapping peak, you'd
% have to know the shape of the derivative without the overlapping peak,
% and that gets more complicated as the derivative order increases.
% 
% To help in making such decisions, here a simple script that creates and
% plots the first five derivatives of two overlapping peaks. You can change
% the peak widths and separation and relative height of peak 2 in lines
% 44-50. In figure(1), you can clearly see that the central maximum of the
% even derivatives gets narrower at higher derivative orders. Nevertheless,
% detecting the overlapping peak visually does not necessarily get any
% easier at the higher orders. The problem is that the shapes get more
% complex and it's hard to remember exactly what the shape would have been
% for only a single peak. The first derivative is the simplest and easiest
% to remember, so one could make the argument that it's easiest to know
% when the peak is not a single peak by looking at the first derivative.
%
% To make this point clearer, look at figure 2. This is the same as
% figure(1) except that the derivatives of the original single peak 1 (in
% blue) are overlayed with the derivatives of the pairs pf peaks (in red).
% Even with this aid, it's not clear that the detection is any easier in
% the higher derivatives. On top of that, if there is any noise at all in
% the signal (line 50), the higher derivatives would need to be smoothed,
% which changes their shape slightly.  Of course, you might not have the
% luxury of having a single peak shape available for comparison, in which
% case you'd have to make do with figure(1) only.
%
% With a resolution of 0.346 (separation equals the half-widths of the
% peaks), the limit of visual detection for the second peak is a height of
% 0.1 of the larger peak. If the peaks are much closer, say position2=5,
% which makes the resolution 0.173, and with the height ratio 0.8, there is
% a change in shape observable in the 4th and 5th derivative orders that is
% not so easily seen in the lower orders, but the signal would need to have
% a very good signal-to-noise ratio in order for the 4th and 5th
% derivatives to be useful.
%
% Tom O'Haver, 2019.
% Requires gaussian.m, derivxy.m and fastsmooth.m in the path, download
% from http://tinyurl.com/cey8rwh
%
format compact
format short g

% Peak parameters you can change
increment=.01; % Difference between adjacent x values (1/sampling rate)
width=2; % Initial peak width (FWHM).
position1=4; % x-axis position of peak 1
position2=5.5;
h1=1; % Peak heights
h2=0.1;
noise=.00000001; % Random white noise standard deviation (optional)

% Calculations
sigma=width/(2*log(2)); % peak width (sigma)
Resolution=(position2-position1)./(4*sigma); % Original 4-sigma resolution
x=-1:increment:13;
y1=h1.*gaussian(x,position1,width); % Peak 1 (fixed height)
y2=h2.*gaussian(x,position2,width);% Peak 2 (variable height)
z=y1+y2;
y=z+noise*randn(size(x)); % Overlapping peak 1 and 2 signal with noise
y=fastsmooth(y,5,4); % Some smoothng requiredin noise>0 in line 50

figure(1)
subplot(3,2,1)
plot(x,y1,'b');
xlabel('x (time)')
ylabel('y');grid
title(['Gaussian peaks, Resolution= ' num2str(Resolution)  ', Height ratio= ' num2str(h2/h1) ])

subplot(3,2,2)
dy=derivxy(x,y);
plot(x,dy);
xlabel('x (time)')
ylabel('dy/dx');grid
title('First derivative')

subplot(3,2,3)
d2=derivxy(x,dy);
plot(x,d2);
xlabel('x (time)')
ylabel('d2y/dx^2');grid
title('Second derivative ')

subplot(3,2,4)
d3=derivxy(x,d2);
plot(x,d3);
xlabel('x (time)')
ylabel('d3y/dx^3');grid
title('Third derivative ')

subplot(3,2,5)
d4=derivxy(x,d3);
plot(x,d4);
xlabel('x (time)')
ylabel('d4y/dx^4');grid
title('Fourth derivative ')

subplot(3,2,6)
d5=derivxy(x,d4);
plot(x,d5);
xlabel('x (time)')
ylabel('d5y/dx^5');grid
title('Fifth derivative ')

figure(2)
subplot(3,2,1)
plot(x,y1,'b',x,y2,'g',x,y,'r');
xlabel('x (time)')
ylabel('y')
title(['Gaussian peaks, Resolution= ' num2str(Resolution)  ', Height ratio= ' num2str(h2/h1) ])

subplot(3,2,2)
dy=derivxy(x,y);
dy1=derivxy(x,y1);
plot(x,dy1,x,dy);
xlabel('x (time)')
ylabel('dy/dx')
title('First derivative')

subplot(3,2,3)
d2=-derivxy(x,dy);
d21=-derivxy(x,dy1);
plot(x,d21,x,d2);
xlabel('x (time)')
ylabel('-d2y/dx^2')
title('Negative of second derivative ')

subplot(3,2,4)
d3=derivxy(x,-d2);
d31=derivxy(x,-d21);
plot(x,d31,x,d3);
xlabel('x (time)')
ylabel('d3y/dx^3')
title('Third derivative ')

subplot(3,2,5)
d4=derivxy(x,d3);
d41=derivxy(x,d31);
plot(x,d41,x,d4);
xlabel('x (time)')
ylabel('d4y/dx^4')
title('Fourth derivative ')

subplot(3,2,6)
d5=derivxy(x,d4);
d51=derivxy(x,d41);
plot(x,d51,x,d5);
xlabel('x (time)')
ylabel('d5y/dx^5')
title('Fifth derivative ')