% Demonstration of effect of wavelet level on signal consisting of three
% noisy unit-height Gaussian or Lorentzaian peaks. Requires gaussian.m or
% lorentzian.m, wdenoise.m, and plotit in the path. Wavelet type is set in
% lines 37, 38. Wavelet level varies from 1 to 11.

format compact
clear
% User-adjustable parameters
increment=.1;
x=5:increment:300;
% Peak parameters
PeakPosition1=75;
PeakWidth1=3; % Try widths from 10 to 50.
PeakPosition2=126;
PeakWidth2=10; % Try widths from 10 to 50.
PeakPosition3=200;
PeakWidth3=30; % Try widths from 10 to 50.
noise=.1; % Lower noise is better for wavelet
OriginalSignalToNoiseRatio=1/noise;
% Restrict measurements to the central region to avoid edge effects.
 plotrange=round(.2*length(x)):round(.8*length(x));
    
for trial=1:40
    level=2+round(trial/4)-1;
    Level(trial)=level;
    % Create a signal with a single Gaussian peak
    z=gaussian(x,PeakPosition1,PeakWidth1)+gaussian(x,PeakPosition2,PeakWidth2)+gaussian(x,PeakPosition3,PeakWidth3); % Alternative shape
    % z=lorentzian(x,PeakPosition,PeakWidth); % Alternative shape
    
    % Add random noise
    y=z+noise*randn(size(x));
    NoiseOnly=noise*randn(size(x));
    STDnoiseonly=std(NoiseOnly);
    SNR=1/STDnoiseonly;
    
    wdenoiseD=wdenoise(y, Level(trial),'Wavelet','coif5'); %  Specify wavelet type here
    DenoisedNoiseOnly=wdenoise(NoiseOnly, Level(trial),'Wavelet','coif5');
    STDDenoisedNoiseOnly=std(DenoisedNoiseOnly);
    SNRDenoised=1/STDDenoisedNoiseOnly;
    
    figure(1)
    clf
    subplot(2,1,1)
    plot(x(plotrange),y(plotrange),'r',x(plotrange),z(plotrange),'b')
    title('Original raw signal')
    xlabel(['Original signal-to-noise ratio = ' num2str(SNR)])
    ylabel('Signal amplitude')
    legend('Observed noisy signal','Original peak')
    axis([min(x) max(x) -0.2 1.2])
    subplot(2,1,2)
    plot(x(plotrange),z(plotrange),x(plotrange),wdenoiseD(plotrange))
    title('Wavelet denoise')
    axis([min(x) max(x) -0.1 1.1])
    xlabel(['Wavelet level = ' num2str(level) '    Signal-to-noise ratio = ' num2str(SNRDenoised)])
    legend('Original peak','Denoised raw signal')
    
    ResidualSTD(trial)=std(z(plotrange)-wdenoiseD(plotrange));
    SNRImprovementFactor(trial)=SNRDenoised/SNR;
    PeakHeight(trial)=max(wdenoiseD);
  
end
% Plot signal quality measurements vs wavelet level
figure(2)
clf
plotit(Level,PeakHeight,2,'o','c');
title('Height of denoised peak (true value = 1.0)')
xlabel('wavelet level')
grid 

figure(3)
clf
plotit(Level,ResidualSTD,2,'o','c');
title('Closeness to underlying noisless peak')
ylabel('Standard Deviation of Residual')
xlabel('wavelet level')
grid

figure(4)
clf
plotit(Level,sqrt(SNRImprovementFactor),1,'o','k');
title('Signal-to-noise ratio improvement')
xlabel('wavelet level')
ylabel('Square root of signal-to-noise ratio improvement factor')
grid