function [Hratio,wb]=Hratio(x,fmax,numfreq)
% function [HIndex]=IndexOfHarmonicity(x,numfreq);
% 
% INPUT, signal, fmax is max frequency that can be considered as main
% frequency = around 1 for stride and around 2 for steps. 
% numfreq = number of superharmonics in analysis
% OUTPUT: Harmoncity Index and mainfrequency
% PLOT: with no output a plot is generated
% CALLS: WndMainFreq.m
% CL 2016


if nargin<3 || isempty(numfreq), numfreq = 10; end

fs = 100;% sample frequency
le=length(x);
temp=[];
   n=10*length(x);
   ws=fix(length(x)/2);
   no=fix(0.99*ws);

[Px,fx]=spectrum2(x,n,no,boxcar(le),100,'mean');
Px = Px/sum(Px); % All power spectral densities were normalized by dividing the power by the sum of the total power spectrum, which equals the variance. 

% select  frequency
wb = WindMainFreq(x,100,fmax) ; % mainfrequency frequency
wbr= wb + 0.055;  % % basis mainfrequency + tail
wbl = wb - 0.055; 

wind=find(fx >=wbl & fx<= wbr); % index van basisfrequentie

PI_0=sum(Px(wind));

temp=PI_0;
for k=1:numfreq;
    w=wb*(k+1);
    wr=w+0.055;
    wl=w-0.055;
    wInd=find(fx >=wl & fx<=wr);
    if nargout ==0;      % controle 
        plot(fx,Px,'r'); hold on;
        set(gca','XLim', [0 6],'FontSize', 8)
        xlabel('Frequency (Hz)');  ylabel('Power');
        plot(wb,max(Px),'r.','MarkerSize', 10),
        plot(wbr,0,'c.'); plot(wbl,0,'c.');
        plot(w,0,'r.','MarkerSize', 10)
        plot(wl,0,'c.'); plot(wr,0,'c.');
    end

    PI=sum(Px(wInd));
    temp=[temp; PI];
    
end

 Hratio = sum(temp(1:2:end-1))/sum(temp(2:2:end));