%% Fourier self-deconvolution
% A tool for peak sharpening by Fourier saelf-deconvolution. Note:To see the graphs on the right, right-click on the right panel and select "Disable synchronous scrolling". Tom O'Haver 4/5/2023

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

%


FrequencySpectra=false; % Check this box to display frequency spectra 
PlotAllSteps=false; % Check this box to plot intermediate steps  
NumPoints=length(mydata);

% Set the x-axis scale expansion beginning and end points (0 - 100%)
startpc=3.1; % Percentage of data points to start data selection
endpc=95.6;% Percentage of data points to end data selection
startpoint=round(startpc*NumPoints./100);
endpoint=round(endpc*NumPoints./100);
x=mydata(startpoint:endpoint,1);
y=mydata(startpoint:endpoint,2);
 
% Deconvolution function shape
PeakShape="Gaussian";
PCGaussian=14;

% dw = Deconvolution width (half-width of deconvolution function df)
dw=1.7 ; % Double-click sliders to change their ranges if necessary
% Percent denominator addition
DA=3;
% Percentage of entire spectrum included for Fourier filter
FrequencyCutoff=14;
% for Fourier filter
CutOffRate=3;

if PeakShape=="Gaussian"
    df=gaussian(x,min(x),dw)+gaussian(x,max(x),dw);
    ydc=ifft(fft(y)./fft(df)).*sum(df);  % Deconvolution by fft/ifft
end

if PeakShape=="Lorentzian"
    df=lorentzian(x,min(x),dw)+lorentzian(x,max(x),dw);
    ydc=ifft(fft(y)./fft(df)).*sum(df);  % Deconvolution by fft/ifft
end

if PeakShape=="Logistic"
    df=logistic(x,min(x),dw)+logistic(x,max(x),dw);
    ydc=ifft(fft(y)./fft(df)).*sum(df);  % Deconvolution by fft/ifft
end

if PeakShape=="GL blend"
    df=GL(x,min(x),dw,PCGaussian)+GL(x,max(x),dw,PCGaussian);
    ydc=ifft(fft(y)./fft(df)).*sum(df);  % Deconvolution by fft/ifft
end

% Apply only denominator addition to reduce noise and ringing
fftc=fft(df);
% Deconvolution by fft/ifft with Denominator addition
ydcDA=ifft(fft(y)./(fftc+DA.*0.01.*max(fftc))).*sum(df);

% Apply only Fourier filtering to reduce noise and ringing
% Low-pass mode
sy=FourierFilter(x,ydc,0,FrequencyCutoff.*0.01.*length(x),CutOffRate,1);

% Apply both denominator addition and Fourier filtering to reduce noise
%  and ringing; Low-pass mode
syDA=FourierFilter(x,ydcDA,0,FrequencyCutoff.*0.01.*length(x),CutOffRate,1);

% Plot all the steps
if PlotAllSteps==true
    figure(1)
    clf
    subplot(2,2,1);
    plot(x,y,'k'); title(['Original data']);grid
    xlabel(['Deconvolution width = ' num2str(dw) '%']  );

    subplot(2,2,3);
    plot(x,ydcDA,'b');grid
    title('FSD, with denominator addition only');
    xlabel(['Denom. Addition = ' num2str(DA) '%']  );

    subplot(2,2,2);
    plot(x,sy,'g');grid
    title('FSD with filtering only');
    xlabel(['FrequencyCutoff = ' num2str(FrequencyCutoff) ]);

    subplot(2,2,4);
    plot(x,syDA,'r');grid
    title('Both filtering and denom. addition');
else
    figure(1)
    clf
    plot(x,y,'k',x,syDA); title('Original data');grid
    legend('Original data','Self-Deconvoluted')
end
AreaRecovery=100.*trapz(x,syDA)/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(ydcDA);
    syc=fyc .* conj(fyc);
    psrange=1:length(fyc)/2;
    psyc=real(syc(psrange));
    loglog(psrange,psyc,'b','linewidth',1)

    hold on
    fyc=fft(sy);
    syc=fyc .* conj(fyc);
    psyc=real(syc(psrange));
    loglog(psrange,psyc,'g','linewidth',1)

    fyc=fft(y);
    syc=fyc .* conj(fyc);
    minspec=min(syc);maxspec=max(syc);
    psyc=real(syc(psrange));
    loglog(psrange,psyc,'k','linewidth',1)

    fyc=fft(syDA);
    syc=fyc .* conj(fyc);
    psyc=real(syc(psrange));
    loglog(psrange,psyc,'r','linewidth',1)

    hold off
    grid
    axis([1 NumPoints/4 minspec maxspec])
    title('Frequency spectra')
    xlabel('Frequency')
    legend('Denominator addition only','Filtering only','Original signal','Both filtering and denom. addition','Location','SouthWest')
else
end

%

function g = gaussian(x,pos,wid)
%  gaussian(x,pos,wid) = gaussian peak centered on pos, half-width=wid
%  x may be scalar, vector, or matrix, pos and wid both scalar
%  T. C. O'Haver, 1988
% Examples: gaussian([0 1 2],1,2) gives result [0.5000    1.0000    0.5000]
% plot(gaussian([1:100],50,20)) displays gaussian band centered at 50 with width 20.
g = exp(-((x-pos)./(0.60056120439323.*wid)) .^2);
end

function g = lorentzian(x,position,width)
% lorentzian(x,position,width) Lorentzian function.
% where x may be scalar, vector, or matrix
% position and width scalar
% T. C. O'Haver, 1988
% Example: lorentzian([1 2 3],2,2) gives result [0.5 1 0.5]
g=ones(size(x))./(1+((x-position)./(0.5.*width)).^2);
end

function g = logistic(x,pos,wid)
% logistic function.  pos=position; wid=half-width (both scalar)
% logistic(x,pos,wid), where x may be scalar, vector, or matrix
% pos=position; wid=half-width (both scalar)
% T. C. O'Haver, 1991
n = exp(-((x-pos)/(.477.*wid)) .^2);
g = (2.*n)./(1+n);
end

function g = GL(x,pos,wid,m)
% function g = GL(x,pos,wid,m)
% Unit-height Gaussian/Lorentzian blend.
% pos=position; wid=half-width;
% m = percent Gaussian character.
%  T. C. O'Haver, 2012
g=(m/100)*gaussian(x,pos,wid)+(1-(m/100))*lorentzian(x,pos,wid);
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
