function [ new_ecg ,fsig,fecg] = QRS_Energy( ecg, fs,highpass ,Aff)
%QRS_ENERGY Summary of this function goes here
%preprocessing
 if highpass,
    ecg=ecg-mean(ecg);           
    %%Filtering Section design and implementation
    Cuttoff=4;
    Wp=(Cuttoff/fs)*2;
    b=fir1(100,Wp,'high');
    ecg = filtfilt(b,1,ecg);        
    Cuttoff=40; Wp=(Cuttoff/fs)*2;
    b=fir1(Cuttoff,Wp,'low');
    ecg = filtfilt(b,1,ecg);  
%     Transition=Cuttoff+8;
%     Wp=(Cuttoff/fs)*2;    % Cuttoff frequency
%     Ws=(Transition/fs)*2; % Transtions phase
%     Rp=0.5; %Ripple
%     RS=20; % attenuation in the stopband
%     [N,Wn] = buttord(Wp,Ws,Rp,RS);
%     [b,a] = butter(N,Wn) ;
%     ecg = filtfilt(b,a,ecg);
%     ecg = (ecg - min(ecg)) / ( max(ecg) - min(ecg) );
 end

%   Detailed explanation goes here
ecg=ecg-mean(ecg);
fecg=ecg;
qWin=round(.1*fs);
lWin=round(1*fs);
global win;
win=floor(lWin/2);
dur=floor(qWin/2);

qE=zeros(length(ecg),1);
lE=zeros(length(ecg),1);
fsig=zeros(length(ecg),1);
% fmax=ones(length(ecg),1);
power=2;
%calculate the Energy for the first possible window
 i=win+1;
qE(i)=sum(abs(ecg(i-dur:i+dur)).^power);
lE(i)=sum(abs(ecg(i-win:i+win)).^power);
fsig(i)=qE(i)/lE(i);
%optimized calculation of the energy for the rest of the singal
for i=win+2:length(ecg)-win-1
    qE(i)=qE(i-1)+(ecg(i+dur)^power)-(ecg(i-1-dur)^power);
    lE(i)=lE(i-1)+(ecg(i+win)^power)-(ecg(i-1-win)^power);      
%     fsig(i)=qE(i)/lE(i);  
%      fsig(i)=log(qE(i)/lE(i));  
    fsig(i)=min(qE(i)/lE(i),1);
%     fsig(i)=max(qE(i)-lE(i),0);
%     fsig(i)=min(qE(i)/0.9*lE(i),1);
end
% for i=5*win+2:length(ecg)-5*win-1
%     fmax(i)=max(ecg(i-5*win:i+5*win));
% end

if highpass
% % %     mfsig=mean(fsig); ifsig=iqr(fsig);
% % %     mfsig=mean(fsig(fsig>mfsig+ifsig)); ifsig=iqr(fsig(fsig>mfsig+ifsig)); 
% % %     fisgmax=max(fsig(fsig<=mfsig+ifsig));
% % %     new_ecg1=(fsig/fisgmax).*ecg;
% % %     new_ecg=(fsig/max(fsig)).*ecg;
% % %     new_ecg= (new_ecg - min(new_ecg1)) / ( max(new_ecg1) - min(new_ecg1) );
% % %     new_ecg=new_ecg-mean(new_ecg);    
% % %     figure(21);
% % %     plot(new_ecg); hold on; plot((fsig/max(fsig)).*ecg,'r'); hold off;



    new_ecg=(fsig/max(fsig)).*ecg;
%         new_ecg=(fsig./fmax).*ecg; new_ecg(isnan(new_ecg))=0; new_ecg(new_ecg==inf)=0; new_ecg(new_ecg==-inf)=0;
%     new_ecg= (new_ecg - min(new_ecg)) / ( max(new_ecg) - min(new_ecg) );
    new_ecg=new_ecg-mean(new_ecg);
    
%     mfsig=mean(fsig); ifsig=iqr(fsig);
%     mfsig=mean(fsig(fsig>mfsig+ifsig)); ifsig=iqr(fsig(fsig>mfsig+ifsig)); 
%     fisgmax=max(fsig(fsig<=mfsig+ifsig));
%     new_ecg= (new_ecg - min(new_ecg)) / ( fisgmax - min(new_ecg) );
%     new_ecg=new_ecg-mean(new_ecg);
else
    new_ecg=(fsig).*ecg;
end

if(Aff)
    figure(10);
    plot(ecg);
    hold on;
    plot(qE/max(qE),'r');
    plot(lE/max(lE),'k');
    hold off;


    figure(11);
    plot(ecg);
    hold on;
    plot(fsig/max(fsig),'r');
    plot((fsig/max(fsig)).*ecg,'k');
    hold off;

    figure(12);
    plot(ecg); hold on; plot(new_ecg,'r'); hold off;
end


