function [SPI, SPI_smooth, H1] = Spectral_Purity_Sliding_Window(PPG,fs,Aff,L)
% Inputs:
% PPG: Input signal (not necessary PPG)
% fs: sampling frequency
% Aff: 1 for display
% Outputs:
% - SPI
% - SPI smoothed (alpha parma in the function)
% - H1: first Hjorth descriptor
% Compute spectral purity index (p. 102, Sornmo Bioelectrical Signal
% Processing)
% Time-domain estimates of spectral moments


%% Miror
PPG = PPG(:);
PPG = [flipud(PPG(1:L+1,:)); PPG];


%% Using sliding window of size L
K = size(PPG,2);

SPI = zeros(size(PPG)); SPI_smooth = zeros(size(PPG));
H1 = zeros(size(PPG)); H2 = zeros(size(PPG));
d1 = zeros(size(PPG)); d2 = zeros(size(PPG));
w0 = zeros(size(PPG)); w2 = zeros(size(PPG)); w4 = zeros(size(PPG));

alpha = 0.95;

for n=L+1:length(PPG(:,1))-1    
    for j=1:K % For each column of the input data matrix
    
    d1(n,j) = PPG(n,j)-PPG(n-1,j);
    d2(n,j) = PPG(n+1,j) -2*PPG(n,j) + PPG(n-1,j);
    
    % zero-, second-, and fourth-order moments
    w0(n,j) = (2*pi/L)* sum(PPG(n-L+1:n,j).^2);
    w2(n,j) = (2*pi/L)* sum(d1(n-L+1:n,j).^2);
    w4(n,j) = (2*pi/L)* sum(d2(n-L+1:n,j).^2);
    
    SPI(n,j) = w2(n,j)^2/(w0(n,j)*w4(n,j));
    
    SPI_smooth(n,j) = alpha*(SPI_smooth(n-1,j))+ (1-alpha)*SPI(n,j);
    
      H1(n,j) = sqrt(w2(n,j)/w0(n,j));
%     H2(n,j) = sqrt((w4(n)/w2(n))-(w2(n)/w0(n)));    
    end
end


%% 
PPG = PPG((L+2):end,:);
SPI = SPI((L+2):end,:);
SPI_smooth = SPI_smooth((L+2):end,:);
H1 = H1((L+2):end,:);

% ----------------------------
if(Aff ==1)

axe1(1) = subplot(211);
plot((1:length(PPG(:,1))),PPG); title('Input signal'); grid on;
axe1(2) = subplot(212);
plot((1:length(SPI_smooth(:,1))),SPI_smooth); grid on; title('SPI');
set(gca,'YLim',[0 1.1]);
linkaxes(axe1,'x');

end

