function challenge(recordName)
%POPRAWIAM 29.05.2014

%--------------------------------------------------------------------------
%-[1] Konfiguracja progdispramu -----------------------------------------------
%--------------------------------------------------------------------------
%tic
	[~,config] = wfdbloadlib;                           % inicjalizacja pakietu wfdb
    if(config.inOctave)
        crash_dumps_octave_core(0);
    end
              
	annName = 'qrs';                                    % ustalenie nazwy anotacji (wymagane przez warunki Challenge)
	
	[siginfo, fs] = wfdbdesc2(recordName);   
                  % wczytanie danych informacyjnych o nagraniu pod zmienną siginfo (dziwna struktura o formacie wfdb)
	description = squeeze(struct2cell(siginfo));        % przekonwertowanie danych informacyjnych na czytelny format tablicy "komórkowej" matlaba
	descriptionNamesTable = description(5,:);           % wyciągnięcie samych nazw sygnałów na wyjściu mamy jednowymiarową tablicę ze stringami zawierającymi nazwy kanałów
    if(~isnan(fs))
        sfreq = fs(1);
    else
      sfreq=sscanf(siginfo(1).SamplingFrequency,'%f');
    end
    %%LOADING MATLAB
    
    if(config.inOctave)
        stringInput =sprintf('wfdb2mat -r %s -H',recordName);
    	system(stringInput);
        load(strcat(recordName, 'm.mat'));
        signal=val(1:end,:)'; %#ok<NODEF>
        tm=(1:size(signal,1))./sfreq; %#ok<NASGU>
        for ii=1:size(signal,2)
            gain=sscanf(siginfo(ii).Gain,'%f');
            base=sscanf(siginfo(ii).Baseline,'%f');
            signal(:,ii)=(signal(:,ii)-base)./gain;
        end
    else
    	[tm,signal] = rdsamp(recordName);                   %#ok<ASGLU> % wczytanie sygnałów dla aktualnie analizowanego przypadku 
                                                        % tm - podstawa czasowa sygnału
                                                        % signal - wszystkie sygnały istniejące dla danego przypadku. Tablica N x m, gdzie N to długość nagrania, a m - ilość dostępnych sygnałów
   
    end
	

   % eeeLista=[];                    % utworzenie pustej tablicy dla indeksow sygnalow EOG, EEG i EMG (o ile istnieja)
%--------------------------------------------------------------------------  
%-[2] Rozpoznanie, któy sygnał jest który ---------------------------------
%--------------------------------------------------------------------------
    ecg_ind = 1;
    ecg_ind = [ecg_ind; get_index(descriptionNamesTable,'ECG'); get_index(descriptionNamesTable,'Ecg'); get_index(descriptionNamesTable,'ecg')];
    ecg_ind = [ecg_ind; get_index(descriptionNamesTable,'EKG'); get_index(descriptionNamesTable,'Ekg'); get_index(descriptionNamesTable,'ekg')];
    ecg_ind = unique(ecg_ind);

    bp_ind = [get_index(descriptionNamesTable,'BP'); get_index(descriptionNamesTable,'Bp'); get_index(descriptionNamesTable,'bp')];
    bp_ind = [bp_ind; get_index(descriptionNamesTable,'ART'); get_index(descriptionNamesTable,'art')];
    bp_ind = unique(bp_ind);
    
	
	
%     eeg_ind = [get_index(descriptionNamesTable,'EEG'); get_index(descriptionNamesTable,'Eeg'); get_index(descriptionNamesTable,'eeg')];
%     eeg_ind = unique(eeg_ind);
%     if(length(eeg_ind)>0) 
%         eeeLista = [eeeLista; eeg_ind];        %kasiamateusz - przypisanie indeksu sygnalu EEG do tablicy 
%     end
%     
%     eog_ind = [get_index(descriptionNamesTable,'EOG'); get_index(descriptionNamesTable,'Eog'); get_index(descriptionNamesTable,'eog')];
%     eog_ind = [eog_ind; get_index(descriptionNamesTable,'eye'); get_index(descriptionNamesTable,'Eye')];
%     eog_ind = unique(eog_ind);
%     if(length(eog_ind)>0)
%         eeeLista = [eeeLista; eog_ind];        %kasiamateusz - przypisanie indeksu sygnalu EOG do tablicy 
%     end
% 
%     emg_ind = [get_index(descriptionNamesTable,'EMG'); get_index(descriptionNamesTable,'Emg'); get_index(descriptionNamesTable,'emg'); get_index(descriptionNamesTable,'leg')];
%     emg_ind = unique(emg_ind);
%     if(length(emg_ind)>0) 
%         eeeLista = [eeeLista; emg_ind];   %kasiamateusz - przypisanie indeksu sygnalu EMG do tablicy 
%     end
%     eeeLista = eeeLista';
 
    

