function features=calculateSSFeatures(ecg)
Fs=300;
[pxx,f]=pwelch(ecg,256,[],[],Fs);
features(1)=real(sum(pxx(f>=5 & f<=15))/sum(pxx(f>=5 & f<=40)));
features(2)=real(sum(pxx(f>=1 & f<=40))/sum(pxx(f<=40)));
% features(3)=mean(ecg);
% features(4)=var(ecg);
features(3)=skewness(ecg);
features(4)=kurtosis(ecg);

fs2=35;
Cuttoff=40; Wp=(Cuttoff/Fs)*2; % Cuttoff frequency to filter ECG
b=fir1(100,Wp); a=1;
ECG = filtfilt(b,a,ecg);
ECG = zscore(ECG);
ECG = resample(ECG,fs2,Fs);
[trans1, ~] = Moving_Average(ECG,5);
[~, SPI_smooth, ~] = PPG_Spectral_Purity_recursive(trans1(:),fs2,0,2*fs2);

L=3; SPI_stat= zeros(15-L+1,2);
ll=floor(length(SPI_smooth)/fs2)-3;

for j=1:ll % To compute mean and max SPI on L sec windows
    % j
    SPI_stat(j,1) = mean(SPI_smooth(1+((j-1)*fs2):(1+(j-1)*fs2)+(L*fs2),1));
    %     SPI_stat(j,2) = mean(SPI_smooth((5*fs2)+((j-1)*fs2):(5*fs2)+((j-1)*fs2)+(L*fs2),2));
    SPI_stat(j,2) = max(SPI_smooth(1+((j-1)*fs2):1+(j-1)*fs2+(L*fs2),1));
    %     SPI_stat(j,4) = max(SPI_smooth((5*fs2)+((j-1)*fs2):(5*fs2)+((j-1)*fs2)+(L*fs2),2));
end
DMAX_II = 0; DMAX_V = 0;
[MAX_II, indII] = max(SPI_stat(:,1)); [MAX_V, indV] = max(SPI_stat(:,2));
if(indII > 2)
    DMAX_II = MAX_II - SPI_stat(indII-2,1); % Reflects an increase of SPI
end
if(indV > 2)
    DMAX_V = MAX_V - SPI_stat(indV-2,2); % Reflects an increase of SPI
end

FEATURES(1) = max(SPI_stat(:,1));
FEATURES(2) = max(SPI_stat(:,2));
FEATURES(3) = DMAX_II;
FEATURES(4) = DMAX_V;



SPI_mean = zeros(8,1);
ll2=ll;
for j=1:ll2 % 3 sec window, 1 sec overlap
    SPI_mean(j,1) = mean(SPI_smooth(1+((j-1)*fs2):1+((j-1)*fs2)+(3*fs2),1));
end

features=[FEATURES(:)' max(SPI_mean) features(:)'];