% Script that compares various types of smoothing to wavelet denoise, for a signal
% consisting of two Gaussian peaks with very different widths. The smooth
% parameters (line 18-23) are chosen to optimize each method. Prints out
% the factor by which the wavelet method is slower and the change in the
% peak properties of the first and last peaks. Requires SegmentedSmooth.m,
% ProcessSignal.m, expfilter.m, hamfilter.m, halfwidth.m., peakfit.m, and
% wdenoise (from the Matlab wavelet toolbox). Compares unweighted sliding
% average, Gaussian weighted, Savitsky-Golay, exponential (RC) filter,
% Hamming filter, and two different wavelets.
% coif2 wavelet.
format short e
format compact
figure(1)
clf
% Create new data file consisting of several peaks of increasing width
x=1:10000; % load one data set created by DemoSegmentedSmooth.m
z=gaussian(x,2000,20)+gaussian(x,8000,1000);
y=z+0.1.*randn(size(x)); % Add extra noise for testing
% Optionally load previous data set
% load TwoGaussians

% Initial settings (user settable)
smoothwidths=10; % unweighted sliding average smooth width
gausswidth=15 ; % Gaussian weighted sliding average smooth width
segwidth=[5 250];
SmoothType=1; % 1=rectangular boxcar sliding average, 3=Gaussian weighted
ends=0;
SGwidth=12; % Savitsky-Golay smooth width
hamwidth=16; % width of Hamming filter in number of points
tau=80;% time constant (RC) of exponential filter
Level=9; % Wavelet level
Peak1Index=2000; % Index of narrowest (first) peak
Peak2Index=8000; % Index of widest (last) peak
w1=1950:2015; % Peak 1 window range
w2=6000:10000; % Peak 2 window range
% Peak properties of first and last peaks (for comparison to denoised data
TrueHeight1=max(z(w1));
TrueHeight2=max(z(w2));
TrueHalfwidth1=halfwidth(x,z,Peak1Index);
TrueHalfwidth2=halfwidth(x,z,Peak2Index);
TrueArea1=trapz(x(w1),z(w1));
TrueArea2=trapz(x(w2),z(w2));
disp('                               Percent errors of peak1|peak2')
disp('Method       Residuals      Height         Width          Area ')

clf
subplot(2,4,1)
plot(x,y,'b',x,z,'r')  % Graph the original signal in blue
thisaxis=axis;
RSDresidual=100.*std((z-y)./max(z));
RSDresidual=round(RSDresidual*100)./100;
xlabel(['Original unsmoothed data. RSD residual =  ' num2str(100.*std((z-y)./max(z))) '%' ])
title('SmoothVsWavelets2Gaussians.m')
legend('Noisy data','Noiseless data')
% Measure percent errors of peak properties of first and last peaks
 Measured1=max(y(w1));
 PercentHeightError1=-100.*(TrueHeight1-Measured1)/TrueHeight1;
 PercentHeightError1=round(PercentHeightError1*100)./100;
 Measured2=max(y(w2));
 PercentHeightError2=-100.*(TrueHeight2-Measured2)/TrueHeight2;
 PercentHeightError2=round(PercentHeightError2*100)./100;
 Measured1=halfwidth(x,y,Peak1Index);
 PercentWidthError1=-100.*(TrueHalfwidth1-Measured1)/TrueHalfwidth1;
 PercentWidthError1=round(PercentWidthError1*100)./100;
 Measured2=halfwidth(x,y,Peak2Index);
 PercentWidthError2=-100.*(TrueHalfwidth2-Measured2)/TrueHalfwidth2;
 PercentWidthError2=round(PercentWidthError2*100)./100;
 Measured1=trapz(x(w1),y(w1));
 PercentAreaError1=-100.*(TrueArea1-Measured1)/TrueArea1;
 PercentAreaError1=round(PercentAreaError1*100)./100;
 Measured2=trapz(x(w2),y(w2));
 PercentAreaError2=-100.*(TrueArea2-Measured2)/TrueArea2;
 PercentAreaError2=round(PercentAreaError2*100)./100;
disp(['Original       ' num2str(RSDresidual) '%   ' [num2str(PercentHeightError1) '%|' num2str(PercentHeightError2)] '%   ' [num2str(PercentWidthError1) '%|' num2str(PercentWidthError2)] '%   ' [num2str(PercentAreaError1) '%|' num2str(PercentAreaError2) '%'] ])

figure(1);subplot(2,4,2)
SmoothType=3; % Gaussian smooth
tic;sy=SegmentedSmooth(y,segwidth,SmoothType,ends);t2=toc;
NumSegments=length(segwidth);
plot(x,z,'b',x,sy,'r')  % Graph the smoothed signal in red
axis(thisaxis); % Use the same plotting axes as before.
RSDresidual=100.*std((z-sy)./max(z));
RSDresidual=round(RSDresidual*100)./100;
xlabel(['RSD residual =  ' num2str(RSDresidual) '%' ])
SegGaussianSpeedPenalty=round(t2/t1); % How much slower this function is
title(['Segmented Gaussian weighted, widths =  ' num2str(segwidth) ])
legend('Noiseless data','Denoised noisy data')
% Measure percent errors of peak properties of first and last peaks
 Measured1=max(sy(w1));
 PercentHeightError1=-100.*(TrueHeight1-Measured1)/TrueHeight1;
 PercentHeightError1=round(PercentHeightError1*100)./100;
 Measured2=max(sy(w2));
 PercentHeightError2=-100.*(TrueHeight2-Measured2)/TrueHeight2;
 PercentHeightError2=round(PercentHeightError2*100)./100;
 Measured1=halfwidth(x,sy,Peak1Index);
 PercentWidthError1=-100.*(TrueHalfwidth1-Measured1)/TrueHalfwidth1;
 PercentWidthError1=round(PercentWidthError1*100)./100;
 Measured2=halfwidth(x,sy,Peak2Index);
 PercentWidthError2=-100.*(TrueHalfwidth2-Measured2)/TrueHalfwidth2;
 PercentWidthError2=round(PercentWidthError2*100)./100;
 Measured1=trapz(x(w1),sy(w1));
 PercentAreaError1=-100.*(TrueArea1-Measured1)/TrueArea1;
 PercentAreaError1=round(PercentAreaError1*100)./100;
 Measured2=trapz(x(w2),sy(w2));
 PercentAreaError2=-100.*(TrueArea2-Measured2)/TrueArea2;
 PercentAreaError2=round(PercentAreaError2*100)./100;
disp(['Gaussian        '  num2str(RSDresidual) '%   ' [num2str(PercentHeightError1) '%|' num2str(PercentHeightError2)] '%   ' [num2str(PercentWidthError1) '%|' num2str(PercentWidthError2)] '%   ' [num2str(PercentAreaError1) '%|' num2str(PercentAreaError2) '%'] ])

figure(1);subplot(2,4,3)
SmoothType=3; % Gaussian smooth
tic;sy=SegmentedSmooth(y,gausswidth,SmoothType,ends);t2=toc;
NumSegments=length(smoothwidths);
plot(x,z,'b',x,sy,'r')  % Graph the smoothed signal in red
axis(thisaxis); % Use the same plotting axes as before.
RSDresidual=100.*std((z-sy)./max(z));
RSDresidual=round(RSDresidual*100)./100;
xlabel(['RSD residual =  ' num2str(RSDresidual) '%' ])
GaussianSpeedPenalty=round(t2/t1); % How much slower this function is
title(['Gaussian weighted segmented, width =  ' num2str(gausswidth) ])
legend('Noiseless data','Denoised noisy data')
% Measure percent errors of peak properties of first and last peaks
 Measured1=max(sy(w1));
 PercentHeightError1=-100.*(TrueHeight1-Measured1)/TrueHeight1;
 PercentHeightError1=round(PercentHeightError1*100)./100;
 Measured2=max(sy(w2));
 PercentHeightError2=-100.*(TrueHeight2-Measured2)/TrueHeight2;
 PercentHeightError2=round(PercentHeightError2*100)./100;
 Measured1=halfwidth(x,sy,Peak1Index);
 PercentWidthError1=-100.*(TrueHalfwidth1-Measured1)/TrueHalfwidth1;
 PercentWidthError1=round(PercentWidthError1*100)./100;
 Measured2=halfwidth(x,sy,Peak2Index);
 PercentWidthError2=-100.*(TrueHalfwidth2-Measured2)/TrueHalfwidth2;
 PercentWidthError2=round(PercentWidthError2*100)./100;
 Measured1=trapz(x(w1),sy(w1));
 PercentAreaError1=-100.*(TrueArea1-Measured1)/TrueArea1;
 PercentAreaError1=round(PercentAreaError1*100)./100;
 Measured2=trapz(x(w2),sy(w2));
 PercentAreaError2=-100.*(TrueArea2-Measured2)/TrueArea2;
 PercentAreaError2=round(PercentAreaError2*100)./100;
disp(['Segmented        '  num2str(RSDresidual) '%   ' [num2str(PercentHeightError1) '%|' num2str(PercentHeightError2)] '%   ' [num2str(PercentWidthError1) '%|' num2str(PercentWidthError2)] '%   ' [num2str(PercentAreaError1) '%|' num2str(PercentAreaError2) '%'] ])

figure(1);subplot(2,4,4)
tic;dny=ProcessSignal(x,y,0,SGwidth,4,0,0,0,0,0,0,0,0);t2=toc;
plot(x,z,'b',x,dny,'r')   % Graph the denoised signal 
axis(thisaxis);
title(['Savitsky-Golay smooth, width = ' num2str(SGwidth) ])
RSDresidual=100.*std((z-dny)./max(z));
RSDresidual=round(RSDresidual*100)./100;
xlabel(['RSD residual =  ' num2str(RSDresidual) '%' ])
SavitskyGolaySpeedPenalty=round(t2/t1); % How much slower this function is
legend('Noiseless data','Denoised noisy data')
% Measure percent errors of peak properties of first and last peaks
 Measured1=max(dny(w1));
 PercentHeightError1=-100.*(TrueHeight1-Measured1)/TrueHeight1;
 PercentHeightError1=round(PercentHeightError1*100)./100;
 Measured2=max(dny(w2));
 PercentHeightError2=-100.*(TrueHeight2-Measured2)/TrueHeight2;
 PercentHeightError2=round(PercentHeightError2*100)./100;
 Measured1=halfwidth(x,dny,Peak1Index);
 PercentWidthError1=-100.*(TrueHalfwidth1-Measured1)/TrueHalfwidth1;
 PercentWidthError1=round(PercentWidthError1*100)./100;
 Measured2=halfwidth(x,dny,Peak2Index);
 PercentWidthError2=-100.*(TrueHalfwidth2-Measured2)/TrueHalfwidth2;
 PercentWidthError2=round(PercentWidthError2*100)./100;
 Measured1=trapz(x(w1),dny(w1));
 PercentAreaError1=-100.*(TrueArea1-Measured1)/TrueArea1;
 PercentAreaError1=round(PercentAreaError1*100)./100;
 Measured2=trapz(x(w2),dny(w2));
 PercentAreaError2=-100.*(TrueArea2-Measured2)/TrueArea2;
 PercentAreaError2=round(PercentAreaError2*100)./100;
disp(['Savitsky-Golay  '  num2str(RSDresidual) '%   ' [num2str(PercentHeightError1) '%|' num2str(PercentHeightError2)] '%   ' [num2str(PercentWidthError1) '%|' num2str(PercentWidthError2)] '%   ' [num2str(PercentAreaError1) '%|' num2str(PercentAreaError2) '%'] ])

figure(1);subplot(2,4,5)
tic;RCf=expfilter(x,y,tau);RCftoc=toc;
plot(x,z,'b',x,RCf,'r')   % Graph the Fourier filter signal 
axis(thisaxis);
title('RC filter (low pass)')
title(['RC filter tau =  ' num2str(tau) ])
RSDresidual=100.*std((z-RCf)./max(z));
RSDresidual=round(RSDresidual*100)./100;
xlabel(['RSD residual =  ' num2str(RSDresidual) '%' ])
expFilterSpeedPenalty=round(RCftoc/t1); % How much slower this function is
legend('Noiseless data','Denoised noisy data')
% Measure percent errors of peak properties of first and last peaks
 Measured1=max(RCf(w1));
 PercentHeightError1=-100.*(TrueHeight1-Measured1)/TrueHeight1;
 PercentHeightError1=round(PercentHeightError1*100)./100;
 Measured2=max(RCf(w2));
 PercentHeightError2=-100.*(TrueHeight2-Measured2)/TrueHeight2;
 PercentHeightError2=round(PercentHeightError2*100)./100;
 Measured1=halfwidth(x,RCf,Peak1Index);
 PercentWidthError1=-100.*(TrueHalfwidth1-Measured1)/TrueHalfwidth1;
 PercentWidthError1=round(PercentWidthError1*100)./100;
 Measured2=halfwidth(x,RCf,Peak2Index);
 PercentWidthError2=-100.*(TrueHalfwidth2-Measured2)/TrueHalfwidth2;
 PercentWidthError2=round(PercentWidthError2*100)./100;
 Measured1=trapz(x(w1),RCf(w1));
 PercentAreaError1=-100.*(TrueArea1-Measured1)/TrueArea1;
 PercentAreaError1=round(PercentAreaError1*100)./100;
 Measured2=trapz(x(w2),RCf(w2));
 PercentAreaError2=-100.*(TrueArea2-Measured2)/TrueArea2;
 PercentAreaError2=round(PercentAreaError2*100)./100;
disp(['RC filter       '  num2str(RSDresidual) '%   ' [num2str(PercentHeightError1) '%|' num2str(PercentHeightError2)] '%   ' [num2str(PercentWidthError1) '%|' num2str(PercentWidthError2)] '%   ' [num2str(PercentAreaError1) '%|' num2str(PercentAreaError2) '%'] ])

figure(1);subplot(2,4,6)
tic;hamf=hamfilter(y,hamwidth);t2=toc;
plot(x,z,'b',x,hamf,'r')   % Graph the denoised signal 
axis(thisaxis);
title(['Hamming filter, width = ' num2str(hamwidth) ])
RSDresidual=100.*std((z-hamf)./max(z));
RSDresidual=round(RSDresidual*100)./100;
xlabel(['RSD residual =  ' num2str(RSDresidual) '%' ])
HammingSpeedPenalty=round(t2/t1); % How much slower this function is
legend('Noiseless data','Denoised noisy data')
% Measure percent errors of peak properties of first and last peaks
 Measured1=max(hamf(w1));
 PercentHeightError1=-100.*(TrueHeight1-Measured1)/TrueHeight1;
 PercentHeightError1=round(PercentHeightError1*100)./100;
 Measured2=max(hamf(w2));
 PercentHeightError2=-100.*(TrueHeight2-Measured2)/TrueHeight2;
 PercentHeightError2=round(PercentHeightError2*100)./100;
 Measured1=halfwidth(x,hamf,Peak1Index);
 PercentWidthError1=-100.*(TrueHalfwidth1-Measured1)/TrueHalfwidth1;
 PercentWidthError1=round(PercentWidthError1*100)./100;
 Measured2=halfwidth(x,hamf,Peak2Index);
 PercentWidthError2=-100.*(TrueHalfwidth2-Measured2)/TrueHalfwidth2;
 PercentWidthError2=round(PercentWidthError2*100)./100;
 Measured1=trapz(x(w1),hamf(w1));
 PercentAreaError1=-100.*(TrueArea1-Measured1)/TrueArea1;
 PercentAreaError1=round(PercentAreaError1*100)./100;
 Measured2=trapz(x(w2),hamf(w2));
 PercentAreaError2=-100.*(TrueArea2-Measured2)/TrueArea2;
 PercentAreaError2=round(PercentAreaError2*100)./100;
disp(['Hamming filter  '  num2str(RSDresidual) '%   ' [num2str(PercentHeightError1) '%|' num2str(PercentHeightError2)] '%   ' [num2str(PercentWidthError1) '%|' num2str(PercentWidthError2)] '%   ' [num2str(PercentAreaError1) '%|' num2str(PercentAreaError2) '%'] ])

figure(1);subplot(2,4,7)
tic;dny=wdenoise(y,Level,'Wavelet','coif2');t2=toc;
plot(x,z,'b',x,dny,'r')   % Graph the denoised signal 
axis(thisaxis);
title(['coif2 wavelet denoise, Level = ' num2str(Level) ])
RSDresidual=100.*std((z-dny)./max(z));
RSDresidual=round(RSDresidual*100)./100;
xlabel(['RSD residual =  ' num2str(RSDresidual) '%' ])
coif2WaveletSpeedPenalty=round(t2/t1); % How much slower this function is
legend('Noiseless data','Denoised noisy data')
% Measure percent errors of peak properties of first and last peaks
 Measured1=max(dny(w1));
 PercentHeightError1=-100.*(TrueHeight1-Measured1)/TrueHeight1;
 PercentHeightError1=round(PercentHeightError1*100)./100;
 Measured2=max(dny(w2));
 PercentHeightError2=-100.*(TrueHeight2-Measured2)/TrueHeight2;
 PercentHeightError2=round(PercentHeightError2*100)./100;
 Measured1=halfwidth(x,dny,Peak1Index);
 PercentWidthError1=-100.*(TrueHalfwidth1-Measured1)/TrueHalfwidth1;
 PercentWidthError1=round(PercentWidthError1*100)./100;
 Measured2=halfwidth(x,dny,Peak2Index);
 PercentWidthError2=-100.*(TrueHalfwidth2-Measured2)/TrueHalfwidth2;
 PercentWidthError2=round(PercentWidthError2*100)./100;
 Measured1=trapz(x(w1),dny(w1));
 PercentAreaError1=-100.*(TrueArea1-Measured1)/TrueArea1;
 PercentAreaError1=round(PercentAreaError1*100)./100;
 Measured2=trapz(x(w2),dny(w2));
 PercentAreaError2=-100.*(TrueArea2-Measured2)/TrueArea2;
 PercentAreaError2=round(PercentAreaError2*100)./100;
disp(['coif2 wavelet   '  num2str(RSDresidual) '%   ' [num2str(PercentHeightError1) '%|' num2str(PercentHeightError2)] '%   ' [num2str(PercentWidthError1) '%|' num2str(PercentWidthError2)] '%   ' [num2str(PercentAreaError1) '%|' num2str(PercentAreaError2) '%'] ])

figure(1);subplot(2,4,8)
tic;dny=wdenoise(y,Level,'Wavelet','db2');t2=toc;plot(x,dny,'r')
plot(x,z,'b',x,dny,'r')   % Graph the denoised signal 
axis(thisaxis);
title(['db2 wavelet denoise, Level = ' num2str(Level) ])
RSDresidual=100.*std((z-dny)./max(z));
RSDresidual=round(RSDresidual*100)./100;
xlabel(['RSD residual =  ' num2str(RSDresidual) '%' ])
db2WaveletSpeedPenalty=round(t2/t1);
legend('Noiseless data','Denoised noisy data')
% Measure percent errors of peak properties of first and last peaks
 Measured1=max(dny(w1));
 PercentHeightError1=-100.*(TrueHeight1-Measured1)/TrueHeight1;
 PercentHeightError1=round(PercentHeightError1*100)./100;
 Measured2=max(dny(w2));
 PercentHeightError2=-100.*(TrueHeight2-Measured2)/TrueHeight2;
 PercentHeightError2=round(PercentHeightError2*100)./100;
 Measured1=halfwidth(x,dny,Peak1Index);
 PercentWidthError1=-100.*(TrueHalfwidth1-Measured1)/TrueHalfwidth1;
 PercentWidthError1=round(PercentWidthError1*100)./100;
 Measured2=halfwidth(x,dny,Peak2Index);
 PercentWidthError2=-100.*(TrueHalfwidth2-Measured2)/TrueHalfwidth2;
 PercentWidthError2=round(PercentWidthError2*100)./100;
 Measured1=trapz(x(w1),dny(w1));
 PercentAreaError1=-100.*(TrueArea1-Measured1)/TrueArea1;
 PercentAreaError1=round(PercentAreaError1*100)./100;
 Measured2=trapz(x(w2),dny(w2));
 PercentAreaError2=-100.*(TrueArea2-Measured2)/TrueArea2;
 PercentAreaError2=round(PercentAreaError2*100)./100;
disp(['db2 wavelet     '  num2str(RSDresidual) '%   ' [num2str(PercentHeightError1) '%|' num2str(PercentHeightError2)] '%   ' [num2str(PercentWidthError1) '%|' num2str(PercentWidthError2)] '%   ' [num2str(PercentAreaError1) '%|' num2str(PercentAreaError2) '%'] ])

disp(' ')