%--------------------------------------------------------------------------
%-[3] Analiza sygnałów EKG ------------------------------------------------
%--------------------------------------------------------------------------
    %granice dla sprawdzania:
    graniceRSslopeEKG = [   0.017, 0.036;... 
                            0.026, 0.036;...
                            0.035, 0.046;...
                            0.017, 0.046;...
                            0.026, 0.046;...  
                            0.017, 0.056;...
                            0.026, 0.056;... 
                            0.035, 0.056;...
                            0.017, 0.15;... 
                            0.035, 0.15];
    graniceRSslopeEKG = floor(graniceRSslopeEKG .* sfreq);    
    graniceAmpEKG= [0.03, 0.97];
    % teraz ustawiam parametry do wyszukiwania pikow R za pomoca algorytmu RS slope
    %dolnaLOW= floor( 0.017 *sfreq);                                     %  dolnaLOW   - lewy zakres szukania minimalnej dlugosci spadku na 17 ms
    %dolnaHIGH= floor(0.033 *sfreq);                                     %  dolnaHIGH  - prawy zakres szukania minimalnej dlugosci spadku na 36 ms
    %gornaLOW= floor(0.036 *sfreq);                                      %  gornaLOW - lewy zakres szukania maksymalnej dlugosci spadku na 40 ms
    %gornaHIGH= floor(0.055 *sfreq);                                     %  gornaHIGH  - prawy zakres szukania maksymalnej dlugosci spadku na 55 ms
    %stepLn = ceil(0.006*sfreq);                                        %  stepLn - rozdzielczosc szukania dobrych zakresow dlugosci spadkow na 2 ms
    lowAmpBest = 0.03;                                                        %#ok<NASGU> % bede szukal spadku RS, tylko wtedy, gdy amplituda piku R jest wyzsza od wartości 95% próbek w sygnale
    highAmpBest = 0.97;                                                     %#ok<NASGU>   % nie ma górnego ograniczenia na wartość amplitudy R (100% rozkładu amplitud próbek)
    window = round(0.250*sfreq);                                        % wybieram dlugosc okna dla ruchomej mediany na 250 ms
    quality=0;                                             %#ok<NASGU> 
    AnnSkeletonECG=[];                                             %#ok<NASGU> 
    dAnnSkeletonECG=[];                                             %#ok<NASGU> 
    ANNourREFstage1=[]; 
    prog1=0;                                             %#ok<NASGU> 
    prog2=0;                                             %#ok<NASGU> 
    lowLengthBest=floor( 0.017 *sfreq);                                             %#ok<NASGU> 
    highLengthBest=floor(0.055 *sfreq);                                             %#ok<NASGU> 
    fullSkelet =[];                                             %#ok<NASGU> 
    dobrySYGNAL=0;                                             %#ok<NASGU> 
    resultHistEKG = 'FAILED';                                             %#ok<NASGU> 

    
    if(~isempty(ecg_ind) )                                               % sprawdzam, czy występuje w ogole kanal EKG. Jesli tak - zaczynam analize 
        for currentEcg=1:length(ecg_ind)
             disp('EKG');
            dobrySYGNAL = movingMedian(signal(:,ecg_ind(currentEcg)),window);               % odejmowanie ruchomej mediany od sygnalu EKG za pomoca funkcji umieszczonej w pliku movingMedian.m
                                                                                %  dobrySygnal - sygnal EKG po odjeciu mediany
                                                                                % uruchamiam algorytm wyszukiwania spadkow RS
            [b, a] = butter(8,40/(sfreq/2));                                    % filtruje sygnał poniżej 40 Hz filtrem Butterwortha 8 rzędu
            dobrySYGNAL = filter(b,a,dobrySYGNAL);                                                                           
            [AnnSkeletonECG,lowLengthBest,highLengthBest, fullSkelet,lowAmpBest,highAmpBest] = rsSlopeDetectionEKG(dobrySYGNAL,sfreq, graniceRSslopeEKG, graniceAmpEKG,1); %#ok<NASGU,ASGLU>
                                                                                %  AnnSkeletonECG - znalezione anotacje EKG ( szkielet 1.0)
                                                                                %  lowLengthBest - dobrana minimalna dlugosc spadku RS (w probkach)
                                                                                %  highLengthBest - dobrana maksymalna dlugosc spadku RS (w probkach)
                                                                                %  fullSkelet - znalezione anotacje EKG wraz z dodatkowymi, nie odrzucownymi podczas czyszczenia smieciami                                          
            quality = signalQuality( fullSkelet,size(dobrySYGNAL,1),sfreq );
            if(quality==0)
                %ANN
            else
                stary = ANNourREFstage1;    
                nowy = AnnSkeletonECG;   
                polaczonySzkielet=joinSkeletons(stary,nowy,sfreq);              % lacze anotacje EKG (stary szkielet) z anotacjami z kanalu BP, uzyskując nowy szkielet (szkielet 2.0)
                ANNourREFstage1 =  polaczonySzkielet;                           %  ANNourREFstage2 - ostateczny szkielet 2.0 
            end 
        end
    end         
        
        

    
