--- a
+++ b/preprocessOfApneaECG/BioSigKit/Algorithms/nqrsdetect.m
@@ -0,0 +1,430 @@
+function QRS=nqrsdetect(S,fs)
+% nqrsdetect - detection of QRS-complexes
+%
+%   QRS=nqrsdetect(S,fs);
+%
+% INPUT
+%   S       ecg signal data
+%   fs      sample rate
+%
+% OUTPUT
+%   QRS     fiducial points of qrs complexes
+%
+%
+% see also: QRSDETECT
+%
+% Copyright (C) 2006 by Rupert Ortner
+%
+%% 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+%% USA
+%% ============== Now Part of BioSigToolkit ============= %%
+%% Modified and imporved by Hooman Sedghamiz, Feb 2018.
+
+S=S(:);
+S=full(S);
+N=round(fs);   %Filter order
+%---------------------------------------
+%Replaces filter bank in [1]
+Bw=5.6;     %filter bandwidth
+Bwn=1/(fs/2)*Bw;    
+M=round((fs/2)/Bw); %downsampling rate
+
+Wn0=Bwn;    %bandwidth of the first filter
+Wn1=[Bwn 2*Bwn];    %bandwidth of the second filter
+Wn2=[2*Bwn 3*Bwn];
+Wn3=[3*Bwn 4*Bwn];
+Wn4=[4*Bwn 5*Bwn];
+
+h0=fir1(N,Wn0); %impulse response of the first filter
+h1=fir1(N,Wn1,'bandpass');
+h2=fir1(N,Wn2,'bandpass');
+h3=fir1(N,Wn3,'bandpass');
+h4=fir1(N,Wn4,'bandpass');
+
+%Polyphase implementation of the filters
+y=cell(1,5);    
+y{1}=polyphase_imp(S,h0,M); %W0 (see [1]) filtered and downsampled signal
+y{2}=polyphase_imp(S,h1,M); %W1 
+y{3}=polyphase_imp(S,h2,M); %W2
+y{4}=polyphase_imp(S,h3,M); %W3
+y{5}=polyphase_imp(S,h4,M); %W4
+%----------------------------------------------
+
+cut=ceil(N/M);  %Cutting off of initial transient because of the filtering
+y1=[zeros(cut,1);y{1}(cut:length(y{1}))];  
+y2=[zeros(cut,1);y{2}(cut:length(y{2}))];
+y3=[zeros(cut,1);y{3}(cut:length(y{3}))];
+y4=[zeros(cut,1);y{4}(cut:length(y{4}))];
+y5=[zeros(cut,1);y{5}(cut:length(y{5}))];
+%----------------------------------------
+
+P1=sum([abs(y2) abs(y3) abs(y4)],2); %see [1] equation (13)
+P2=sum([abs(y2) abs(y3) abs(y4) abs(y5)],2);
+P4=sum([abs(y3) abs(y4) abs(y5)],2);
+
+FL1=MWI(P1); %Feature 1 according to Level 1 in [1]
+FL2=MWI(P2); %Feature 2 according to Level 2
+FL4=MWI(P4); %Feature 4 according to Level 4
+%--------------------------------------
+%Level 1 [1]
+d=sign(diff(FL1));
+d1=[0;d];
+d2=[d;0];
+f1=find(d1==1);
+f2=find(d2==-1);
+EventsL1=intersect(f1,f2); %Detected events
+%-------------------------------------------------------
+%Level 2 [1]
+meanL1=sum(FL2(EventsL1),1)/length(EventsL1);
+NL=meanL1-meanL1*0.1;   %Start Noise Level
+SL=meanL1+meanL1*0.1;   %Start Signal Level
+threshold1=0.08;    %Threshold detection block 1
+threshold2=0.7;     %Threshold detection block 2
+[~,~,DS1,Class1]=detectionblock(FL2,EventsL1,NL,SL,threshold1);
+[~,~,DS2,Class2]=detectionblock(FL2,EventsL1,NL,SL,threshold2);
+%---------------------------------------------------
+ClassL3=[]; 
+for i=1:length(EventsL1)
+    C1=Class1(i);
+    C2=Class2(i);
+    if C1==1
+        if C2==1
+            ClassL3=[ClassL3 1];   %Classification as Signal
+        else
+            delta1=(DS1(i)-threshold1)/(1-threshold1);
+            delta2=(threshold2-DS2(i))/threshold2;
+            if delta1>delta2
+                ClassL3=[ClassL3 1]; %Classification as Signal
+            else
+                ClassL3=[ClassL3 0];  %Classification as Noise
+            end
+        end
+    else
+        if C2==1
+            ClassL3=[ClassL3 1]; %Classification as Signal
+        else
+            ClassL3=[ClassL3 0];  %Classification as Noise
+        end
+    end
+end
+SignalL3=EventsL1(find(ClassL3));   %Signal Level 3
+NoiseL3=EventsL1(find(ClassL3==0)); %Noise Level 3
+%--------------------------------------------
+%Level 4 [1]
+threshold=0.3;
+VSL=(sum(FL4(SignalL3),1))/length(SignalL3); 
+VNL=(sum(FL4(NoiseL3),1))/length(NoiseL3);   
+SL=(sum(FL4(SignalL3),1))/length(SignalL3);   %Initial Signal Level
+NL=(sum(FL4(NoiseL3),1))/length(NoiseL3);      %Initial Noise Level
+SignalL4=[];
+NoiseL4=[];
+DsL4=[];   %Detection strength Level 4
+for i=1:length(EventsL1)
+    Pkt=EventsL1(i);    
+    if ClassL3(i)==1;   %Classification after Level 3 as Signal
+       SignalL4=[SignalL4,EventsL1(i)]; 
+       SL=history(SL,FL4(Pkt));       
+       Ds=(FL4(Pkt)-NL)/(SL-NL);       %Detection strength
+       if Ds<0
+           Ds=0;
+       elseif Ds>1
+           Ds=1;
+       end
+       DsL4=[DsL4 Ds];
+    else        %Classification after Level 3 as Noise
+       Ds=(FL4(Pkt)-NL)/(SL-NL);   
+       if Ds<0
+           Ds=0;
+       elseif Ds>1
+           Ds=1;
+       end
+       DsL4=[DsL4 Ds];
+       if Ds>threshold          %new classification as Signal 
+           SignalL4=[SignalL4,EventsL1(i)];  
+           SL=history(SL,FL4(Pkt));     
+       else                      %new classification as Noise
+           NoiseL4=[NoiseL4,EventsL1(i)];
+           NL=history(NL,FL4(Pkt));      
+       end
+   end
+end
+%------------------------------------------------
+%Level 5  
+%if the time between two RR complexes is too long => go back and check the
+%events again with lower threshold
+SignalL5=SignalL4;
+NoiseL5=NoiseL4;
+periods=diff(SignalL4);
+M1=100;
+a=1;
+b=1/(M1)*ones(M1,1);
+meanperiod=filter(b,a,periods); %mean of the RR intervals
+SL=sum(FL4(SignalL4))/length(SignalL4);
+NL=sum(FL4(NoiseL4))/length(NoiseL4);
+threshold=0.2;
+for i=1:length(periods)
+    if periods(i)>meanperiod*1.5     %if RR-interval is to long
+        intervall=SignalL4(i):SignalL4(i+1);
+        critical=intersect(intervall,NoiseL4);   
+        for j=1:length(critical)
+            Ds=(FL4(critical(j))-NL)/(SL-NL); 
+            if Ds>threshold         %Classification as Signal
+                SignalL5=union(SignalL5,critical(j));   
+                NoiseL5=setxor(NoiseL5,critical(j));
+            end
+        end
+    end
+end
+%---------------------------------------------------
+%Umrechnung auf Originalsignal (nicht downgesamplet)
+Signaln=conversion(S,FL2,SignalL5,M,N,fs);
+%----------------------------------------------------
+%Level 6 If interval of two RR-complexes <0.24 => go back and delete one of them
+height=FL2(SignalL5);   
+Signal=Signaln;
+temp=round(0.1*fs);
+difference=diff(Signaln);  %Difference between two signal points
+k=find(difference<temp);
+for i=1:length(k)
+    pkt1=SignalL5(k(i));
+    pkt2=SignalL5(k(i)+1);
+    verg=[height(k(i)),height(k(i)+1)]; 
+    [x,j]=max(verg);    
+    if j==1
+        Signal=setxor(Signal,Signaln(k(i)+1)); %Deleting first Event
+    else
+        Signal=setxor(Signal,Signaln(k(i))); %Deleting second Event
+    end
+end
+QRS=Signal;
+
+
+%-------------------------------------------------------------------
+%-------------------------------------------------------------------
+%-------------------------------------------------------------------
+%subfunctions
+
+function y=MWI(S)
+
+% MWI - Moving window integrator, computes the mean of two samples
+%   y=MWI(S)
+%
+% INPUT
+%   S       Signal
+%
+% OUTPUT
+%   y       output signal
+a=[0;S];
+b=[S;0];
+c=[a,b];
+y=sum(c,2)/2;
+y=y(1:length(y)-1);
+%------------------------------------------------
+function y=polyphase_imp(S,h,M)
+
+% polyphase_imp - polyphase implementation of decimation filters [2]
+%   y=polyphase_imp(S,h,M)
+%
+% INPUT
+%   S       ecg signal data
+%   h       filter coefficients
+%   M       downsampling rate
+%
+% OUTPUT
+%   y       filtered signal
+%
+
+%Determining polyphase components ek
+e=cell(M,1);
+l=1;
+m=mod(length(h),M);
+while m>0
+    for n=1:ceil(length(h)/M)
+        el(n)=h(M*(n-1)+l);
+    end
+    e{l}=el;
+    l=l+1;
+    m=m-1;
+end
+clear el;
+for i=l:M
+    for n=1:floor(length(h)/M)
+        el(n)=h(M*(n-1)+i);      
+    end
+    e{i}=el;
+end
+%Filtering
+max=ceil((length(S)+M)/M); 
+Sdelay=S;
+for i=1:M
+    Sd=downsample(Sdelay,M);
+    a=filter(e{i},1,Sd);     
+    if length(a)<max
+        a=[a;zeros(max-length(a),1)]; 
+    end
+    w(:,i)=a;
+    Sdelay=[zeros(i,1);S];
+end
+y=sum(w,2);
+%----------------------------------------------------------
+function [Signal,Noise,VDs,Class]=detectionblock(mwi,Events,NL,SL,threshold)
+
+% detectionblock - computation of one detection block 
+%
+%   [Signal,Noise,VDs,Class]=detectionblock(mwi,Events,NL,SL,threshold)
+%
+% INPUT
+%   mwi         Output of the MWI
+%   Events      Events of Level 1 (see [1])
+%   NL          Initial Noise Level
+%   SL          Initial Signal Level
+%   threshold   Detection threshold (between [0,1])
+%
+% OUTPUT
+%   Signal      Events which are computed as Signal
+%   Noise       Events which are computed as Noise
+%   VDs         Detection strength of the Events
+%   Class       Classification: 0=noise, 1=signal
+
+Signal=[];
+Noise=[];
+VDs=[];
+Class=[];
+sumsignal=SL;
+sumnoise=NL;
+for i=1:length(Events)
+    P=Events(i);
+    Ds=(mwi(P)-NL)/(SL-NL); %Detection strength
+    if Ds<0
+        Ds=0;
+    elseif Ds>1
+        Ds=1;
+    end
+    VDs=[VDs Ds];
+    if Ds>threshold     %Classification as Signal
+        Signal=[Signal P];
+        Class=[Class;1];
+        sumsignal=sumsignal+mwi(P);
+        SL=sumsignal/(length(Signal)+1);    %Updating the Signal Level
+    else        %Classification as Noise
+        Noise=[Noise P];
+        Class=[Class;0];
+        sumnoise=sumnoise+mwi(P);
+        NL=sumnoise/(length(Noise)+1);  %Updating the Noise Level
+    end
+end
+%------------------------------------------------------------
+function [pnew]=conversion(S,FL2,pold,M,N,fs)
+
+% conversion - sets the fiducial points of the downsampled Signal on the
+% samplepoints of the original Signal
+% 
+%   [pnew]=conversion(S,FL2,pold,M,N,fs)
+%
+% INPUT
+%   S           Original ECG Signal
+%   FL2         Feature of Level 2 [1]
+%   pold        old fiducial points
+%   M           M downsampling rate
+%   N           filter order
+%   fs          sample rate
+%
+% OUTPUT
+%   pnew        new fiducial points
+%
+
+Signaln=pold;    
+P=M;
+Q=1;
+FL2res=resample(FL2,P,Q);       %Resampling
+nans1=isnan(S);
+nans=find(nans1==1);
+S(nans)=mean(S);    %Replaces NaNs in Signal
+for i=1:length(Signaln)
+    Signaln1(i)=Signaln(i)+(M-1)*(Signaln(i)-1);    
+end
+%------------------- Sets the fiducial points on the maximum of FL2
+Signaln2=Signaln1;  
+Signaln2=Signaln2';     
+int=2*M;    %Window length for the new fiducial point
+range=1:length(FL2res);
+for i=1:length(Signaln2)
+    start=Signaln2(i)-int/2;
+    if start<1
+        start=1;
+    end
+    stop=Signaln2(i)+int/2;
+    if stop>length(FL2res)
+        stop=length(FL2res);
+    end
+    intervall=start:stop;       %interval
+    FL2int=FL2res(intervall);
+    pkt=find(FL2int==max(FL2int));  %Setting point on maximum of FL2
+    if length(pkt)==0   % if pkt=[];
+        pkt=Signaln2(i)-start;
+    else
+        pkt=pkt(1); 
+    end
+    delay=N/2+M;
+    Signaln3(i)=pkt+Signaln2(i)-int/2-delay;    %fiducial points according to FL2
+end
+%Sets the fiducial points on the maximum or minimum
+%of the signal
+Bw=5.6;   
+Bwn=1/(fs/2)*Bw;
+Wn=[Bwn 5*Bwn];
+N1=32;
+b=fir1(N1,Wn,'bandpass');
+Sf=filtfilt(b,1,S);     %Filtered Signal with bandwidth 5.6-28 Hz
+beg=round(1.5*M);
+fin=1*M;
+for i=1:length(Signaln3)
+    start=Signaln3(i)-beg;
+    if start<1
+        start=1;
+    end
+    stop=Signaln3(i)+fin;
+    if stop>length(Sf)
+        stop=length(Sf);
+    end
+    intervall=start:stop;   %Window for the new fiducial point
+    Sfint=abs(detrend(Sf(intervall),0));
+    pkt=find(Sfint==max(Sfint));    %Setting point on maximum of Sfint
+    if length(pkt)==0   %if pkt=[];
+        pkt=Signaln3(i)-start;
+    else
+        pkt=pkt(1); 
+    end
+    pkt=pkt(1);
+    Signaln4(i)=pkt+Signaln3(i)-beg-1;
+end
+Signal=Signaln4';   %New fiducial points according to the original signal
+
+cutbeginning=find(Signal<N);    %Cutting out the first points because of initial transient of the filter in polyphase_imp
+fpointsb=Signal(cutbeginning);
+cutend=find(Signal>length(S)-N); %Cutting out the last points
+fpointse=Signal(cutend);
+pnew=setxor(Signal,[fpointsb;fpointse]);
+%-------------------------------------------
+function yn=history(ynm1,xn)
+
+% history - computes y[n]=(1-lambda)*x[n]+lambda*y[n-1]
+%
+%   yn=history(ynm1,xn)
+
+lambda=0.8; %forgetting factor
+yn=(1-lambda)*xn+lambda*ynm1;
+
+