Skip to content

Commit

Permalink
res_rtp_asterisk: Fix standard deviation calculation
Browse files Browse the repository at this point in the history 
For some input to the standard deviation algorithm extremely large,
and wrong numbers were being calculated.

This patch uses a new formula for correctly calculating both the
running mean and standard deviation for the given inputs.

ASTERISK-29364 #close

Change-Id: Ibc6e18be41c28bed3fde06d612607acc3fbd621f
  • Loading branch information
@kharwell @gtjoseph
kharwell authored and gtjoseph committed Apr 1, 2021
1 parent c4a376a commit 0fc906a
Showing 1 changed file with 70 additions and 102 deletions.
172 changes: 70 additions & 102 deletions res/res_rtp_asterisk.c
View file
Expand Up @@ -384,7 +384,7 @@ struct ast_rtp {
unsigned int txcount; /*!< How many packets have we sent? */
unsigned int txoctetcount; /*!< How many octets have we sent? (txcount*160)*/
unsigned int cycles; /*!< Shifted count of sequence number cycles */
double rxjitter; /*!< Interarrival jitter at the moment in seconds */
double rxjitter; /*!< Interarrival jitter at the moment in seconds to be reported */
double rxtransit; /*!< Relative transit time for previous packet */
struct ast_format *lasttxformat;
struct ast_format *lastrxformat;
Expand Down Expand Up @@ -511,34 +511,36 @@ struct ast_rtcp {
unsigned int reported_jitter; /*!< The contents of their last jitter entry in the RR */
unsigned int reported_lost; /*!< Reported lost packets in their RR */

double reported_maxjitter;
double reported_minjitter;
double reported_normdev_jitter;
double reported_stdev_jitter;
unsigned int reported_jitter_count;

double reported_maxlost;
double reported_minlost;
double reported_normdev_lost;
double reported_stdev_lost;

double rxlost;
double maxrxlost;
double minrxlost;
double normdev_rxlost;
double stdev_rxlost;
unsigned int rxlost_count;

double maxrxjitter;
double minrxjitter;
double normdev_rxjitter;
double stdev_rxjitter;
unsigned int rxjitter_count;
double maxrtt;
double minrtt;
double normdevrtt;
double stdevrtt;
unsigned int rtt_count;
double reported_maxjitter; /*!< Maximum reported interarrival jitter */
double reported_minjitter; /*!< Minimum reported interarrival jitter */
double reported_normdev_jitter; /*!< Mean of reported interarrival jitter */
double reported_stdev_jitter; /*!< Standard deviation of reported interarrival jitter */
unsigned int reported_jitter_count; /*!< Reported interarrival jitter count */

double reported_maxlost; /*!< Maximum reported packets lost */
double reported_minlost; /*!< Minimum reported packets lost */
double reported_normdev_lost; /*!< Mean of reported packets lost */
double reported_stdev_lost; /*!< Standard deviation of reported packets lost */
unsigned int reported_lost_count; /*!< Reported packets lost count */

double rxlost; /*!< Calculated number of lost packets since last report */
double maxrxlost; /*!< Maximum calculated lost number of packets between reports */
double minrxlost; /*!< Minimum calculated lost number of packets between reports */
double normdev_rxlost; /*!< Mean of calculated lost packets between reports */
double stdev_rxlost; /*!< Standard deviation of calculated lost packets between reports */
unsigned int rxlost_count; /*!< Calculated lost packets sample count */

double maxrxjitter; /*!< Maximum of calculated interarrival jitter */
double minrxjitter; /*!< Minimum of calculated interarrival jitter */
double normdev_rxjitter; /*!< Mean of calculated interarrival jitter */
double stdev_rxjitter; /*!< Standard deviation of calculated interarrival jitter */
unsigned int rxjitter_count; /*!< Calculated interarrival jitter count */

double maxrtt; /*!< Maximum of calculated round trip time */
double minrtt; /*!< Minimum of calculated round trip time */
double normdevrtt; /*!< Mean of calculated round trip time */
double stdevrtt; /*!< Standard deviation of calculated round trip time */
unsigned int rtt_count; /*!< Calculated round trip time count */

/* VP8: sequence number for the RTCP FIR FCI */
int firseq;
Expand Down Expand Up @@ -3317,49 +3319,32 @@ static unsigned int ast_rtcp_calc_interval(struct ast_rtp *rtp)
return interval;
}

/*! \brief Calculate normal deviation */
static double normdev_compute(double normdev, double sample, unsigned int sample_count)
static void calc_mean_and_standard_deviation(double new_sample, double *mean, double *std_dev, unsigned int *count)
{
normdev = normdev * sample_count + sample;
sample_count++;
double delta1;
double delta2;

/*
It's possible the sample_count hits the maximum value and back to 0.
Set to 1 to prevent the divide by zero crash if the sample_count is 0.
*/
if (sample_count == 0) {
sample_count = 1;
}

return normdev / sample_count;
}
/* First convert the standard deviation back into a sum of squares. */
double last_sum_of_squares = (*std_dev) * (*std_dev) * (*count ?: 1);

static double stddev_compute(double stddev, double sample, double normdev, double normdev_curent, unsigned int sample_count)
{
/*
for the formula check http://www.cs.umd.edu/~austinjp/constSD.pdf
return sqrt( (sample_count*pow(stddev,2) + sample_count*pow((sample-normdev)/(sample_count+1),2) + pow(sample-normdev_curent,2)) / (sample_count+1));
we can compute the sigma^2 and that way we would have to do the sqrt only 1 time at the end and would save another pow 2 compute
optimized formula
*/
#define SQUARE(x) ((x) * (x))

stddev = sample_count * stddev;
sample_count++;
if (++(*count) == 0) {
/* Avoid potential divide by zero on an overflow */
*count = 1;
}

/*
It's possible the sample_count hits the maximum value and back to 0.
Set to 1 to prevent the divide by zero crash if the sample_count is 0.
* Below is an implementation of Welford's online algorithm [1] for calculating
* mean and variance in a single pass.
*
* [1] https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
*/
if (sample_count == 0) {
sample_count = 1;
}

return stddev +
( sample_count * SQUARE( (sample - normdev) / sample_count ) ) +
( SQUARE(sample - normdev_curent) / sample_count );
delta1 = new_sample - *mean;
*mean += (delta1 / *count);
delta2 = new_sample - *mean;

#undef SQUARE
/* Now calculate the new variance, and subsequent standard deviation */
*std_dev = sqrt((last_sum_of_squares + (delta1 * delta2)) / *count);
}

static int create_new_socket(const char *type, int af)
Expand Down Expand Up @@ -4434,7 +4419,6 @@ static void calculate_lost_packet_statistics(struct ast_rtp *rtp,
unsigned int expected_packets;
unsigned int expected_interval;
unsigned int received_interval;
double rxlost_current;
int lost_interval;

/* Compute statistics */
Expand Down Expand Up @@ -4464,6 +4448,13 @@ static void calculate_lost_packet_statistics(struct ast_rtp *rtp,
/* Update RTCP statistics */
rtp->rtcp->received_prior = rtp->rxcount;
rtp->rtcp->expected_prior = expected_packets;

/*
* While rxlost represents the number of packets lost since the last report was sent, for
* the calculations below it should be thought of as a single sample. Thus min/max are the
* lowest/highest sample value seen, and the mean is the average number of packets lost
* between each report. As such rxlost_count only needs to be incremented per report.
*/
if (lost_interval <= 0) {
rtp->rtcp->rxlost = 0;
} else {
Expand All @@ -4478,16 +4469,9 @@ static void calculate_lost_packet_statistics(struct ast_rtp *rtp,
if (lost_interval > rtp->rtcp->maxrxlost) {
rtp->rtcp->maxrxlost = rtp->rtcp->rxlost;
}
rxlost_current = normdev_compute(rtp->rtcp->normdev_rxlost,
rtp->rtcp->rxlost,
rtp->rtcp->rxlost_count);
rtp->rtcp->stdev_rxlost = stddev_compute(rtp->rtcp->stdev_rxlost,
rtp->rtcp->rxlost,
rtp->rtcp->normdev_rxlost,
rxlost_current,
rtp->rtcp->rxlost_count);
rtp->rtcp->normdev_rxlost = rxlost_current;
rtp->rtcp->rxlost_count++;

calc_mean_and_standard_deviation(rtp->rtcp->rxlost, &rtp->rtcp->normdev_rxlost,
&rtp->rtcp->stdev_rxlost, &rtp->rtcp->rxlost_count);
}

static int ast_rtcp_generate_report(struct ast_rtp_instance *instance, unsigned char *rtcpheader,
Expand Down Expand Up @@ -5423,7 +5407,6 @@ static void calc_rxstamp(struct timeval *tv, struct ast_rtp *rtp, unsigned int t
double prog;
int rate = ast_rtp_get_rate(rtp->f.subclass.format);

double normdev_rxjitter_current;
if ((!rtp->rxcore.tv_sec && !rtp->rxcore.tv_usec) || mark) {
gettimeofday(&rtp->rxcore, NULL);
rtp->drxcore = (double) rtp->rxcore.tv_sec + (double) rtp->rxcore.tv_usec / 1000000;
Expand Down Expand Up @@ -5458,11 +5441,8 @@ static void calc_rxstamp(struct timeval *tv, struct ast_rtp *rtp, unsigned int t
if (rtp->rtcp && rtp->rxjitter < rtp->rtcp->minrxjitter)
rtp->rtcp->minrxjitter = rtp->rxjitter;

normdev_rxjitter_current = normdev_compute(rtp->rtcp->normdev_rxjitter,rtp->rxjitter,rtp->rtcp->rxjitter_count);
rtp->rtcp->stdev_rxjitter = stddev_compute(rtp->rtcp->stdev_rxjitter,rtp->rxjitter,rtp->rtcp->normdev_rxjitter,normdev_rxjitter_current,rtp->rtcp->rxjitter_count);

rtp->rtcp->normdev_rxjitter = normdev_rxjitter_current;
rtp->rtcp->rxjitter_count++;
calc_mean_and_standard_deviation(rtp->rxjitter, &rtp->rtcp->normdev_rxjitter,
&rtp->rtcp->stdev_rxjitter, &rtp->rtcp->rxjitter_count);
}
}

Expand Down Expand Up @@ -5778,7 +5758,6 @@ static int update_rtt_stats(struct ast_rtp *rtp, unsigned int lsr, unsigned int
unsigned int rtt_lsw;
unsigned int lsr_a;
unsigned int rtt;
double normdevrtt_current;

gettimeofday(&now, NULL);
timeval2ntp(now, &msw, &lsw);
Expand Down Expand Up @@ -5815,16 +5794,8 @@ static int update_rtt_stats(struct ast_rtp *rtp, unsigned int lsr, unsigned int
rtp->rtcp->maxrtt = rtp->rtcp->rtt;
}

normdevrtt_current = normdev_compute(rtp->rtcp->normdevrtt,
rtp->rtcp->rtt,
rtp->rtcp->rtt_count);
rtp->rtcp->stdevrtt = stddev_compute(rtp->rtcp->stdevrtt,
rtp->rtcp->rtt,
rtp->rtcp->normdevrtt,
normdevrtt_current,
rtp->rtcp->rtt_count);
rtp->rtcp->normdevrtt = normdevrtt_current;
rtp->rtcp->rtt_count++;
calc_mean_and_standard_deviation(rtp->rtcp->rtt, &rtp->rtcp->normdevrtt,
&rtp->rtcp->stdevrtt, &rtp->rtcp->rtt_count);

return 0;
}
Expand All @@ -5836,7 +5807,6 @@ static int update_rtt_stats(struct ast_rtp *rtp, unsigned int lsr, unsigned int
static void update_jitter_stats(struct ast_rtp *rtp, unsigned int ia_jitter)
{
double reported_jitter;
double reported_normdev_jitter_current;

rtp->rtcp->reported_jitter = ia_jitter;
reported_jitter = (double) rtp->rtcp->reported_jitter;
Expand All @@ -5849,9 +5819,9 @@ static void update_jitter_stats(struct ast_rtp *rtp, unsigned int ia_jitter)
if (reported_jitter > rtp->rtcp->reported_maxjitter) {
rtp->rtcp->reported_maxjitter = reported_jitter;
}
reported_normdev_jitter_current = normdev_compute(rtp->rtcp->reported_normdev_jitter, reported_jitter, rtp->rtcp->reported_jitter_count);
rtp->rtcp->reported_stdev_jitter = stddev_compute(rtp->rtcp->reported_stdev_jitter, reported_jitter, rtp->rtcp->reported_normdev_jitter, reported_normdev_jitter_current, rtp->rtcp->reported_jitter_count);
rtp->rtcp->reported_normdev_jitter = reported_normdev_jitter_current;

calc_mean_and_standard_deviation(reported_jitter, &rtp->rtcp->reported_normdev_jitter,
&rtp->rtcp->reported_stdev_jitter, &rtp->rtcp->reported_jitter_count);
}

/*!
Expand All @@ -5861,11 +5831,10 @@ static void update_jitter_stats(struct ast_rtp *rtp, unsigned int ia_jitter)
static void update_lost_stats(struct ast_rtp *rtp, unsigned int lost_packets)
{
double reported_lost;
double reported_normdev_lost_current;

rtp->rtcp->reported_lost = lost_packets;
reported_lost = (double)rtp->rtcp->reported_lost;
if (rtp->rtcp->reported_jitter_count == 0) {
if (rtp->rtcp->reported_lost_count == 0) {
rtp->rtcp->reported_minlost = reported_lost;
}
if (reported_lost < rtp->rtcp->reported_minlost) {
Expand All @@ -5874,9 +5843,9 @@ static void update_lost_stats(struct ast_rtp *rtp, unsigned int lost_packets)
if (reported_lost > rtp->rtcp->reported_maxlost) {
rtp->rtcp->reported_maxlost = reported_lost;
}
reported_normdev_lost_current = normdev_compute(rtp->rtcp->reported_normdev_lost, reported_lost, rtp->rtcp->reported_jitter_count);
rtp->rtcp->reported_stdev_lost = stddev_compute(rtp->rtcp->reported_stdev_lost, reported_lost, rtp->rtcp->reported_normdev_lost, reported_normdev_lost_current, rtp->rtcp->reported_jitter_count);
rtp->rtcp->reported_normdev_lost = reported_normdev_lost_current;

calc_mean_and_standard_deviation(reported_lost, &rtp->rtcp->reported_normdev_lost,
&rtp->rtcp->reported_stdev_lost, &rtp->rtcp->reported_lost_count);
}

/*! \pre instance is locked */
Expand Down Expand Up @@ -6449,7 +6418,6 @@ static struct ast_frame *ast_rtcp_interpret(struct ast_rtp_instance *instance, s
}
update_jitter_stats(rtp, report_block->ia_jitter);
update_lost_stats(rtp, report_block->lost_count.packets);
rtp->rtcp->reported_jitter_count++;

if (rtcp_debug_test_addr(addr)) {
ast_verbose(" Fraction lost: %d\n", report_block->lost_count.fraction);
Expand Down

0 comments on commit 0fc906a

Please sign in to comment.