[Asterisk-code-review] res_rtp_asterisk: Fix standard deviation calculation (asterisk[18])
George Joseph
asteriskteam at digium.com
Thu Apr 1 08:43:09 CDT 2021
George Joseph has submitted this change. ( https://gerrit.asterisk.org/c/asterisk/+/15673 )
Change subject: res_rtp_asterisk: Fix standard deviation calculation
......................................................................
res_rtp_asterisk: Fix standard deviation calculation
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
---
M res/res_rtp_asterisk.c
1 file changed, 67 insertions(+), 99 deletions(-)
Approvals:
Joshua Colp: Looks good to me, but someone else must approve
George Joseph: Looks good to me, approved; Approved for Submit
diff --git a/res/res_rtp_asterisk.c b/res/res_rtp_asterisk.c
index 4fd803b..607fd91 100644
--- a/res/res_rtp_asterisk.c
+++ b/res/res_rtp_asterisk.c
@@ -384,7 +384,7 @@
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;
@@ -511,34 +511,36 @@
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_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;
- double reported_minlost;
- double reported_normdev_lost;
- double reported_stdev_lost;
+ 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;
- double maxrxlost;
- double minrxlost;
- double normdev_rxlost;
- double stdev_rxlost;
- unsigned int rxlost_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;
- 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 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;
@@ -3317,49 +3319,32 @@
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;
+ /* First convert the standard deviation back into a sum of squares. */
+ double last_sum_of_squares = (*std_dev) * (*std_dev) * (*count ?: 1);
+
+ if (++(*count) == 0) {
+ /* Avoid potential divide by zero on an overflow */
+ *count = 1;
}
- return normdev / sample_count;
-}
-
-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++;
-
/*
- 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)
@@ -4434,7 +4419,6 @@
unsigned int expected_packets;
unsigned int expected_interval;
unsigned int received_interval;
- double rxlost_current;
int lost_interval;
/* Compute statistics */
@@ -4464,6 +4448,13 @@
/* 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 {
@@ -4478,16 +4469,9 @@
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,
@@ -5423,7 +5407,6 @@
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;
@@ -5458,11 +5441,8 @@
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);
}
}
@@ -5778,7 +5758,6 @@
unsigned int rtt_lsw;
unsigned int lsr_a;
unsigned int rtt;
- double normdevrtt_current;
gettimeofday(&now, NULL);
timeval2ntp(now, &msw, &lsw);
@@ -5815,16 +5794,8 @@
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;
}
@@ -5836,7 +5807,6 @@
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;
@@ -5849,9 +5819,9 @@
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);
}
/*!
@@ -5861,11 +5831,10 @@
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) {
@@ -5874,9 +5843,9 @@
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 */
@@ -6449,7 +6418,6 @@
}
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);
--
To view, visit https://gerrit.asterisk.org/c/asterisk/+/15673
To unsubscribe, or for help writing mail filters, visit https://gerrit.asterisk.org/settings
Gerrit-Project: asterisk
Gerrit-Branch: 18
Gerrit-Change-Id: Ibc6e18be41c28bed3fde06d612607acc3fbd621f
Gerrit-Change-Number: 15673
Gerrit-PatchSet: 5
Gerrit-Owner: Kevin Harwell <kharwell at digium.com>
Gerrit-Reviewer: Friendly Automation
Gerrit-Reviewer: George Joseph <gjoseph at digium.com>
Gerrit-Reviewer: Joshua Colp <jcolp at sangoma.com>
Gerrit-MessageType: merged
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.digium.com/pipermail/asterisk-code-review/attachments/20210401/75780719/attachment-0001.html>
More information about the asterisk-code-review
mailing list