% Script to compare segmented-smooth-based peak finder findpeaksSG to the
% wavelet denoise-based peak finder findpeaksSGw. Requires findpeaksSG.m
% findpeaksSGw.m, gaussian.m (or other peak shape) function, and wdenoise.m
% from the Wavelet toolbox. Creates 100 noisy signals with 3 peaks with
% different widths, then measures the peak positions and heights, and
% computes the percent relative standard deviation of 100 measurements. The
% segmented peak finder findpeaksSG, with 3-segment peak detection
% parameters, reliably detects amd measures all three peaks. The
% denoise-based method often gives better slightly precision but sometimes
% misses peaks. Both methods give better precision for the wider peak
% (because it has more data points).
clear
clf
x=1:1000;
format compact
format short g
p=[200 400 800]; % Peak positions of the 3 peaks.
h=[1 1 1]; % Peak heights of the 3 peaks.
w=[10 50 200]; % Peak widths of the 3 peaks.
Noise=.2; % Standard deviation of the random noise
disp(['Standard deviation of noise: ' num2str(Noise) ])
% Define the peak detection parameters for each segment
SlopeThresholds=[.0001 .000001 .0000001];
AmpThresholds=[.2 .2 .2];
smoothwidths=[41 51 61]; % Smoothing for peak detection only
peakgroups=[10 50 250];
level=7;
NumTrials=200; % Number of independent signals measured
y=h(1).*gaussian(x,p(1),w(1))+h(2).*gaussian(x,p(2),w(2))+h(3).*gaussian(x,p(3),w(3));
% 
tic
for trial=1:NumTrials
    % Create the noisy signal
    ny(:,trial)=y+Noise.*randn(size(x));
    % Detect and measure peaks in the signal based on raw unsmoothed signal
    Ps=findpeaksSG(x,ny(:,trial),SlopeThresholds,AmpThresholds,smoothwidths,peakgroups,3);
    sizePs=size(Ps);
    NumPeaksS(trial)=sizePs(1);
    figure(1)
    if trial==1 % Plot the first signal
        plot(x,ny,'c'),
    end
    hold on
    for peak=1:NumPeaksS(trial)
        text(Ps(peak,2),Ps(peak,3),'.') % Number the peaks
        % Record the measured peak parameters of each peak in each trial
        PeakPositionsS(trial,peak)=Ps(peak,2);
        PeakHeightsS(trial,peak)=Ps(peak,3); 
    end
end
t1=toc;
hold off
disp('Segmented smoothing ')
xlabel('Each detected peak position and height is marked with a dot')
ylabel('Y')
title('Segmented smoothing')
disp(['Average number of peaks: ' num2str(mean(NumPeaksS)) ])
disp('Peak position and height errors and relative standard deviations')
disp('                 Peak 1       Peak 2        Peak 3')
PositionsErrorS=100.*(p-mean(PeakPositionsS(:,1:3)))./p;
disp(['% Position Error  ' num2str(PositionsErrorS) ])
PositionsRSDS=100.*stdiqr(PeakPositionsS(:,1:3))./mean(PeakPositionsS(:,1:3));
disp(['% Position RSD    ' num2str(PositionsRSDS) ])
HeightsErrorS=100.*(h-mean(PeakHeightsS(:,1:3)))./h;
disp(['% Height Error    ' num2str(HeightsErrorS) ])
HeightsRSDS=100.*stdiqr(PeakHeightsS(:,1:3))./mean(PeakHeightsS(:,1:3));
disp(['% Height RSD      ' num2str(HeightsRSDS) ])

