Remove outliers from HRV signal

I tried many ways of removing outliers from my signal, but it seems specific to me, so I opened this topic.The thing is, when I try to remove the artifacts from my ECG signal either by removing the values from file or from that specific vector, the outliers appear on the other place on signal. I need the "clean" signal without these outliers in order to get the right values of heart rate variability (HRV) parameters. My code is below, and the current signal is in attachment.

    clc; 
    clear all;
    close all;
    ecg=load('ecg.txt');
    f_s=250;
    N=length(ecg);
    t=[0:N-1]/f_s; %time period(total sample/Fs )

    % Peak detection algorithm 
    hh=ecg_smooth;
     j=[];           
    time=0;          
    th=0.3*max(hh);  %thresold setting at 30 percent of maximum value

    for i=2:N1-1 % length selected for comparison  
        % deopping first ie i=1:N-1  point because hh(1-1) 
       % in the next line  will be zero which is not appreciable in matlab 
        if((hh(i)>hh(i+1))&&(hh(i)>hh(i-1))&&(hh(i)>th))  
    % condition, i should be> then previous(i-1),next(i+1),threshold point;
            j(i)=hh(i);                                   
    %if condition satisfy store hh(i)in place of j(i)value which is initially 0;

            time(i)=[i-1]/250;           %position stored where peak value met;              

        end
    end
     j(j==0)=[];               % neglect all zeros from array;
     time(time==0)=[];     % neglect all zeros from array;
    m=(time)';               % converting rows in column;
    k=length(m);
    figure;
    plot(t,hh);            %x-axis time, y-smooth signal value;
    hold on;                 % hold the plot and wait for next instruction;
    plot(time,j,'*r'); title('PEAK POINTS DETECTED IN ECG SIGNAL')   

    %x-axis time, y-axis peak value,r=marker;
    xlabel('time')
    ylabel('amplitude')
    hold off                
    % to remove unwanted zeros from variable j and time 
    rr2=m(2:k);     %second array from 2nd to last point;
    rr1=m(1:k-1);   %first array from 1st to 2nd last point;
    % rr2 & rr1 is of equal length now;
    rr3=rr2-rr1;

    time1 = time(2:end);

    ix = { 1:400, 401:900, 901:1300, 1301:numel(time1) };
    colors = 'rymg';
    figure; cla;
    for i=1:4
        plot(time1(ix{i}),rr3(ix{i}),'Color',colors(i)); hold on
    end

    xlabel('Vreme [s]'),ylabel('Interval [ms]'), title('Difference between R-R intervals')
    end

[Image Analyst]

It seems weird that if you remove outliers from certain locations (elements) of your signal that they would magically appear at other locations in your signal.

I haven't run your code yet but would have liked to see plots of the signal, before and after outlier removal/moving, where you show how the outliers somehow mystically get transported to new locations. Maybe you can attach those screenshots before others open this post, to help them understand.

[miki90]

These are plots which were problematic in the other signal which file was too big to be attached. It turned out that I only needed to remove a bit greater intervals to remove the outliers. That, unfortunately, didn't resolve my problem of getting the right results when computing RMSSD and SDNN paramters.