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

%

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

% Filter type: Low, high, band, notch, or comb.
filtermode="Band-pass"; 
% Center frequency for band pass and notch filters, % of max
CenterFrequency=14;
% Filter width 
filterwidth=12; 
% Filter Shape: 1=Gaussian; higher => more rectangular
filtershape=7; 
% Set linear or log scales on x and y axes
plotmode=1; 
% Choose between frequency and period x scale
xmode=0; 
% Fixed or variable y scale on output plot
ymode=0; 
% Optionally play the output wavefrom as a sound.
PlayOutputAsSound=false;

NumPoints=length(data);
startpoint=round(startpc*NumPoints./100);
endpoint=round(endpc*NumPoints./100);
xvector=data(startpoint:endpoint,1);
yvector=data(startpoint:endpoint,2);
samplingtime=max(xvector)-min(xvector);
centerfrequency=(length(xvector)/2)*CenterFrequency/200;
filterWidth=(length(xvector)/2)*filterwidth/400;

fy=fft(yvector);
lft1=1:(length(fy)/2);
lft2=(length(fy)/2+1):length(fy);
% Compute filter shape.
if strcmp(filtermode,'Band-pass')
       ffilter1=shape(lft1,centerfrequency+1,filterWidth,filtershape);
       ffilter2=shape(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
       ffilter=[ffilter1,ffilter2];
end
if strcmp(filtermode,'High-pass')
       centerfrequency=length(xvector)/2;
       ffilter1=shape(lft1,centerfrequency+1,filterWidth,filtershape);
       ffilter2=shape(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
       ffilter=[ffilter1,ffilter2];
end
if strcmp(filtermode,'Low-pass')
       centerfrequency=0;
       ffilter1=shape(lft1,centerfrequency+1,filterWidth,filtershape);
       ffilter2=shape(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
       ffilter=[ffilter1,ffilter2];
end
if strcmp(filtermode,'Band-reject (notch)')
       ffilter1=shape(lft1,centerfrequency+1,filterWidth,filtershape);
       ffilter2=shape(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
       ffilter=1-[ffilter1,ffilter2]; 
end
if strcmp(filtermode,'Comb pass')
    n=2;
    ffilter1=shape(lft1,centerfrequency+1,filterWidth,filtershape);
    ffilter2=shape(lft2,length(fy)-centerfrequency+1,filterwidth,filtershape);
    while n<50
       ffilter1=ffilter1+shape(lft1,n*(centerfrequency+1),filterWidth,filtershape);
       ffilter2=ffilter2+shape(lft2,length(fy)-n*(centerfrequency+1),filterWidth,filtershape);
       n=n+1;
    end
       ffilter=[ffilter1,ffilter2];
end
if strcmp(filtermode,'Comb notch')
    n=2;
    ffilter1=shape(lft1,centerfrequency+1,filterWidth,filtershape);
    ffilter2=shape(lft2,length(fy)-centerfrequency+1,filterWidth,filtershape);
    while n<30
       ffilter1=ffilter1+shape(lft1,n*(centerfrequency+1),filterWidth,filtershape);
       ffilter2=ffilter2+shape(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));

subplot(3,1,1)  % Plot original signal in top plot
plot(xvector,yvector,'b');
title('Fourier Filter Tool')
axis([xvector(1) xvector(length(xvector)) min(yvector) max(yvector)]);  

subplot(3,1,3)  % Plot filtered signal in lower plot
plot(xvector,ry,'r');
grid
if ymode
   axis([xvector(1) xvector(length(xvector)) min(ry) max(ry)]);  
else
  axis([xvector(1) xvector(length(xvector)) min(yvector) max(yvector)]);    
end
% title([ filtermode ' mode:  Freq= ' num2str(centerfrequency./range(xvector)) '    Period= ' num2str(range(xvector)./centerfrequency) '   Width= ' num2str(filterwidth./range(xvector)) '     Shape=  ' num2str(filtershape)])
title('Output of '+filtermode+' filter')

% xlabel(['OUTPUT: Filtered signal    x=time   Width = ' num2str(halfwidth(xvector,ry))  '    Area = ' num2str(trapz(xvector,ry)) '    STD = ' num2str(std(ry))]);
xlabel(['Filter Width = ' num2str(filterwidth)  '    Area = ' num2str(trapz(xvector,ry)) '    STD = ' num2str(std(ry))]);


if ymode
    ylabel('Fixed Y scale.')
else
    ylabel('Variable Y scale.')
end
subplot(3,1,2)    % Plot power spectrum (blue) and filter (red) in middle plot
py=fy .* conj(fy); % Compute power spectrum
plotrange=2:length(fy)/2;
if xmode
   f=range(xvector)./(plotrange-1);
else
   f=((plotrange-1)./range(xvector));
end
rpypr=real(py(plotrange));
switch plotmode
  case 1
    plot(f,rpypr,f,max(rpypr).*ffilter(plotrange),'r')
    ylabel('Linear y')
  case 2
    semilogx(f,rpypr,f,max(rpypr).*ffilter(plotrange),'r')
    ylabel('Linear y')
  case 3
    semilogy(f,rpypr,f,max(rpypr).*ffilter(plotrange),'r')
    ylabel('Log y')
  case 4
    loglog(f,rpypr,f,max(rpypr).*ffilter(plotrange),'r')
    ylabel('Log y')
    otherwise
end
grid
% title('Output of '+filtermode+' filter')
title('Fourier spectra:  BLUE = Input signal    RED = '+filtermode+' filter')
if xmode
    xlabel('x=Time (e.g. seconds)')
else
    xlabel('x=Frequency (e.g. cycles/second)')
end
% Expand frequency axis if filter is in the lower half of frequency axis
if centerfrequency+filterwidth/2<length(fy)/4
    axis([0 max(f)/2 min(py) 1.1*max(rpypr)]) 
else
    axis([0 max(f) min(py) 1.1*max(rpypr)])
end

if PlayOutputAsSound==1
    sound(ry)
end
PlayOutputAsSound=0;

function g = shape(x,pos,wid,n)
%  shape(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: shape([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

function r = range(arr)
r = max(arr) - min(arr);
end