Olav's 2010-09-15 version

d0aa7efd · Jon Kåre Hellan · 5c672be5 · d0aa7efd
Commit d0aa7efd authored Apr 28, 2011 by Jon Kåre Hellan
Show whitespace changes
Inline Side-by-side

Showing with 327 additions and 29 deletions

qstream qstream +327 -29

No files found.
--- a/qstream
+++ b/qstream
 #!/usr/bin/perl

 # compute packet stream quality statistics for UDP/RTP from network or captured files(pcap)
+# Base Olav Kvittem
+# Mpeg support Odd Rune Mykkeltveit Lykkbo,  2009-2010
+# Mos values computation Gurvinder Singh, 2010-05-07

 # use Socket;
 use IO::Socket qw(:DEFAULT :crlf);
@@ -41,6 +44,8 @@ $usage="$0 [option]... [file...|ip/|:port]...
 -log log gaps longer than 1 (use with crude)
 -rtp look forinto rtp headers (default just udp)
 -mpeg decode mpeg-2 transport stream and look for sequence and mediatime
+-rtt  roundtrip  - round trip time in ms (for MOS)
+-codec codec - codec G.711, G.723.1, G.729A, iLBC (G.711)
 -pcap listen to all multicast groups with pcap at the same time
 -period seconds - period in seconds between stats output
 -nperiod number of periodes
@@ -58,11 +63,15 @@ $usage="$0 [option]... [file...|ip/|:port]...
 -debug  - print debugging info
 -help/h this message
 \n";
-@opts=('list', 'fullformat', 'net', 'rtp', 'mpeg', 'pcap', 'crude', 'bins=s', 'log', 'nperiod=s', 'period=s', 'last=s', 'packets=s', 'src=s', 'flow_key=s', 'flow_no=s', 'flow_min=s', 'nohead', 'dump=s', 'format=s', 'id=s', 'sum=s', 'ttl', 'log', 'verbose', 'v', 'help', 'h', 'debug');
+@opts=('list', 'fullformat', 'net', 'rtp', 'mpeg', 'pcap', 'crude', 'bins=s', 'log', 'nperiod=s', 'period=s', 'last=s', 'packets=s', 'src=s', 'flow_key=s', 'flow_no=s', 'flow_min=s', 'nohead', 'dump=s', 'format=s', 'id=s', 'sum=s', 'ttl', 'rtt=s', 'codec=s', 'log', 'verbose', 'v', 'help', 'h', 'debug');
 &NGetOpt(@opts) || die $usage;
 die $usage if $opt_h or $opt_help;
 $continous= ! $opt_packets;
 $opt_v=$opt_v or $opt_verbose;
