function [ F ] = frequency_features( ecg, fs )

% 
% Copyright (C) 2017 
% Shreyasi Datta
% Chetanya Puri
% Ayan Mukherjee
% Rohan Banerjee
% Anirban Dutta Choudhury
% Arijit Ukil
% Soma Bandyopadhyay
% Rituraj Singh
% Arpan Pal
% Sundeep Khandelwal
% 
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.


%FREQUENCY_FEATURES Summary of this function goes here
%   Detailed explanation goes here

% 2 second window 
    winlen = fs*2;  
    winnum = 2*floor(length(ecg)/winlen) - 1;
    window_spec = [];
    if (winnum<=2)
        F = zeros(1, 12);
    else
    
    for j = 1:winnum
        temp = ecg((j-1)*winlen/2+1:(j-1)*winlen/2+winlen);
        z = abs(fft(temp));
        z = z(1:end/2);
        pow = z.*z;
        bin = (1:length(z)).';
        if (j~=1)
        % frequency range 0-100, 100-200, 200-300, 300-400; 400-500 Hz
        freq_range = [100;150];
        %freq_range = [50;100;150;200;250];
        bin_range = freq_range*length(z)/fs;
        % calculate total spectral power
        PT = sum(pow);
       % PT = 1;
        % calculate normalized power in each of the 4 spetrams
        P1 = sum(pow(1:bin_range(1)))/PT;
        P2 = sum(pow(bin_range(1):bin_range(2)))/PT;
        b = 1:length(z);
        centroid = sum(b*z)/sum(z);
        rolloff = 0.85*sum(z);
        flux = norm(z - z_prev);
        kur = kurtosis(temp);
        sk = skewness(temp);
        % define the window
        wintype = 'hamming';
        winlen = length(temp);
        winamp = [0.5,1]*(1/winlen);
        window_spec = [window_spec; centroid rolloff flux kur];
       
        end
        z_prev = z;
    end
    
    % mean window spec feature
    
    spec1 = mean(window_spec(:,1:end-1));
    kur1 = window_spec(:,end);
    ttt = find(isnan(kur1)==1);
    kur1(ttt) = [];
   kur_mean = mean(kur1);
   % window_spec_mean = mean(window_spec);
   window_spec_mean = [spec1 kur_mean]; 

level = 5;
[c,l] = wavedec(ecg,level,'db4');

%a5 = appcoef(c,l,'db4');
%det5 = detcoef(c,l,level);
%det4 = detcoef(c,l,level-1);
det3 = detcoef(c,l,level-2);
%--------------------------------------------------------------------------
f7 = log2(var(det3));                                       %7
%--------------------------------------------------------------------------
% xx = downsample(ecg, 8);
% [f14, ~, f15] = Spectral(xx);                               %13 14
%--------------------------------------------------------------------------
[lp,g] = lpc(double(ecg),10);
f161821232425 = lp([2 4 7 9 10 11]);                        %16 18 21 23 24 25
%--------------------------------------------------------------------------
N = length(ecg);
xdft = fft(ecg);
% if rem(N,2) ~= 0
%     N = N - 1;
% end
xdft = xdft(1:N/2+1);
psdx = (1/(2*pi*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);

f = fs*(0:(N/2))/N;
frequency_centroid = (sum(f'.*(psdx.^2)))/(sum(psdx.^2));   %26

F = [window_spec_mean f7 f161821232425 frequency_centroid];
    end