%--------------------------------------------------------------------------
%-[3] Analiza sygnału BP --------------------------------------------------
%--------------------------------------------------------------------------

    window = floor(0.35*sfreq);                                            % wybieram dlugosc okna dla ruchomej mediany na 350 ms                                                                           %  dobreBP - sygnal BP po odjeciu mediany i odwroceniu czasu
    dolnaLOW= floor( 0.090 *sfreq);                                 %  dolnaLOW   - lewy zakres szukania minimalnej dlugosci spadku na 90 ms
    dolnaHIGH= floor(0.090 *sfreq);                                 %  dolnaHIGH  - prawy zakres szukania minimalnej dlugosci spadku na 110 ms
    gornaLOW= floor(0.140 *sfreq);                                  %  usgornaLOW - lewy zakres szukania maksymalnej dlugosci spadku na 140 ms
    gornaHIGH= floor(0.160 *sfreq);                                 %  gornaHIGH  - prawy zakres szukania maksymalnej dlugosci spadku na 160 ms
    stepLn = round(0.02*sfreq);                                    %  stepLn - rozdzielczosc szukania dobrych zakresow dlugosci spadkow na 5 ms
    lowg = 0.10;                                                     % bede szukal spadku RS, tylko wtedy, gdy amplituda piku R jest wyzsza od wartości 90% próbek w sygnale
    highg = 0.90;                                                       % nie ma górnego ograniczenia na wartość amplitudy R (100% rozkładu amplitud próbek)
    offset=0;                                             %#ok<NASGU> 
    qualityBP=0;                                             %#ok<NASGU> 
    dobreBPmedian = 0;                                             %#ok<NASGU> 
    AnnSkeletonBP = [];                                             %#ok<NASGU> 
    dAnnSkeletonBP =  [];                                               %#ok<NASGU> % jesli nie bylo sygnalu BP -  szkielet 2.0 bedzie po prostu kopia szkieletu 1.0
    ANNourREFstage2=ANNourREFstage1;                                           
    resultHistBP = 'FAILED';                                               %#ok<NASGU> 

    if(~isempty(bp_ind) )                                                % sprawdzam, czy występuje w ogole kanal cisnienia BP. Jesli tak - zaczynam analize
         for currentBP=1:length(bp_ind)      
             
            disp('BP');
            dobreBPmedian = movingMean(signal(:,bp_ind(currentBP)),window);            % odejmowanie ruchomej mediany od sygnalu BP za pomoca funkcji umieszczonej w pliku movingMedian.m                   
            dobreBP = flipud(dobreBPmedian);                                % Odwracamy czas (-t), by patrzec na wzrosty. 
                                                             % uruchamiam algorytm wyszukiwania spadkow RS na sygnale BP
            
            [b, a] = butter(8,40/(sfreq/2));  
            dobreBP = filter(b,a,dobreBP);
            [annsBP,lowLengthBestBP,highLengthBestBP, fullSkeletBP] = rsSlopeDetection(dobreBP,sfreq, dolnaLOW, dolnaHIGH, gornaLOW, gornaHIGH, stepLn, lowg, highg,0); %#ok<ASGLU>
                                                                            %  annsBP - znalezione anotacje BP
                                                                            %  lowLengthBestBP - dobrana minimalna dlugosc spadku krzywej BP (w probkach)
                                                                            %  highLengthBestBP - dobrana maksymalna dlugosc spadku krzywej BP (w probkach)
                                                                            %  fullSkeletBP - znalezione anotacje pikow BP wraz z dodatkowymi, nie odrzucownymi podczas czyszczenia smieciami            
            annsBP_notmoved = size(dobreBP,1) - annsBP;                     % Odwracamy z powrotem czas dla anotacji (-t) --> (t),  mieć je ustawione w dobrym kierunku. 
           
            if(~isempty(ecg_ind))                                       % Jesli byl sygnal EKG, licze na jego podstawie srednie przesuniecie pomiedzy pikami R w sygnale EKG, a maksimami wychylen krzywej cisnienia
                if(size(ANNourREFstage1,1)<3)
                    offset = floor(0.28 * sfreq);
                else
                    offset = floor(bpOffsetDetection(ANNourREFstage1,annsBP_notmoved,sfreq));
                end
            else                                                            % Jesli nie bylo kanalu EKG - ustalam przeuniecie na 280 ms
                offset = floor(0.28 * sfreq);                                      %  offset - srednie przesuniecie pomiedzy pikami R w EKG a pikami krzywej BP
            end
            AnnSkeletonBP = annsBP_notmoved - round(offset);                % Przesuwam anotacje dla krzywej BP by otrzymac ostateczna liste anotacji                            
            AnnSkeletonBP = sort(AnnSkeletonBP);  
            
            qualityBP = signalQuality( fullSkeletBP,size(dobreBP,1),sfreq );
            if(qualityBP==0)
              %  AnnSkeletonBP=[];                                             %#ok<NASGU> 
            else
                stary = ANNourREFstage2;              
                nowy = AnnSkeletonBP; 
                polaczonySzkielet=joinSkeletons(stary,nowy,sfreq);              % lacze anotacje EKG (stary szkielet) z anotacjami z kanalu BP, uzyskując nowy szkielet (szkielet 2.0)
                ANNourREFstage2 =  polaczonySzkielet;                           %  ANNourREFstage2 - ostateczny szkielet 2.0 
            end    
         end
    end
           