disp(' ')
disp('Wavelet denoising method')
tic
for trial=1:NumTrials
    % Detect and measure peaks in the signal
    % CHOOSE A WAVELET from the four following lines
    % wdy=ny; if trial==1;disp('No denoising');end % No denoising;
    % wdy=wdenoise(ny(:,trial);if trial==1;disp('Default denoising');end % Default denoising
    wdy=wdenoise(ny(:,trial),level,'DenoisingMethod','BlockJS'); if trial==1;disp('BlockJS denoising');end
    % wdy=wdenoise(ny(:,trial),level,'Wavelet','coif2'); if trial==1;disp('coif2 denoising');edn
    
    Pw=findpeaksSGw(x,wdy,SlopeThresholds,AmpThresholds,peakgroups,level);
    sizePw=size(Pw);
    NumPeaksW(trial)=sizePw(1);
    figure(2)
    if trial==1 % Plot the first signal
        plot(x,ny(:,trial),'c'),
    end
    hold on
    for peak=1:NumPeaksW(trial)
        text(Pw(peak,2),Pw(peak,3),num2str('.')) % Number the peaks
        % Record the measured peak parameters of each peak in each trial
        PeakPositionsW(trial,peak)=Pw(peak,2);
        PeakHeightsW(trial,peak)=Pw(peak,3);
    end
end
t2=toc;
hold off
xlabel('Each detected peak position and height is marked with a dot')
ylabel('Y')
title('Wavelet denoising.')
disp(['Wavelet level: ' num2str(mean(level)) ])
disp(['Average number of peaks: ' num2str(mean(NumPeaksW)) ])
disp('Peak position and height errors and relative standard deviations')
disp('                 Peak 1       Peak 2        Peak 3')
PositionsErrorW=100.*(p-mean(PeakPositionsW(:,1:3)))./p;
disp(['% Position Error  ' num2str(PositionsErrorW) ])
PositionsRSDW=100.*stdiqr(PeakPositionsW(:,1:3))./mean(PeakPositionsW(:,1:3));
disp(['% Position RSD    ' num2str(PositionsRSDW) ])
HeightsErrorW=100.*(h-mean(PeakHeightsW(:,1:3)))./h;
disp(['% Height Error    ' num2str(HeightsErrorW) ])
HeightsRSDW=100.*stdiqr(PeakHeightsW(:,1:3))./mean(PeakHeightsW(:,1:3));
disp(['% Height RSD      ' num2str(HeightsRSDW) ])
TimeRatio=t2/t1;
disp(' ')
MeanAbsolutePositionErrorRatio=mean(abs(PositionsErrorS))./mean(abs(PositionsErrorW));
MeanAbsoluteHeightErrorRatio=mean(abs(HeightsErrorS))./mean(abs(HeightsErrorW));
if MeanAbsolutePositionErrorRatio>1
    disp(['Wavelet method better for accuracy of peak position by a factor of ' num2str(MeanAbsolutePositionErrorRatio) ] )
else
    disp(['Segmented smooth method better for accuracy of peak position by a factor of ' num2str(1/MeanAbsolutePositionErrorRatio) ] )
end
if MeanAbsoluteHeightErrorRatio>1
    disp(['Wavelet method better for accuracy of peak height by a factor of ' num2str(MeanAbsoluteHeightErrorRatio) ] )
else
    disp(['Segmented smooth method better for accuracy of peak height by a factor of ' num2str(1/MeanAbsoluteHeightErrorRatio)])
end

MeanAbsolutePositionRSDRatio=mean(abs(PositionsRSDS))./mean(abs(PositionsRSDW));
MeanAbsoluteHeightRSDRatio=mean(abs(HeightsRSDS))./mean(abs(HeightsRSDW));
if MeanAbsolutePositionRSDRatio>1
    disp(['Wavelet method better for precision of peak position by a factor of ' num2str(MeanAbsolutePositionRSDRatio)])
else
    disp(['Segmented smooth method better for precision of peak height by a factor of ' num2str(1/MeanAbsolutePositionRSDRatio)])
end
if MeanAbsoluteHeightRSDRatio>1
    disp(['Wavelet method better for precision of peak height by a factor of ' num2str(MeanAbsoluteHeightRSDRatio)])
else
    disp(['Segmented smooth method better for precision of peak height by a factor of ' num2str(1/MeanAbsoluteHeightRSDRatio)])
end

disp('------------------------------------------------------------')