D  Reference delay computation algorithm  p. 345

function ref_jb(channel,jb_adaptation_lookback,delay_delta_max,target_loss)
% channel         = file name of the channel
% lookback        = look back factor when estimating the max jitter
%                   buffer level [number of frames]
% delay_delta_max = max timescaling related modification (%) of the
%                   delay
% target_loss     = target late loss (%)
% example syntax:
% ref_jb('channel_1.dat',200,15,0.5);
framelength = 20;
% this value sets the speech data in each RTP packet to 20 ms. For 2 speech
% frames/RTP packet the value would be 40 ms. 
jitter_est_window=50; 
% Sets the jitter estimation window in number of frames
delay_delta_max_ms = framelength*delay_delta_max*0.01;
% Sets the maximum allowed time scaling
tscale = 1; 
% Scale factor of delay data
% In this case the files are assumend to be ascii files with one delay
% entry per line, the entries are in ms, a negative value denotes
% a packet loss.
x = load(channel);
x =x';
% remove packet losses
% remove inital startup empty frames
ix = find(x > 0);
x(1:ix(1)-1) = x(ix(1));
% remove packet losses (replace with nearby delay values)
ix = find(x < 0);
packet_loss = length(ix)/length(x)*100;
for n=1:length(ix)
    if (ix(n) > 1)
        x(ix(n)) = x(ix(n)-1);
    end;
end;
% convert timescale to ms
x = x*tscale;
L = length(x);
T = 1:L;
% estimate min and max TX delay, estimate a delta_delay
for n=1:L
    ix = [max(1,n-jitter_est_window):n];
    max_delay(n) = max(x(ix));
    min_delay(n) = min(x(ix));   
    delta_delay(n) = max_delay(n)-min_delay(n);
end
% compute the target max jitter buffer level with some slow adaptation
% downwards, just to mimick how a jitter buffer might behave
for n=1:L
    ix = [max(1,n-jb_adaptation_lookback):n];
    jb(n) = max(delta_delay(ix));
    % The timescaling is not allowed to adjust the jitterbuffer target max level
    % too fast.
    if n == 1
        jb_ = jb(n);
    end  
    delta = abs(jb_-jb(n));  
    if delta < delay_delta_max_ms;
        jb_ = jb(n);
    else
        if (jb(n) < jb_)
            jb_ = jb_-delay_delta_max_ms;
        else
            jb_ = jb_+delay_delta_max_ms;
        end
        jb(n) = jb_;
    end    
    % jitter buffer target max level can only assume an integer number of frames
    jbq(n) = ceil(jb(n)/framelength)*framelength;
    % compute estimated delay
    del(n) = jbq(n)+min_delay(n);    
end
if target_loss > 0
    % decrease the max jitter buffer leve until a target late loss has been
    % reached.
    late_loss = length(find(del < x))/L*100.0;
    jbq_save = jbq; % as the max level is increased until the late loss > target one
    % must be able to revert back to the previous data
    while late_loss < target_loss 
        jbq_save = jbq;
        jbq = min(max(jbq)-framelength,jbq);
        del = jbq+min_delay;
        late_loss = length(find(del < x))/L*100.0;
    end 
    jbq = jbq_save;
    del = jbq+min_delay;
end
jdel = max(0,del-x);
%Calculate and plot the CDF of the reference buffer.
figure(1);plot(T,jbq,T,del,T,x);
[n,x] = hist(jdel,140); y = cumsum(n);y = y/max(y)*100;
figure(2);plot(x,y);axis([0 200 0 100]);ylabel('%');xlabel('ms');title('CDF of packet delay in JB');