function dataout=intuitive_lowpass(datain,minwavelength,varargin)
% Function by Mike Wolovick, updated on 10/18/2011
% updated again on 2/7/2023

% This function lowpass filters an input in 1D by convolving it with a 
% gaussian weighting function.  

% The standard deviation of the gaussian is one-half the minwavelength. ie,
% full width at half max (FWHM) is roughly minwavelength.

% Minwavelength must be in terms of number of samples of datain. 

% A linear best-fit is removed before the data are filtered, then added
% back in at the end. 

% Edge effects are minimized by adding a buffer to either side of the input
% equal to the value at the endpoints.  This is done after the best-fit is
% removed.  IE, the boundary condition is an extrapolation of the endpoint
% values based on the best-fit slope of the entire input. 

% The function assumes that the input is evenly spaced.

% MOD on 10/18/2011:
% The function now offers an option to control how BC's are dealt with.
% The default is the old method described above, the alternate is an
% extrapolation of the endpoints based on the local slope.  A global
% bestfit is removed from the data in either case; this only affects the
% extrapolation method.

% MOD on 2/7/2023:  
% Adding the 'quasilocal' BC option.  Instead of extrapolating based on
% the slope at the endpoints, this option extrapolates based on a linear
% best-fit near the end points.  The best-fit is computed over a region
% comprising one wavelength from the boundary.  
% In addition, I have added informative comments throughout the code
% because apparently I was shit at commenting my code early in grad school.

% Syntax:
% output=intuitive_lowpass(input,minwavelength);
% or
% output=intuitive_lowpass(input,minwavelength,'method');  'method' can be
% either 'global' (default) or 'local'

%% Check Inputs:

% Check size and orientation of datain:
rowvec=0;
if isscalar(datain)==1
    error('Data Input must be a vector or a matrix.')
elseif size(minwavelength,1)~=1 || size(minwavelength,2)~=1
    error('Minimum Wavelength must be a scalar.')
elseif length(varargin)>1
    error('Function cannot have more than 3 inputs.')
elseif size(datain,1)==1
    rowvec=1;
    datain=datain';
end

% Check BC style:
if isempty(varargin)==0
    if strcmp(varargin,'local')==0 && strcmp(varargin,'quasilocal')==0 && strcmp(varargin,'global')==0
        error('method must be a string equal to "global", "local", or "quasilocal".')
    elseif strcmp(varargin,'local')==1
        bcmethod='local';
    elseif strcmp(varargin,'quasilocal')==1
        bcmethod='quasilocal';
    else
        bcmethod='global';
    end
else
    bcmethod='global';
end

%% Work:

% Compute length of input data:
L=length(datain);

% Compute integer wavelength:
integer_wavelength=max(1,round(minwavelength));

% Compute length of expanded dataset:
Lbig=L+6*integer_wavelength;

% Compute and remove linear best-fit:
A=[ones(L,1),linspace(1,L,L)'];
equation=(A'*A)\A'*datain;
bestline=equation(1)+equation(2)*linspace(1,L,L)';
datain=datain-bestline;

% Create expanded data vector:
xs=linspace(1,Lbig,Lbig)';
dataspread=zeros(Lbig,1);
dataspread(3*integer_wavelength+1:end-3*integer_wavelength)=datain;

% Check BC method:
if strcmp(bcmethod,'global')==1
    % Extrapolate endpoints based on overall best-fit slope:
    dataspread(1:3*integer_wavelength)=datain(1);
    dataspread(end-3*integer_wavelength+1:end)=datain(end);
elseif strcmp(bcmethod,'local')==1
    % Extrapolate based on local slope at endpoints:
    d_start=datain(1)-datain(2);
    dataspread(1:3*integer_wavelength)=datain(1)+d_start*linspace(3*integer_wavelength,1,3*integer_wavelength)';
    d_end=datain(end)-datain(end-1);
    dataspread(end-3*integer_wavelength+1:end)=datain(end)+d_end*linspace(1,3*integer_wavelength,3*integer_wavelength)';
else % quasilocal
    % Compute number of samples for the local fit:
    numendpointsamples=max(2,integer_wavelength);
    % Compute first local best-fit:
    A=[ones(numendpointsamples,1),xs(3*integer_wavelength+[1:numendpointsamples])];
    equation1=(A'*A)\A'*dataspread(3*integer_wavelength+[1:numendpointsamples]);
    % Fill values with the first local fit:
    dataspread(1:3*integer_wavelength)=equation1(1)+equation1(2)*xs(1:3*integer_wavelength);
    % Compute second local best-fit:
    A=[ones(numendpointsamples,1),xs(end-3*integer_wavelength+[-numendpointsamples+1:0])];
    equation2=(A'*A)\A'*dataspread(end-3*integer_wavelength+[-numendpointsamples+1:0]);
    % Fill values with the second local fit:
    dataspread(end-3*integer_wavelength+1:end)=equation2(1)+equation2(2)*xs(end-3*integer_wavelength+1:end);
end

% Compute filter function:
xcent=Lbig/2;
stdev=.5*minwavelength;
filter=exp(-.5*((xs-xcent)/stdev).^2);
filter=filter/sum(filter);

% Compute FFT of filter and data:
filtertrans=fft(filter);
datatrans=fft(dataspread);

% Perform convolution by multiplying in the frequency domain:
dataout=fftshift(ifft(filtertrans.*datatrans));

% Crop to central data region:
dataout=dataout(3*integer_wavelength:end-3*integer_wavelength-1);

% Add the linear best-fit back in:
dataout=dataout+bestline;

% Adjust to final orientation:
if rowvec==1
    dataout=dataout';
end

