% Script that compares various types of smoothing to wavelet denoise, for a signal
% consisting of several Gaussian peaks whose widths increase gradually over
% the duration of the signal. parameters (line 25-36) are chosen to
% optimize each method. Prints out the "residual", whr RSD of the
% difference between the true noisless signal and the denoised noisy
% signal, and the percent change in the peak height, width, and area of
% the first and last peaks. Requires fastsmooth.m, 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 and db2.

format short g
format compact
figure(1)
clf
% load data file consisting of several peaks of increasing width
load WideningPeaks.mat % load one data set created by DemoSegmentedSmooth.m
% This file contains noisless signal z

% Initial settings (user settable)
noise=.001;
d2=derivxy(x,derivxy(x,z)); % Noiseless second derivative of z
nz=z+noise.*randn(1,length(x)); % z signal with noise added
y=derivxy(x,derivxy(x,nz)); % Second derivatievve of noisy z
z=d2; % Reassign z

% Initial settings (user settable)
NumSegments=20;
ends=0;
SGwidth=50; % Savitsky-Golay smooth width
hamwidth=32; % width of Hamming filter in number of points
tau=400;% time constant (RC) of exponential filter
Level=4; % Wavelet level
% Define segmented smooth width vector based on peak width range
startw=11;
endw=50;
wstep=(endw-startw)/NumSegments;
smoothwidths=startw:wstep:endw;
gausswidth=50; % Gaussian weighted sliding average smooth width

clf
subplot(2,4,1)
plot(x,y,'r',x,z,'b')  % Graph the original signal in blue

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('SmoothVsWaveletsDerivative2Segmented.m')
legend('Noiseless derivative','Noisy derivative')

figure(1);subplot(2,4,2)
SmoothType=1; % 1=rectangular boxcar sliding average, 3=Gaussian weighted
tic;sy=fastsmooth(y,gausswidth,SmoothType,ends);t1=toc;
plot(x,z,'b',x,sy,'r')  % Graph the smoothed signal in red
thisaxis=axis; % Use the same plotting axes as before.
RSDresidual=100.*std((z-sy)./max(z));
RSDresidual=round(RSDresidual*100)./100;
xlabel(['RSD residual =  ' num2str(RSDresidual) '%' ])
title(['Fixed Gaussian smooth, average width =  ' num2str(mean(smoothwidths)) ])
legend('Noiseless data','Denoised noisy data')

figure(1);subplot(2,4,3)
SmoothType=3; % 1=rectangular boxcar sliding average, 3=Gaussian weighted
tic;sy=SegmentedSmooth(y,smoothwidths,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(['Segmented Gaussian smooth, average width =  ' num2str(mean(gausswidth)) ])
legend('Noiseless data','Denoised noisy data')

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')

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')

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')

figure(1);subplot(2,4,7)
tic;dny=wdenoise(y,Level,'Wavelet','coif2','ThresholdRule','soft');t2=toc;
plot(x,z,'b',x,dny,'r')   % Graph the denoised signal 
axis(thisaxis);
title(['coif2 soft 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')

figure(1);subplot(2,4,8)
tic;dny=wdenoise(y,Level,'Wavelet','bior5.5','ThresholdRule','soft');t2=toc;plot(x,dny,'r')
plot(x,z,'b',x,dny,'r')  % Graph the denoised signal 
axis(thisaxis);
title(['bior5.5 soft 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')

disp(' ')