%Estimates the bandwidth (bw) or sigma of a Gaussian kernel for kernel
%density estimation (KDE) using a total variation (TV) constraint.
%Input arguments: histoInput - a normalised histogram of the data.
% noDataPoints - number of datapoints
% lambda - the parameter in Lagrangian minimization (free parameter)
% sampleRate - sample rate of the original data ("time" domain).
function [bw]=TVEstBwExhaustivePub(histoInput,noDataPoints,lambda,sampleRate)
    nF=1:round(noDataPoints/2);
    OmegaS=2*pi/sampleRate;
    fourierIndCTFT=linspace(-OmegaS/(2),OmegaS/(2),noDataPoints);
    fourierIndCTFTOrigin=linspace(0,OmegaS/2,length(nF));
    histoFFT=fft(histoInput,noDataPoints);
    tv1=zeros(1,length(nF))*NaN;
    residual=zeros(1,length(nF))*NaN;
    for ff=nF(2:end)
       %do KDE with the kernel of the bandwidth and the data
       curEst=myKDEFourier(fourierIndCTFTOrigin(ff)^-1,histoFFT,sampleRate,fourierIndCTFT);
       %normalise
       curEst=curEst/sum(curEst);
       %compute TV for the estimate
       tv1(ff)=myTV(curEst);
       residual(ff)=sum(abs(histoInput-curEst));
    end
    TVconstraint=residual+lambda*tv1;
    [~,ind] = min(TVconstraint);
    bw=fourierIndCTFTOrigin(ind)^-1;
end

function out=myTV(in)
    out=sum(abs(in(2:end)-in(1:end-1)));
end

function out=myKDEFourier(sigma,data,sampleRate,fourierIndCTFT)
    kernFreqThreshFFT=sampleRate^-1*(sigma*sqrt(2*pi))*ifftshift(pdf('norm',fourierIndCTFT,0,sigma^-1));
    out=abs(ifft(data.*kernFreqThreshFFT,length(data)));
end