+my $codec = $opt_codec || 'G.711';
+my $rtt = $opt_rtt || 10;		# ms round trip time of the connection
+
+
 my @streams=@ARGV;

 if ($opt_bins){
@@ -123,7 +132,7 @@ $SIG{KILL} = sub { $uninterrupted=0; # return if $nkill++ < 1;
 $SIG{ALRM} = sub { 
 #    die "No packets received" if $n_packets < 1;
 #    $uninterrupted=0;
-    &display_stats();
+    &display_stats if ! $endstream;
 #	die "End after alarm\n"; 
 #	return(0);
    $endstream=1;
@@ -854,14 +863,14 @@ sub display_stats{

 		    print  "      flu-    pcr_jitter(ms)             dwell_time(ms)   " if $opt_mpeg;

-		    print "        thrust(bps)            source\n";                      
+		    print "             load(bps)            source\n";                      
 		    printf "%-8s %6s %6s     %4s%4s %4s %4s %4s   %6s %4s %4s %4s %4s",
 		    "time","span", "setup",
 		    "numb ", "avg", "sdv", "min", "max",
 		    "numb", "avg", "sdv", "min", "max";
 		    printf "  %5s %5s %5s %5s %5s %5s %5s",
 		    "stamp", "dup", "late", "lost", "gaps", "avg", "sdv", if $jitter_stat;
-		    printf " %6s %4s %4s %4s %4s %4s ", "numb", "avg", "sdv", "min", "max", "rfc" if $jitter_stat;
+		    printf " %6s %4s %4s %4s %4s %4s %5s ", "numb", "avg", "sdv", "min", "max", "rfc", "MOS" if $jitter_stat;


 		    printf " %5s ", "ency", "" if $opt_mpeg;
@@ -878,9 +887,15 @@ sub display_stats{
 		if (! $jitter_stat) { $minjitter{$f}=0; $maxjitter{$f}=0;}
 		if ($ngap{$f} < 1 ){
 		    $mingap{$f}=0; $maxgap{$f}=0;
-		    printf "Resetting maxbps\n";
+		    printf "Resetting maxbps\n" if $opt_debug;
 		    $max100ms{$f}=0; $maxbps{$f}=0;
 		}
+		my $mos;
+		if ($jitter_stat){
+		    $mos=&mos_r($codec, $rtt,  avg( $sumgap{$f}, $ngap{$f}) / 1000,
+					avg($sumjitter{$f}, $njitter{$f})/1000, 
+					100*$nloss{$f}/$npkt{$f} );
+		}

 		printf "%8s %6s %6s %8s", 
 		$time, &lesbars($span,3), &lesbars($setuptime{$f},3), &lesbar($npkt{$f},8);
@@ -920,18 +935,18 @@ sub display_stats{
 		lesbar( (10**3) * $pts_stats{$f}{'max_dwell'}),
 		if $opt_mpeg;

-		printf " %5s %5s %5s %5s %5s %5s %5s",
+		printf " %5s %5s %5s %5s %5s %5s %5s ",
 		&lesbar($ntimeerr{$f}), &lesbar($dup{$f}),  &lesbar($late{$f}),  &lesbar($lost{$f}), 
 		&lesbar($nloss{$f}),  &lesbar( avg( $sumloss{$f}, $nloss{$f}) ),  
-		&lesbar( sdv($lost{$f}, $sumloss{$f},  $ssumloss{$f}) )
+		&lesbar( sdv($lost{$f}, $sumloss{$f},  $ssumloss{$f}))
 		    if $jitter_stat;

-		printf " %6s %4.1f %4.1f %4.1f %4.1f %4.1f ",
+		printf " %6s %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f ",
 		&lesbar($njitter{$f}),
 		avg($sumjitter{$f}, $njitter{$f})/1000,
 		sdv($njitter{$f}, $sumjitter{$f}, $ssjitter{$f})/1000,
 		$minjitter{$f}/1000, $maxjitter{$f}/1000,
-		$est_jitter{$f}/1000 
+			 $est_jitter{$f}/1000, $mos
 		    if $jitter_stat;

 		printf " %5s %5s %5s %5s  %s\n",
@@ -1020,8 +1035,8 @@ sub display_stats{
 				  &lesbars(sdv($ngap{$f}, $sumgap{$f},$ssgap{$f})/1000),
 				  &lesbars($mingap{$f}/1000),  &lesbars($maxgap{$f}/1000) );
 		}
-		printf "%8s %8d %4d %4d %5s ",
-		$time, $npkt{$f}, $dup{$f}+$ntimeerr{$f}, $late{$f}, &lesbar($lost{$f});
+		printf "%8s %8d %4d %4d %4s",
+		$time, $npkt{$f}, $dup{$f}+$ntimeerr{$f}, $late{$f}, &lesbar($lost{$f},4);

 		if($opt_mpeg){

@@ -1213,21 +1228,24 @@ sub pcr_update_stats
    my ($f, $tspack) = @_;

 # We might have wrapped around, which creates weird results.
-# this is just a work-around...
-#	if (($pcr_stats{$f}{'last_pcr'} > $tspack->pcr_s())  and
-#		($pcr_stats{$f}{'last_pcr'} > 47000) and
-#		($tspack->pcr_s() < 10000)) {
-#		$pcr_stats{$f}{'initialized'} = 0;
-#	}
+# this is just a work-around... - reset pcr_checks if large values
+
+	if (($pcr_stats{$f}{'last_pcr'} > $tspack->pcr_s())  and
+		($pcr_stats{$f}{'last_pcr'} > 47000) and
+		($tspack->pcr_s() < 10000)) {
+	    $pcr_stats{$f}{'initialized'} = 0;
+	}
 #
    if ($pcr_stats{$f}{'initialized'}) { 
 	$d_pcr = (($tspack->pcr_s() - $pcr_stats{$f}{'last_pcr'}) + $PCR_MAX) % $PCR_MAX;
-#Time between our local clocks
+        # Time between our local clocks
 	$d_t = time() - $pcr_stats{$f}{'real_last_time'};

 	$pcr_stats{$f}{'real_last_time'} = time();
 #		++$pcr_stats{$f}{'num'};
-
+        if($d_pcr < 0){ # resequenced
+	    $late{$f}++;
+	}
 	my $sample = abs($d_pcr - $d_t);
 	my $k = $pcr_stats{$f}{'num'};

@@ -1340,7 +1358,12 @@ sub compare_cc_dists
 	$metric = exp($diff/1000);
    }

+    if ($metric){
 	return 10/$metric;
+    } else {
+	warn "Metric invalid $metric" if $opt_debug;
+	return -1;
+    }
 }


@@ -1444,6 +1467,281 @@ sub largest_deviation



+#########################################################################
+#
+# compute the MOS value from R in the E-model
+# Gurvinder Singh, 2010-05-07
+
+
+sub init_mos{
+     $Nc = -70;	# sum of all circuit noise powers, referred to the 0 dBr point
+     $Ps = 35;	# room noise Ps at the sender side
+     $Pr = 35;	# room noise Ps at the receiver side
+     $slr = 8;	# Sender Loudness Rating
+     $rlr = 2;	# Receiver Loudness Rating
+     $Ds = 3;	# D-Value of Telephone, Sender Side
+     $Dr = 3;	# D-Value of Telephone, Receiver Side
+     $lstr = 18;	# Listener Sidetone Rating
+     $stmr = 15;	# Sender Sidetone Masking Rating
+     $Nfor = -64;	# Noise Floor at the Receiver Side
+     $telr = 65;	# Talker Echo Loudness Rating
+     $wepl = 110;	# Weighted Echo Path Loss
+}
+
+# Subtorutine to convert the natural log value to log base 10
+    sub log10 {
+	my $n = shift;
+	return log($n)/log(10);
+    }
+    
+# Subroutine to calculate the time related variable values to be used
+# in calculating the impairements caused by the delay (Id)
+    sub calc_timevar
+    {
+	my ($l_pl, $l_jitter, $l_rtt) = @_;
+	my $Ta = 2 * $l_pl + $l_jitter + $l_rtt;
+	my $T = $Ta;
+	my $Tr = 2 * $Ta;
+	
+	return $Ta, $T, $Tr;
+    }
+
+# Subroutine to calculate the Mean Opinion Score (MOS) value from the given 
+# R value obtained from E-model (ITU-T Rec. G.107). An estimate of the result
+# of a subjective test. Higher values indicate better quality
+sub calc_mos
+{
+    my ($r) = @_;
+    my $mos;
+    
+    if ($r < 0) {
+	$mos = 1;
+    } elsif ($r > 100) {
+	$mos = 4.5;
+    } else {
+	$mos = 1 + 0.035*$r + $r * ($r - 60) * (100 - $r) * 7 * (10 ** (-6));
+	if ($mos < 1) {
+	    $mos = 1;
+	}
+    }
+
+    printf "Mos value is: %.02f \n",$mos if $opt_debug;
+    return $mos;
+}
+
+# Subroutine to calculate the Ro, which represents in principle the basic
+# signal-to-noise ratio, including noise sources such as circuit noise and
+# room noise. T
+sub calc_Ro
+{
+    my $Nfo = $Nfor + $rlr;
+    my $olr = $slr + $rlr;
+    my $Pre =35 + 10 * log10(1 + (10 ** (-0.8)));
+    my $Nor = $rlr - 121 + $Pre + 0.008 * (($Pre - 35)**2);
+    my $Nos = $Ps - $slr - $Ds - 100 + 0.004 *(($Ps - $olr - $Ds - 14)**2);
+    my $No = 10 * log10((10**($Nc/10)) + (10**($Nos/10)) +
+			(10**($Nor/10)) + (10**($Nfo/10)));
+
+    my $Ro = 15 - (1.5 * ($slr + $No));
+    print "Ro is: $Ro No is: $No\n"  if $opt_debug;
+    return $Ro, $No;
+}
+
+# Subroutine to calculate Is, where Is is a combination of all impairments
+# which occur more or less simultaneously with the voice signal.
+sub calc_Is
+{
+    my ($qdu, $Ro, $No, $T) = @_;
+
+    ######### Calculate the Iolr ##########
+    my $olr = $slr + $rlr;
+    my $xolr = $olr + 0.2*(64 + $No - $rlr);
+    my $Iolr = 20 * (((1 + ($xolr/8)**8)**(1/8)) - ($xolr/8));
+
+    ######### Calculate the Ist ##########
+    my $stmro = -10 * log10((10**(-$stmr/10))+(exp(-$T/10) *
+					       (10**(-$telr/10))));
+
+    my $Ist = (12 * ((1+((($stmro - 13)/6)**8))**(1/8))) -
+	(28 * ((1+((($stmro + 1)/19.4)**35))**(1/35))) -
+	(13 * ((1+((($stmro - 3)/33)**13))**(1/13))) + 29;
+
+    ######### Calculate the Iq ##########
+    my $Q = 37;
+    $Q -= 15 * log10($qdu) if $qdu > 0;
+    my $G = 1.07 + (0.258 * $Q) + (0.0602 * ($Q**2));
+    my $Z = 1.53 - $G/40;	# 46/30 = 1.53
+    my $Y = (($Ro - 100)/15) + 5.47 - $G/9;	# 46/8.4 = 5.47
+    my $Iq = 15 * log10(1 + 10**$Y + 10 **$Z);
+
+    ######### Calculate the Is ##########
+    my $Is = $Iolr + $Ist + $Iq;
+    print "Is is: $Is, Iolr is: $Iolr, Ist is: $Ist, Iq is : $Iq\n" if $opt_debug;
+
+    return $Is, $Ist;
+}
+
+# Subroutine to calculate the Id, which represents the impairments caused
+# by the delay
+sub calc_Id
+{
+    my ($Ro, $No, $Ist, $Ta, $T, $Tr) = @_;
+
+    ######### Calculate the Idte ##########
+    my $Idte;
+    my $terv;
+    if ($T < 1) {
+	$Idte = 0;
+    } else {
+	$terv = $telr - 40 * log10((1 + ($T/10))/(1 + ($T/150))) +
+	    6 * exp(-0.3 * ($T**2));
+	if ($stmr < 9) {
+	    $terv += $Ist/2;
+	}
+	
+	my $Re = 80 + 2.5 * ($terv - 14);
+	my $Roe = -1.5 * ($No - $rlr);
+	$Idte = (($Roe - $Re)/2 +
+		 sqrt(((($Ro - $Re)**2)/4) + 100) - 1)*
+		 (1 - exp(-$T));
+    }
+
+    if ($stmr > 20) {
+	my $Idtes = sqrt($Idte**2 + $Ist**2);
+	$Idte = $Idtes;
+    }
+
+    ######### Calculate the Idle ##########
+    my $Rle = 10.5 * ($wepl + 7) * (($Tr + 1)**(-1/4));
+    my $Idle = ($Ro - $Rle)/2 + sqrt(((($Ro - $Rle)**2)/4) + 169);
+
+    ######### Calculate the Idd ##########
+    my $Idd;
+    if ($Ta > 100) {
+	my $X = log10($Ta/100) / log10(2);
+	$Idd = 25 * (((1+ ($X**6))**(1/6)) - (3 * ((1+ (($X/3)**6))**(1/6))) + 2);
+    } else {
+	$Idd = 0;
+    }
+    
+    ######### Calculate the Id ##########
+    my $Id = $Idte + $Idle +$Idd;
+    print"Id is: $Id Idte is: $Idte Idle is: $Idle Idd is: $Idd\n"  if $opt_debug;
+
+    return $Id;
+}
+
+# Subroutine to calculate the impairement value depending
+# upon the codec used in the communication
+# Ref: Short- and Long-Term Packet Loss Behavior: Towards 
+# Speech Quality Prediction for Arbitrary Loss Distributions
+# by Alexander Raake
+sub calc_codecvar
+{
+    my $l_codec = shift;
+    
+    # Initialize the impairment values
+    my $Ie = 0;	# A scalar number allocated to a specific type of
+    # impairment (codec), indicating the anticipated
+    # incremental value of impairment (decrease of the
+    # transmission rating factor R) resulting from the
+    # type of impairment (codec). 
+
+    my $Bpl = 0;	# Bpl is the robustness of a codec to random packet
+    # loss. 
+
+    my $Pc = 0.3;	# conditional packet loss probability in case of
+    # consecutive packets loss
+
+    if  ($l_codec eq "G.711") {
+	$Ie = 0;
+	$Bpl = 25.1;
+    } elsif ($l_codec eq "G.723.1") {
+	$Ie = 15;
+	$Bpl = 16.1;
+    } elsif ($l_codec eq "G.729A") {
+	$Ie = 11;
+	$Bpl = 19;
+    } elsif ($l_codec eq "iLBC") {
+	$Ie = 11;
+	$Bpl = 32;
+    }
+
+    return $Ie, $Bpl, $Pc;
+}
+
+# Subroutine to calculate the Ie-eff, which represents impairments caused
+# by low bit rate codecs and  impairment due to packet-losses of random
+# distribution
+sub calc_Ieff
+{
+    my ($Ie, $Bpl, $Pc, $Ppl) = @_;
+    my $pl = $Ppl / 100;
+    my $burstR = (1 - $pl) / (1 - $Pc);
+    my $Ieff = $Ie + (95 - $Ie)*($Ppl/(($Ppl/$burstR) + $Bpl));
+    
+    print"Ieff is: $Ieff\n" if $opt_debug;
+    return $Ieff;
+}
+
+# Subroutine to calculate the R value using E-model (ITU-T Rec. G.107)
+# R is the result of the E-model, estimating users' satisfaction for
+# transmission planning purposes.
+sub calc_r
+{
+    my ($Ro, $Is, $Id, $Ieff, $A) = @_;
+    my $r = $Ro - $Is - $Id - $Ieff + $A;
+    printf"R is: %.02f \n",$r if $opt_debug;
+    return $r;
+}
+
+# Analogue/Digital and Digital/Analogue conversions involve a distortion of
+# the signal, perceived as a quantization noise. qdu is defined as the
+# quantization noise resulting from a complete A/D encoding and D/A decoding
+# following the 8 bit A-law (or mu-law), according to ITU-T Recommendation
+# G.711. Note, that the qdu concept is no longer used for coders other than
+# G.711. If an impairment factor Ie is used for a piece of equipment, then the
+# qdu value for that same piece of equipment must not be used.
+# my $qdu = 1;		# Number of Quantization Distortion Units
+
+
+
+# Advantage factor, depends upon the communication device used
+# This factor enables the planner to take into account the fact that users
+# may accept some decrease in quality for access advantage, e.g. mobility
+# or connections into hard-to-reach regions
+# my $A = 0;	# For fixed wired phones it is 0 (Ref. ITU-T G.107)
+
+sub mos_r{
+# Time values related to the monitored connection
+#    my $codec = $opt_codec || 'G.711';
+#    my $rtt = $opt_rtt || 10;		# ms round trip time of the connection
+#    my $pl;		# RTP packet length in ms
+#    my $jitter;	        # Jitter value for this connection
+#    my $Ppl;		# Percentage of the packet loss
+    my ($codec, $rtt, $pl, $jitter, $Ppl)=@_;
+
+    &init_mos;
+
+    my $qdu = 0;		# Number of Quantization Distortion Units
+    my $A = 0;	# For fixed wired phones it is 0 (Ref. ITU-T G.107)
+
+    my @time_ret = calc_timevar($pl, $jitter, $rtt);
+    my @codec_ret = calc_codecvar($codec);
+    my @ro_ret = calc_Ro;
+    my @Is_ret = calc_Is($qdu, $ro_ret[0], $ro_ret[1],$time_ret[1]);
+    my $Id_ret = calc_Id($ro_ret[0], $ro_ret[1], $Is_ret[1], $time_ret[0],
+			 $time_ret[1], $time_ret[2]);
+    my $Ieff_ret = calc_Ieff($codec_ret[0], $codec_ret[1], $codec_ret[2], $Ppl);
+    my $r_ret = calc_r($ro_ret[0], $Is_ret[0], $Id_ret, $Ieff_ret, $A);
+    my $mos_ret = calc_mos($r_ret);
+
+    return($mos_ret);
+} # of mos_r
+
+
+################################################################################
+
 package mpegtspacket;

 ################