%% Data Smoothing
% with median filter, recursive sliding average, Savitsky-Golay, Fourier low-pass, or wavelet denoise (which requires the Matlab  Wavelet Toolbox). To see the graphs on the right, right-click on the right panel and select "Disable synchronous scrolling".
% Tom O'Haver 4/18/2023

 file= uigetfile('*.csv;*.xlsx');% Click this button to load data
mydata=xlsread(file); disp(file) %                       from disk in xlsx or csv format. 

%

PlotBeforeAndAfter= true; % Check this box to display smoothed signal as well as original.
FrequencySpectra=true; % Check this box to display frequency spectra

% Set the x-axis scale expansion beginning and end points (0 - 100%)
startpc=1; % Percentage of data points to start data selection
endpc=97.6; % Percentage of data points to end data selection

RemoveSpikes=false; % Check this box to use median filter to remove spikes
SpikeWidth=6; % Estimated spike width, in data points

SmoothType ="Sliding average"; % Select smoothing algrorithm

% For sliding average:
SmoothWidth=9; % Smooth width for sliding average
NumPasses=4; % Number of sliding average passes (1-5)
SmoothEnds= true; % Click this box to smooth the ends of the data record.

% For Fourier low-pass filter:
FrequencyCutoff=3.5; % Fourier filter cutoff frequency
CutOffRate=6; % Fourier filter cutoff rate

% For Savitsky-Golay smooth only
SGSmoothWidth=63; % Frame length for Savitsky-Golay
PolynomialOrder=5; % Polynomial order must be less than the frame length

% For wavelet denoising only (requires Matlab Wavelet Toolox be installed)
WaveletType="bior5.5";
Level=4; % The higher the level, the lower is the frequency divider between signal and noise
% -------------------------------------------------------------------------------
NumPoints=length(mydata);
startpoint=round(startpc*NumPoints./100);
endpoint=round(endpc*NumPoints./100);
x=mydata(startpoint:endpoint,1);
rawy=mydata(startpoint:endpoint,2);

if RemoveSpikes==true
    y=medianfilter(rawy,SpikeWidth);
else
    y=rawy;
end

if SmoothType=="Sliding average"
    % Recursive sliding average smoothing
    sy=fastsmooth(y,SmoothWidth,NumPasses,SmoothEnds);
end

if SmoothType=="Fourier low-pass filter"
    % Percentage of entire spectrum included for Fourier filter
    sy=FourierFilter(x,y,0,FrequencyCutoff.*0.01.*length(x),CutOffRate,1);
end

if SmoothType=="Wavelet Denoise"
    if WaveletType=="BlockJS"
        sy=wdenoise(y,Level,'DenoisingMethod','BlockJS');
    end
    if WaveletType=="bior5.5"
        sy=wdenoise(y,Level,'Wavelet','bior5.5','ThresholdRule','soft');
    end
    if WaveletType=="coif2"
        sy=wdenoise(y,Level,'Wavelet','coif2','ThresholdRule','soft');
    end
    if WaveletType=="sym8"
        sy=wdenoise(y,Level,'Wavelet','sym8','ThresholdRule','soft');
    end
    if WaveletType=="db4"
        sy=wdenoise(y,Level,'Wavelet','db4','ThresholdRule','soft');
    end
end

if SmoothType=="Savitsky-Golay"
    sy=savitzkyGolayFilt(y,PolynomialOrder,0,SGSmoothWidth);
end

figure(1)
%clf
if PlotBeforeAndAfter==true
    plot(x,y,'k',x,sy,'r'); title(SmoothType);grid
    legend('Original data','Smoothed data')
else
    plot(x,sy,'r'); title(SmoothType);grid
    legend('Smoothed data')
end

AreaRecovery=100.*trapz(x,sy)/trapz(x,y);
xlabel(['Area recovery: ' num2str(AreaRecovery) ' %'])