%--------------------------------------------------------------------------
%-[3] Analiza GQRS --------------------------------------------------
%--------------------------------------------------------------------------
    ANNourREFstage3=ANNourREFstage2;  
    if(~isempty(ecg_ind))                                               % sprawdzam, czy występuje w ogole kanal EKG. Jesli tak - zaczynam analize 
        for currentEcg=1:length(ecg_ind)
            disp('GQRS');
         
            
            dobrySYGNAL = movingMedian(signal(:,ecg_ind(currentEcg)),window);               % odejmowanie ruchomej mediany od sygnalu EKG za pomoca funkcji umieszczonej w pliku movingMedian.m
            [b, a] = butter(8,40/(sfreq/2));                                    % filtruje sygnał poniżej 40 Hz filtrem Butterwortha 8 rzędu
            dobrySYGNAL = filter(b,a,dobrySYGNAL);  

            

            newrecordName = strcat(recordName,'new');
            mat2wfdb(dobrySYGNAL,newrecordName,sfreq,[],'mV');
            gqrs(newrecordName,[],[],ecg_ind(currentEcg),[],'qra');          %tutaj jest recordName a nie dobrySYGNAL, poniewaz gqrs przyjmuje stringa oznaczajacego nazwe pliku
            ANNfromGQRS = rdann(newrecordName,'qra');
            delete(strcat(newrecordName,'.dat'));
            delete(strcat(newrecordName,'.hea'));
            delete(strcat(newrecordName,'.qra'));
                        
            stary = ANNourREFstage2;    
            nowy = ANNfromGQRS;   
            polaczonySzkielet=joinSkeletons(stary,nowy,sfreq);              % lacze anotacje EKG (stary szkielet) z anotacjami z kanalu BP, uzyskując nowy szkielet (szkielet 2.0)
            ANNourREFstage3 =  polaczonySzkielet;                           %  ANNourREFstage2 - ostateczny szkielet 2.0 

        end
    end         
            

      