% Optional Fourier spectrum plotted if FrequencySpectra checkbox is checked
if FrequencySpectra==true
    figure(2)
    %clf
    fyc=fft(sy);
    syc=fyc .* conj(fyc);
    minspec=min(syc);maxspec=max(syc);
    psrange=1:length(fyc)/2;
    psyc=real(syc(psrange));
    loglog(psrange,psyc,'r','linewidth',1)
    grid
    axis([1 NumPoints/4 minspec maxspec])
    title('Frequency spectra')
    xlabel('Frequency')
    legend('Smoothed data',"Location",'southwest')
    if PlotBeforeAndAfter==true
        fyc=fft(y);
        syc=fyc .* conj(fyc);
        minspec=min(syc);maxspec=max(syc);
        psrange=1:length(fyc)/2;
        psyc=real(syc(psrange));
        hold on
        loglog(psrange,psyc,'k','linewidth',1)
        hold off
        legend('Original data',' Smoothed data',"Location",'southwest')
    end
else
end

%

function ry=FourierFilter(xvector,yvector,centerfrequency,filterwidth,filtershape,filtermode)
% function ry=FourierFilter(xvector,yvector,centerfrequency,filterwidth,filtershape,filtermode)
% Computes fourier filter for signal yvector. Centerfrequency and
% filterwidth are the center frequency and width of the pass band, in
% harmonics, input arguments set the intital values of center frequency,
% filter width, shape, and filtermode (1-6):
% 1 = 'Band-pass'
% 2 = 'Lowpass'
% 3 = 'Highpass',
% 4 = 'Band-reject (notch)',
% 5 = 'Comb pass')
% 6 = 'Comb notch'
% 'filtershape' determines the sharpness of the cut-off.
%
% Adjust x and y vector shape to 1 x n (rather than n x 1)
xvector=reshape(xvector,1,length(xvector));
yvector=reshape(yvector,1,length(yvector));
fy=fft(yvector);
lft1=[1:(length(fy)/2)];
lft2=[(length(fy)/2+1):length(fy)];
% Compute filter shape.
if filtermode==1 % 'Band-pass'
    ffilter1=shapefunction(lft1,centerfrequency+1,filterwidth,filtershape);
    ffilter2=shapefunction(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
    ffilter=[ffilter1,ffilter2];
end
if filtermode==2 % 'High-pass'
    centerfrequency=length(xvector)/2;
    ffilter1=shapefunction(lft1,centerfrequency+1,filterwidth,filtershape);
    ffilter2=shapefunction(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
    ffilter=[ffilter1,ffilter2];
end
if filtermode==3 % 'Low-pass'
    centerfrequency=0;
    ffilter1=shapefunction(lft1,centerfrequency+1,filterwidth,filtershape);
    ffilter2=shapefunction(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
    ffilter=[ffilter1,ffilter2];
end
if filtermode==4 % 'Band-reject (notch)'
    ffilter1=shapefunction(lft1,centerfrequency+1,filterwidth,filtershape);
    ffilter2=shapefunction(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
    ffilter=1-[ffilter1,ffilter2];
end
if filtermode==5 % 'Comb pass'
    n=2;
    ffilter1=shapefunction(lft1,centerfrequency+1,filterwidth,filtershape);
    ffilter2=shapefunction(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
    while n<50
        ffilter1=ffilter1+shapefunction(lft1,n*(centerfrequency+1),filterwidth,filtershape);
        ffilter2=ffilter2+shapefunction(lft2,length(fy)-n*(centerfrequency+1),filterwidth,filtershape);
        n=n+1;
    end
    ffilter=[ffilter1,ffilter2];
end
if filtermode==6 % 'Comb notch'
    n=2;
    ffilter1=shapefunction(lft1,centerfrequency+1,filterwidth,filtershape);
    ffilter2=shapefunction(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
    while n<50
        ffilter1=ffilter1+shapefunction(lft1,n*(centerfrequency+1),filterwidth,filtershape);
        ffilter2=ffilter2+shapefunction(lft2,length(fy)-n*(centerfrequency+1),filterwidth,filtershape);
        n=n+1;
    end
    ffilter=1-[ffilter1,ffilter2];
end
if length(fy)>length(ffilter), ffilter=[ffilter ffilter(1)];end
ffy=fy.*ffilter;  % Multiply filter by Fourier Transform of signal
ry=real(ifft(ffy));

    function g = shapefunction(x,pos,wid,n)
        %  shapefunction(x,pos,wid,n) = peak centered on x=pos, half-width=wid
        %  x may be scalar, vector, or matrix, pos and wid both scalar
        %  Shape is Lorentzian (1/x^2) when n=0, Gaussian (exp(-x^2))
        %  when n=1, and becomes more rectangular as n increases.
        %  Example: shapefunction([1 2 3],1,2,1) gives result [1.0000    0.5000    0.0625]
        if n==0
            g=ones(size(x))./(1+((x-pos)./(0.5.*wid)).^2);
        else
            g = exp(-((x-pos)./(0.6.*wid)) .^(2*round(n)));
        end
    end
end

function SmoothY=fastsmooth(Y,w,type,ends)
% fastsmooth(Y,w,type,ends) smooths vector Y with smooth
%  of width w. Version 3.0, October 2016.
% The argument "type" determines the smooth type:
%   If type=1, rectangular (sliding-average or boxcar)
%   If type=2, triangular (2 passes of sliding-average)
%   If type=3, pseudo-Gaussian (3 passes of sliding-average)
%   If type=4, pseudo-Gaussian (4 passes of same sliding-average)
%   If type=5, multiple-width (4 passes of different sliding-average)
% The argument "ends" controls how the "ends" of the signal
% (the first w/2 points and the last w/2 points) are handled.
%   If ends=0, the ends are zero.  (In this mode the elapsed
%     time is independent of the smooth width). The fastest.
%   If ends=1, the ends are smoothed with progressively
%     smaller smooths the closer to the end. (In this mode the
%     elapsed time increases with increasing smooth widths).
% fastsmooth(Y,w,type) smooths with ends=0.
% fastsmooth(Y,w) smooths with type=1 and ends=0.
%
% Examples:
% fastsmooth([1 1 1 10 10 10 1 1 1 1],3)= [0 1 4 7 10 7 4 1 1 0]
% fastsmooth([1 1 1 10 10 10 1 1 1 1],3,1,1)= [1 1 4 7 10 7 4 1 1 1]
%
% x=1:100;
% y=randn(size(x));
% plot(x,y,x,fastsmooth(y,5,3,1),'r')
% xlabel('Blue: white noise.    Red: smoothed white noise.')
%
% Demonstration scripts and functions: SmoothExperiment.m,
% SmoothWidthTest.m, smoothdemo.m, SmoothOptimization.m,
% DemoSegmentedSmooth.m, DeltaTest.m, SmoothVsCurvefit.m
%
% Related functions: SegmentedSmooth.m, medianfilter.m, killspikes.m
% Download from
% https://terpconnect.umd.edu/~toh/spectrum/functions.html

% Copyright (c) 2012,2018 Thomas C. O'Haver
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
% copies of the Software, and to permit persons to whom the Software is
% furnished to do so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
% THE SOFTWARE.
if nargin==2, ends=0; type=1; end
if nargin==3, ends=0; end
switch type
    case 1
        SmoothY=sa(Y,w,ends);
    case 2
        SmoothY=sa(sa(Y,w,ends),w,ends);
    case 3
        SmoothY=sa(sa(sa(Y,w,ends),w,ends),w,ends);
    case 4
        SmoothY=sa(sa(sa(sa(Y,w,ends),w,ends),w,ends),w,ends);
    case 5
        SmoothY=sa(sa(sa(sa(Y,round(1.6*w),ends),round(1.4*w),ends),round(1.2*w),ends),w,ends);
end
end

function SmoothY=sa(Y,smoothwidth,ends)
w=round(smoothwidth);
SumPoints=sum(Y(1:w));
s=zeros(size(Y));
halfw=round(w/2);
L=length(Y);
for k=1:L-w
    s(k+halfw-1)=SumPoints;
    SumPoints=SumPoints-Y(k);
    SumPoints=SumPoints+Y(k+w);
end
s(k+halfw)=sum(Y(L-w+1:L));
SmoothY=s./w;
% Taper the ends of the signal if ends=1.
if ends==1
    startpoint=(smoothwidth + 1)/2;
    SmoothY(1)=(Y(1)+Y(2))./2;
    for k=2:startpoint
        SmoothY(k)=mean(Y(1:(2*k-1)));
        SmoothY(L-k+1)=mean(Y(L-2*k+2:L));
    end
    SmoothY(L)=(Y(L)+Y(L-1))./2;
end
end

function y=savitzkyGolayFilt(x,N,DN,F,W,DIM)
%savitzkyGolayFilt Savitzky-Golay Filtering.
%   savitzkyGolayFilt(X,N,DN,F) filters the signal X using a Savitzky-Golay
%   (polynomial) filter.  The polynomial order, N, must be less than the
%   frame size, F, and F must be odd.  DN specifies the differentiation
%   order (DN=0 is smoothing). For a DN higher than zero, you'll have to
%   scale the output by 1/T^DN to acquire the DNth smoothed derivative of
%   input X, where T is the sampling interval. The length of the input X
%   must be >= F.  If X is a matrix, the filtering is done on the columns
%   of X.
%
%   Note that if the polynomial order N equals F-1, no smoothing
%   will occur.
%
%   savitzkyGolayFilt(X,N,DN,F,W) specifies a weighting vector W with
%   length F containing real, positive valued weights employed during the
%   least-squares minimization. If not specified, or if specified as
%   empty, W defaults to an identity matrix.
%
%   savitzkyGolayFilt(X,N,DN,F,[],DIM) or savitzkyGolayFilt(X,N,DN,F,W,DIM)
%   operates along the dimension DIM.
%
%   See also savitzkyGolay, FILTER, sgolayfilt

%   References:
%     [1] Sophocles J. Orfanidis, INTRODUCTION TO SIGNAL PROCESSING,
%              Prentice-Hall, 1995, Chapter 8.

%   Author(s): R. Losada
%   Copyright 1988-2004 The MathWorks, Inc.
%   $Revision: 1.11.4.4 $  $Date: 2009/08/11 15:47:54 $

error(nargchk(4,6,nargin,'struct'));

% Check if the input arguments are valid
if round(F) ~= F, error(generatemsgid('MustBeInteger'),'Frame length must be an integer.'), end
if rem(F,2) ~= 1, error(generatemsgid('SignalErr'),'Frame length must be odd.'), end
if round(N) ~= N, error(generatemsgid('MustBeInteger'),'Polynomial order must be an integer.'), end
if N > F-1, error(generatemsgid('InvalidRange'),'The Polynomial order must be less than the frame length.'), end
if DN > N, error(generatemsgid('InvalidRange'),'The Differentiation order must be less than or equal to the Polynomial order.'), end

if nargin < 5 || isempty(W)
    % No weighting matrix, make W an identity
    W = ones(F,1);
else
    % Check for right length of W
    if length(W) ~= F, error(generatemsgid('InvalidDimensions'),'The weight vector must be of the same length as the frame length.'),end
    % Check to see if all elements are positive
    if min(W) <= 0, error(generatemsgid('InvalidRange'),'All the elements of the weight vector must be greater than zero.'), end
end

if nargin < 6, DIM = []; end

% Compute the projection matrix B
pp = fix(-F./2):fix(F./2);
B = savitzkyGolay(pp,N,DN,pp,W);

if ~isempty(DIM) && DIM > ndims(x)
    error(generatemsgid('InvalidDimensions'),'Dimension specified exceeds the dimensions of X.')
end

% Reshape X into the right dimension.
if isempty(DIM)
    % Work along the first non-singleton dimension
    [x, nshifts] = shiftdim(x);
else
    % Put DIM in the first dimension (this matches the order
    % that the built-in filter function uses)
    perm = [DIM,1:DIM-1,DIM+1:ndims(x)];
    x = permute(x,perm);
end

if size(x,1) < F, error(generatemsgid('InvalidDimensions'),'The length of the input must be >= frame length.'), end

% Preallocate output
y = zeros(size(x));

% Compute the transient on (note, this is different than in sgolayfilt,
% they had an optimization leaving out some transposes that is only valid
% for DN==0)
y(1:(F+1)/2-1,:) = fliplr(B(:,(F-1)/2+2:end)).'*flipud(x(1:F,:));

% Compute the steady state output
ytemp = filter(B(:,(F-1)./2+1),1,x);
y((F+1)/2:end-(F+1)/2+1,:) = ytemp(F:end,:);

% Compute the transient off
y(end-(F+1)/2+2:end,:) = fliplr(B(:,1:(F-1)/2)).'*flipud(x(end-(F-1):end,:));

% Convert Y to the original shape of X
if isempty(DIM)
    y = shiftdim(y, -nshifts);
else
    y = ipermute(y,perm);
end
end