%--------------------------------------------------------------------------
%-[8] ZABEZPIECZNIE --------------------------------------
%--------------------------------------------------------------------------    
    if(length(ANNourREFstage3)<5)
        disp('siara');
        sygnalDoGQRS=0;
        
        allSignals = size(signal,2);
        
        lowAmpBest = 0.03;                                                        %#ok<NASGU> % bede szukal spadku RS, tylko wtedy, gdy amplituda piku R jest wyzsza od wartości 95% próbek w sygnale
        highAmpBest = 0.97;                                         %#ok<NASGU>                         
        window = round(0.250*sfreq);                                        % wybieram dlugosc okna dla ruchomej mediany na 250 ms
        quality=0;                                     %#ok<NASGU> 
        AnnSkeletonECG=[];                                     %#ok<NASGU> 
        dAnnSkeletonECG=[];                                     %#ok<NASGU> 
        ANNourREFstage1=[];                                     
        lowLengthBest=floor( 0.017 *sfreq);                                             %#ok<NASGU> 
        highLengthBest=floor(0.055 *sfreq);                                             %#ok<NASGU> 
        fullSkelet =[];                                     %#ok<NASGU> 
        dobrySYGNAL=0;                                     %#ok<NASGU> 

        for currentEcg=1:allSignals
            dobrySYGNAL = movingMedian(signal(:,currentEcg),window);               % odejmowanie ruchomej mediany od sygnalu EKG za pomoca funkcji umieszczonej w pliku movingMedian.m
                                                                               %  dobrySygnal - sygnal EKG po odjeciu mediany
                                                                                % uruchamiam algorytm wyszukiwania spadkow RS
            [b, a] = butter(8,40/(sfreq/2));                                    % filtruje sygnał poniżej 40 Hz filtrem Butterwortha 8 rzędu
            dobrySYGNAL = filter(b,a,dobrySYGNAL);          
            
            [AnnSkeletonECG,lowLengthBest,highLengthBest, fullSkelet,lowAmpBest,highAmpBest] = rsSlopeDetectionEKG(dobrySYGNAL,sfreq, graniceRSslopeEKG, graniceAmpEKG,1); %#ok<NASGU,ASGLU>

            quality = signalQuality( fullSkelet,size(dobrySYGNAL,1),sfreq );
            if(quality==0)
                %disp('ECG fuj');
            else
                %disp('ECG ok');
                if(sygnalDoGQRS==0)
                     sygnalDoGQRS=currentEcg;
                end
                stary = ANNourREFstage1;    
                nowy = AnnSkeletonECG;   
                polaczonySzkielet=joinSkeletons(stary,nowy,sfreq);              % lacze anotacje EKG (stary szkielet) z anotacjami z kanalu BP, uzyskując nowy szkielet (szkielet 2.0)
                ANNourREFstage1 =  polaczonySzkielet;                           %  ANNourREFstage2 - ostateczny szkielet 2.0 
            end 
        end
                 

         %-[3] Analiza sygnału BP --------------------------------------------------

        window = floor(0.35*sfreq);                                            % wybieram dlugosc okna dla ruchomej mediany na 350 ms                                                                           %  dobreBP - sygnal BP po odjeciu mediany i odwroceniu czasu
        dolnaLOW= floor( 0.090 *sfreq);                                 %  dolnaLOW   - lewy zakres szukania minimalnej dlugosci spadku na 90 ms
        dolnaHIGH= floor(0.090 *sfreq);                                 %  dolnaHIGH  - prawy zakres szukania minimalnej dlugosci spadku na 110 ms
        gornaLOW= floor(0.160 *sfreq);                                  %  usgornaLOW - lewy zakres szukania maksymalnej dlugosci spadku na 140 ms
        gornaHIGH= floor(0.160 *sfreq);                                 %  gornaHIGH  - prawy zakres szukania maksymalnej dlugosci spadku na 160 ms
        stepLn = round(0.02*sfreq);                                    %  stepLn - rozdzielczosc szukania dobrych zakresow dlugosci spadkow na 5 ms
        lowg = 0.10;                                                     % bede szukal spadku RS, tylko wtedy, gdy amplituda piku R jest wyzsza od wartości 90% próbek w sygnale
        highg = 0.90;                                                       % nie ma górnego ograniczenia na wartość amplitudy R (100% rozkładu amplitud próbek)
        offset=0;                                     %#ok<NASGU> 
        qualityBP=0;                                     %#ok<NASGU> 
        dobreBPmedian = 0;                                     %#ok<NASGU> 
        AnnSkeletonBP = [];                                     %#ok<NASGU> 
        dAnnSkeletonBP =  [];                                       %#ok<NASGU>                                            % jesli nie bylo sygnalu BP -  szkielet 2.0 bedzie po prostu kopia szkieletu 1.0
        ANNourREFstage2=ANNourREFstage1;
        resultHistBP = 'FAILED';                                      %#ok<NASGU> 

       for currentBP=1:allSignals    
            dobreBPmedian = movingMean(signal(:,currentBP),window);            % odejmowanie ruchomej mediany od sygnalu BP za pomoca funkcji umieszczonej w pliku movingMedian.m                   
            dobreBP = flipud(dobreBPmedian);                                % Odwracamy czas (-t), by patrzec na wzrosty. 
                                                             % uruchamiam algorytm wyszukiwania spadkow RS na sygnale BP

            [b, a] = butter(8,40/(sfreq/2));  
            dobreBP = filter(b,a,dobreBP);
            [annsBP,lowLengthBestBP,highLengthBestBP, fullSkeletBP] = rsSlopeDetection(dobreBP,sfreq, dolnaLOW, dolnaHIGH, gornaLOW, gornaHIGH, stepLn, lowg, highg,0); %#ok<ASGLU>
                                                                            %  annsBP - znalezione anotacje BP
                                                                            %  lowLengthBestBP - dobrana minimalna dlugosc spadku krzywej BP (w probkach)
                                                                            %  highLengthBestBP - dobrana maksymalna dlugosc spadku krzywej BP (w probkach)
                                                                            %  fullSkeletBP - znalezione anotacje pikow BP wraz z dodatkowymi, nie odrzucownymi podczas czyszczenia smieciami            
            annsBP_notmoved = size(dobreBP,1) - annsBP;                     % Odwracamy z powrotem czas dla anotacji (-t) --> (t),  mieć je ustawione w dobrym kierunku. 

            if(length(ecg_ind)>0 )                                       %#ok<ISMT> % Jesli byl sygnal EKG, licze na jego podstawie srednie przesuniecie pomiedzy pikami R w sygnale EKG, a maksimami wychylen krzywej cisnienia
                if(size(ANNourREFstage1,1)<3)
                    offset = floor(0.28 * sfreq);
                else
                    offset = floor(bpOffsetDetection(ANNourREFstage1,annsBP_notmoved,sfreq));
                end
            else                                                            % Jesli nie bylo kanalu EKG - ustalam przeuniecie na 280 ms
                offset = floor(0.28 * sfreq);                                      %  offset - srednie przesuniecie pomiedzy pikami R w EKG a pikami krzywej BP
            end
 

            qualityBP = signalQuality( fullSkeletBP,size(dobreBP,1),sfreq );
            if(qualityBP==0)
                %AnnSkeletonBP=[];     
                %disp('BP fuj');
            else
                %disp('BP ok');
                AnnSkeletonBP = annsBP_notmoved - round(offset);                % Przesuwam anotacje dla krzywej BP by otrzymac ostateczna liste anotacji                            
                AnnSkeletonBP = sort(AnnSkeletonBP); 
                stary = ANNourREFstage2;              
                nowy = AnnSkeletonBP; 
                polaczonySzkielet=joinSkeletons(stary,nowy,sfreq);              % lacze anotacje EKG (stary szkielet) z anotacjami z kanalu BP, uzyskując nowy szkielet (szkielet 2.0)
                ANNourREFstage2 =  sort(polaczonySzkielet);                           %  ANNourREFstage2 - ostateczny szkielet 2.0 
            end    
       end
        %GQRS FAILSAFE
        ANNourREFstage3=ANNourREFstage2;   %#ok<NASGU>

        
        

        dobrySYGNAL = movingMedian(signal(:,ecg_ind(currentEcg)),window);               % odejmowanie ruchomej mediany od sygnalu EKG za pomoca funkcji umieszczonej w pliku movingMedian.m
        [b, a] = butter(8,40/(sfreq/2));                                    % filtruje sygnał poniżej 40 Hz filtrem Butterwortha 8 rzędu
        dobrySYGNAL = filter(b,a,dobrySYGNAL);  



        newrecordName = strcat(recordName,'new');
        mat2wfdb(dobrySYGNAL,newrecordName,sfreq,[],'mV');
        gqrs(newrecordName,[],[],ecg_ind(currentEcg),[],'qra');          %tutaj jest recordName a nie dobrySYGNAL, poniewaz gqrs przyjmuje stringa oznaczajacego nazwe pliku
        ANNfromGQRS = rdann(newrecordName,'qra');
        delete(strcat(newrecordName,'.dat'));
        delete(strcat(newrecordName,'.hea'));
        delete(strcat(newrecordName,'.qra'));

        stary = ANNourREFstage2;    
        nowy = ANNfromGQRS;   
        polaczonySzkielet=joinSkeletons(stary,nowy,sfreq);              % lacze anotacje EKG (stary szkielet) z anotacjami z kanalu BP, uzyskując nowy szkielet (szkielet 2.0)
        ANNourREFstage3 =  polaczonySzkielet;                           %  ANNourREFstage2 - ostateczny szkielet 2.0 
 
    end
    
    
    
    
    
    ANNost = ANNourREFstage3;      
    ANNost = sort(ANNost);
    
    
    %ostatni FAILSAFE
    if(length(ANNost)<5)
        disp('ostatni FAILSAFE');
        %        signalFAILSAFE = movingMedian(signal(:,1),window); 
        %        annFAILSAFE = QRSParameters(signalFAILSAFE,0.1,0.9,sfreq*0.014,sfreq*0.08,0);
        %        ANNourREFstage3 = sort(annFAILSAFE(:,1));
        %        [ANNostAMP,indexesAMP] = sort(annFAILSAFE(:,2),'descend');
        %        ANNost = ANNourREFstage3(indexesAMP);
        %        ileDobrych = (length(tm)/sfreq)/0.9;
        %        if(length(ANNostAMP)>ileDobrych)
        %            ANNost = ANNost(1:ileDobrych);
        %            ANNost = sort(ANNost);     
        %        end      
       gqrs(recordName,[],[],ecg_ind(currentEcg),[],'qra');          %tutaj jest recordName a nie dobrySYGNAL, poniewaz gqrs przyjmuje stringa oznaczajacego nazwe pliku
       ANNfromGQRS = rdann(recordName,'qra');
       ANNost=ANNfromGQRS;
    end
    disp('writing ANNs');
    %ANNost = finalFantasy(ANNost(:,1),ANNost(:,1),sfreq);               %  funkcja uzupelniania i czysczenia ostatecznych anotacji
    dANNost = [sfreq; diff(ANNost)]; 
    ANNost(dANNost < 0.08*sfreq) = []; 
    ANNost(ANNost<=0)=[];
    dANNost = diff(ANNost);  %#ok<NASGU>
    ANNostWrite = floor(ANNost);     
                                        %  ANNost - ostateczny szkielet gotowy do ostatecznego wygladzania (post-processing)
                                                     %  dANNost - ostateczna lista interwalow RR dla ostatecznego szkieletu                                       
  
     wrann(recordName,annName,ANNostWrite,'N',0,0,0);                              % zapisz wyniki pracy (anotacje ostateczne ) do pliku wyjsciowego o rozszerzeniu annName

end


function ind = get_index(description,pattern)
    M = length(description);
    tmp_ind = strfind(description,pattern);
    ind = [];
    for m = 1:M
        if(~isempty(tmp_ind{m}))
            ind(end+1) = m; %#ok<AGROW>
        end
    end
end


