#!/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);
use IO::Socket::INET;
use IO::Socket::Multicast;
use IO::Select;
use Net::Pcap;
use NetPacket::Ethernet qw(:strip);
use NetPacket::IP;
use NetPacket::IP qw(:protos);
use NetPacket::UDP;
use Net::RTP;
use Net::RTP::Packet;
use Time::HiRes qw/ sleep time tv_interval gettimeofday/;
use Data::Dumper;
#use bigint qw/ hex oct/;
use Math::BigInt;
use Math::BigFloat;
use Carp;

use constant MTU => 1580;


#use mpegpackets;


$DEFAULT_PORT = 1955;
my $timeout=1 ; # frac seconds to wait for receiving packets

require "newgetopt.pl";
$usage="$0 [option]... [file...|ip/|:port]...
-list list flows in files
-dump dump data part to file
-format [full|pretty] print full numbers foror short pretty numbers
-net open network stream rather than file(s)
-crude the log is from crude
-bins bins for gaps in ms : b1,b2,,,
-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
-last seconds - exit after seconds
-packets packets - exit after packets
-src source - limit stats to given source
-flow_key flow key - limit stats to this flow(src_ip:src_port->dest_ip:dest_port)
-flow_no flow no - limit stats to this flow no from -list
-flow_min numb  - discard flows with less packets that numb
-nohead - no headings for batch use
-id name - report name instead of flow_key in flows
-sum name - make a summary of all flows to name rather than per flow
-ttl - print ttl values for stream
-verbose - more details packet distances and sizes
-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', '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){
    @bins=split(/,/, $opt_bins);
}

my $format = $opt_format || "pretty";


my %hdr=();
my $err="";
my $us=0; # current packet time in us
my $nkill=0;
my $nint=0;
my $n_packets=0;
my $endstream=0;

#MPEG-TS measurement globals (most are reset in init)
my $CHECK_NUM_TSP = 2;
#globals for mpeg ts 
my @cc = ();
my $discontinuities = 0;
my %last_ccs = ();
my $packets_between_pcr = 0;
#my $PTS_MAX = 95443.7176888888888888
my $PTS_MAX = 95443;
my $PCR_MAX = 95443;
# pcr_stats {filename/handle}
my %pcr_stats = ();
#PTS-stats {filename/handle}
my %pts_stats = ();
my %mpeg_stats = ();
#PTS stuff
#my $real_time = time();
my $late_video_packets = 0;

#CC distribution
my $num_cc = 0;
my %cc_dist = ();
my %pid_next_cc = ();

my %dist1 = ();
my %dist2 = ();

$opt_period = $opt_last if not $opt_period;
#MPEG-TS END




$SIG{USR2} = sub {confess "Caught by SIGUSR2"; };


$SIG{INT} = sub { $uninterrupted=0; # return if $nint++  < 1;
		  &display_stats();  die "End after interrupt.\n";exit(0) };
$SIG{KILL} = sub { $uninterrupted=0; # return if $nkill++ < 1;
		   &display_stats();  die "End after kill.\n";exit(0) };
$SIG{ALRM} = sub { 
#    die "No packets received" if $n_packets < 1;
#    $uninterrupted=0;
    &display_stats if ! $endstream;
#	die "End after alarm\n"; 
#	return(0);
    $endstream=1;
    print "endstream=1" if $opt_debug;
};

if ($opt_dump){
    open DUMP, ">$opt_dump";
}

# $wanted_flow = $opt_flow_key;

if ($opt_pcap){ # listen in parallell
    alarm($opt_last) if $opt_last; 
    $uninterrupted=1;
    $endstream=0;
    eat_pcap_streams(@streams); 
    &display_stats();
} else { # listen serially
    foreach $id (@streams) {
	alarm($opt_last) if $opt_last; 
	$uninterrupted=1;
	$endstream=0;

	init_file($id);
	eat_stream($id);

	if (! ( $opt_list || $opt_sum)){
	    &display_stats() if !$endstream;
	}
	last if ! $uninterrupted;
    }
}
&list_flows if $opt_list;
&display_stats() if $opt_sum;
&display_bins if $opt_bins;

exit(0);


#------------------------------------------------------------------------------

# unpack uri's to ip-address

sub prepare_ip {
    my $uri=shift;
    $uri =~ s/^(udp|rtp)://;
    $uri =~ s#//##;
    die "Unsupported streaming $1" if $uri =~/^(\w:)/;
    if ($multicast = ($uri =~ /@/)){
	($source, $uri)=split(/@/, $uri);
    }

    my ($host, $port)=split(/[\/:]/, $uri);
    die $usage . "undefined address inn $host" unless (defined $host);
    $port = $DEFAULT_PORT unless (defined $port);
    my $address=$host;
    if ($host !~ /^\d+\.\d+\.\d+\.\d+$/){ # hostname
	my $padr=gethostbyname($host) || die "Gethostbyname - $host - $padr : $!";
	$address=Socket::inet_ntoa($padr);
    }   
    my @nib=split(/\./,$address);
    $multicast = $multicast || $nib[0] >= 224 && $nib[0] <= 240;

    return ($address, $port, $multicast);

}


sub eat_stream {
    my $id=shift;
    $tstart{$id}=[gettimeofday];

    if ($opt_net && $opt_rtp){
	# Create RTP socket
	my $elapsed=0;
	my $t0=0.0;
	$n_packets=0;
	my ($address, $port, $multicast) = prepare_ip($id);

	my $rtp = new Net::RTP( LocalPort=>$port,
				LocalAddr=>$address ) || die "Failed to create RTP socket: $!";
	my $select=IO::Select->new($rtp->fileno);

	my $tjoined=[gettimeofday];
	if ($opt_align){ # align to $opt_period
	}

	while ($uninterrupted && !($opt_packets && $n_packets > $opt_packets)) {

	    my $packet;
	    last if ! $uninterrupted || $endstream;
	    if ($select->can_read($timeout) && ($packet = $rtp->recv())){
		my $ssrc = $packet->ssrc();
		$cur_id=$packet->source_ip . ", $address:$port";

		if ( !$file{$cur_id}){
		    push(@myflows, $cur_id);
		    $file{$cur_id}=$id;
		}
		$n_packets++;
		next if $opt_src and $opt_src ne $cur_id;

		$tc=[gettimeofday]; # current flow time
		if( $t0){
		    $elapsed=tv_interval ( $t0, $tc);
		} else {
		    $t0 = $tc;
		    $elapsed=0.0;
		    &init_stats($cur_id);
		}
		$us=$elapsed*10**6;
		
		$tjoined{$cur_id} = $tjoined if !$tjoined{$cur_id};
		# No stats for that SSRC yet?
		&pkt_stats($cur_id, $packet->encode, $elapsed * 10**6, $packet->size());
		if ($opt_dump){
		    print DUMP $packet;
		}
	    } else { # anticipate interrupted by alarm	
#		warn "Failed to receive $n_packets packets from $address:$port : $!";
#		next; #### exit
	    }
	}
#	&display_stats();

#	$mc->shutdown; # gives error message


    } elsif ($opt_net) {			# udp stream

	# Create UDP socket
	my $elapsed=0;
	my $t0=0.0;
        $n_packets=0;
	my ($address, $port, $multicast) = prepare_ip($id);
	my $mc;

	if ( $multicast){
	    $mc = IO::Socket::Multicast->new(LocalPort=>$port, LocalAddr=>$address, ReuseAddr=>1) ||
		die "Could not connect to port $port : $!";
	    $mc->mcast_add($address) || die "could not join $address:$port : $!";
	} else { # unicast
	    $mc = IO::Socket::INET->new( LocalAddr=> $address, LocalPort => $port);
	}
	my $tjoined=[gettimeofday];

	my $select=IO::Select->new($mc->fileno);

	if ($opt_align){ # align to $opt_period
	}

	while ($uninterrupted && ! $endstream && !($opt_packets && $n_packets > $opt_packets)){
	    
	    if ($select->can_read($timeout) && $mc->recv($packet, 1560)) {
		last if ! $uninterrupted || $endstream;
		$n_packets++;
		my $ssrc = $mc->peerhost;
		my $cur_id = $mc->peerhost . ":" . $mc->peerport  . "->$address:$port";
		if ( !$file{$cur_id}){
		    push(@myflows, $cur_id);
		    $file{$cur_id}=$id;
		}
		
		next if $opt_src and $opt_src ne $cur_id;

		$tc=[gettimeofday]; # current flow time
		if( $t0){
		    $elapsed=tv_interval ( $t0, $tc);
		} else {
		    $t0 = $tc;
		    $elapsed=0.0;
		    &init_stats($cur_id);
		}
		$us=$elapsed*10**6;
		
		$tjoined{$cur_id} = $tjoined if !$tjoined{$cur_id};
		# No stats for that SSRC yet?
		&pkt_stats($cur_id, $packet, $elapsed * 10**6, length($packet));
		if ($opt_dump){
		    print DUMP $packet;
		}
		$endstream = tv_interval($t0, $tc) >= $opt_last if $opt_last;
	    } else { # assume that we got an interrupt
#		die "Failed to receive packet: $!"; 
	    }
	}
#	&display_stats();
#        $mc->mcast_drop($address);

    } else {			# files
	my $f=$id;
	my $pcap; 
	my $tmp_file;

	my $pcap_file=$f;
	if ($f =~ /\.gz$|\.Z/){ # gzipped ?
	    $tmp_file="/tmp/pcap$$";
	    $pcap_file=$tmp_file;
	    system( "/bin/zcat $f > $tmp_file") ; # rc code ? || die "Could not unpack $f to $tmp_file : $!";
	}
	if ($opt_crude){
	    my $flow_key;
	    open CRUDE, "<$pcap_file" || die "Could not open $f : $err";
	    while(<CRUDE>){  # crude.sourceforge.net
		my $tx0;
		
		if (($fid, $seq, $src, $dst, $tx, $rx, $size)=
		    /ID=(\d+)\s+SEQ=(\d+)\s+SRC=([\d.:]+)\s+DST=([\d.:]+)\s+Tx=([\d.,]+)\s+Rx=([\d.,]+)\s+SIZE=(\d+)/){

		    $flow_key="$src->$dst";
		    $flow_pkts{$flow_key}++;
		    $tc=[int($rx), ($tx-int($rx))*10^6] ; # tval
		    if ($flow_pkts{$flow_key} <= 1) { # new flow
			push(@flows, $flow_key);  # 
			$flow_no++; 
			$file{$flow_key}=$f;
			$source{$flow_key}=$src;
			$rx0{$flow_key}=$rx;
			$tx0{$flow_key}=$tx;
			$tjoined{$flow_key}=$tc;
		    } 
		    
		    next if $opt_list;  # just count packets per flow
		    next if $opt_flow_key && $flow_key !~ $wanted_flow;
		    next if $opt_flow_no && $opt_flow_no !~ /\b$flow_no\b/;
		    
		    if ($flow_pkts{$flow_key} <= 1) { #
			&init_stats($flow_key);
			push (@myflows, $flow_key);
		    }

		    $seq_num{$flow_key}=$seq;
		    $ptx{$flow_key}=$tx{$flow_key}; # previous
		    $tx{$flow_key}=$tx;
		    $rx{$flow_key}=$rx;
		    $mtime{$flow_key}=($tx-$tx0{$flow_key})*10**6;
		    $us=($rx-$rx0{$flow_key})*10**6; # us
		    &pkt_stats($flow_key, 0, $us, $size);
		}
	    }
	    
	} else {
	    if ($pcap=Net::Pcap::open_offline($pcap_file, \$err)){ 
		#	my $select=IO::Select->new(Net::Pcap::fileno($pcap));
		
		until ( ! $uninterrupted || $endstream || ($opt_packets && $n_packets > $opt_packets) 
			|| ! (my $pkt=Net::Pcap::next($pcap, \%hdr) ) ) {
		    eat_pcap($f, $pcap, $pkt);
		    $n_packets++;
		}
		printf "%s\n", $ttl_log if  $opt_ttl;

	    } else { # open failed 
		if ($format ne "full"){ # always put out data when format full
		    die "Could not open $f : $err";
		}
	    }
	}
	undef $t0;
	unlink $tmpfile if $tmp_file;
    } # of files
} # of eat_stream


sub eat_pcap_streams{ # open groups and listen on pcap
    my @streams=@_; 
    my $filter="( ";
    my $mc = IO::Socket::Multicast->new(LocalPort=>$DEFAULT_PORT, ReuseAddr=>1) ||
	die "Could not connect to port $port : $!";

    foreach $id (@streams){ # listen to all groups
	my ($address, $port, $multicast) = prepare_ip($id);
	$filter .= " or " if $filter =~ /\(\s+\(/;
	if ($multicast){
	    $mc->mcast_add($address) || die "could not join $address:$port : $!";
	    $filter .= "(dst $address and port $port)"; 	
	} else {
	    $filter .= "(src $address and port $port)"; 
	}	
    }
    $tjoined{$f}=[gettimeofday];
    $filter .= " )";
#    $filter .= " and ip multicast";
#    $filter ="ip";

    my ($net, $netmask, $filter_t, $err, $dev);
#    printf $filter;
    $dev = Net::Pcap::lookupdev(\$err) || die "Could not lookupdev $dev : $err : $!";
    my $r = Net::Pcap::lookupnet( $dev, \$net, \$netmask, \$err);
    die "Failed lookupnet for $dev : $r : $err : $!" if $r != 0;

    $pcap=Net::Pcap::open_live( $dev, 100, 1, 0,\$err) || die "Could not open $dev : $err";
    $r=Net::Pcap::compile($pcap, \$filter_t, $filter, 0, $netmask);
    die "Could not compile filter : $r : $filter : $!)" if $r != 0;
    $r=Net::Pcap::setfilter($pcap, $filter_t);
    die "Could not set pcap filter : $r : $filter : $! " if $r != 0 ;
    my $select=IO::Select->new( Net::Pcap::fileno($pcap));

    while ($select->can_read($timeout) && (my $pkt=Net::Pcap::next($pcap, \%hdr)) ) { # get all packets
	last if ! $uninterrupted || $endstream;
	eat_pcap("pcap", $pcap, $pkt);
	$n_packets++;
	last if  $opt_packets && $n_packets > $opt_packets;
    }
    printf "%s\n", $ttl_log if  $opt_ttl;


}


sub eat_pcap { # process pcap packets
    my ($f, $pcap, $pkt)=@_;
    my $ip=null;
    my $rtp = new Net::RTP::Packet(); # spare rtp obj
    if (Net::Pcap::datalink($pcap) == 1){ # ethernet
	$ip = NetPacket::IP->decode(eth_strip($pkt));
    } elsif (Net::Pcap::datalink($pcap) == 113) { # DLT_LINUX_SLL
	my ($head, $data)=unpack('a16a*', $pkt);  
	$ip = NetPacket::IP->decode($data);
    } else {
	die "Invalid link layer type : ".Net::Pcap::datalink($pcap);
    }
    next if $ip->{proto} != IP_PROTO_UDP;
    my $udp = NetPacket::UDP->decode($ip->{data});
    my $dlen = $ip->{len} - $ip->{hlen}*4;
    my $flow_key;

    if ($opt_sum){
	$flow_key=$opt_sum;
    } else {
	$flow_key=sprintf "%s:%d->%s:%d", 
        $ip->{src_ip}, $udp->{src_port}, $ip->{dest_ip}, $udp->{dest_port};
    }
    $flow_pkts{$flow_key}++;
    if ($flow_pkts{$flow_key} <= 1) { # new flow
	push(@flows, $flow_key);  # 
	$flow_no++; 
	$file{$flow_key}=$f;
	$source{$flow_key} =sprintf "%s:%d", $ip->{src_ip}, $udp->{src_port};
    }
    if ($opt_ttl){
	if ($pttl{$flow_key} && ($opt_log ||  ($pttl{$flow_key} !=  $ip->{ttl}))){ # change in ttl
	    printf "%s\n",  $ttl_log;
	    $nttl{$flow_key}=0;		
	}
	$pttl{$flow_key}= $ip->{ttl};
	$nttl{$flow_key}++;
	($sequence_number{$flow_key}, $tx_time_seconds{$flow_key}, $tx_time_useconds{$flow_key}, $flow_id{$flow_key})= unpack('NNNN', $udp->{data});
	$ttl_log=sprintf "ID=%u SEQ=%u SRC=%s DST=%s Tx=%u.%06u  Rx=%u.%06u TTL=%u",  
	$flow_id{$flow_key}, $sequence_number{$flow_key}, 
	$ip->{src_ip}, $ip->{dest_ip}, $tx_time_seconds{$flow_key}, 
	$tx_time_useconds{$flow_key}, $hdr{tv_sec}, $hdr{tv_usec}, $pttl{$flow_key};
	$ttl_log.= sprintf(" No=%08d\n",  $nttl{$flow_key},) if !$opt_log;

    }
    next if $opt_list;  # just count packets per flow
    next if $opt_flow_key && $flow_key !~ $wanted_flow;
    next if $opt_flow_no && $opt_flow_no !~ /\b$flow_no\b/;

    if ($flow_pkts{$flow_key} <= 1) { #
	&init_stats($flow_key);
	push (@myflows, $flow_key);
#    printf "$flow_key\n";
    }
    $tc= [$hdr{tv_sec}, $hdr{tv_usec}];
    if (! defined($t0)) {
	$t0=$tc;
    }
    $tjoined{$f}=$tc if !$tjoined{$f};
    $tjoined{$flow_key} = $tjoined{$f} if !exists($tjoined{$flow_key});

    $us = tv_interval ( $t0, $tc) * 10**6 ; # us relative packet time 	  

# fixate on first flow to appear
#    if ( $first_port) {
#	next if $udp->{dest_port} != $first_port;
#    } else {
#	$first_port =  $udp->{dest_port};
#    }

    &pkt_stats($flow_key, $udp->{data}, $us, $dlen);
    
}

sub init_stats {
    my $f=shift;

    $pinterval=0;
    $npkt{$f}=$sumbyte{$f}=$ssbyte{$f}=0;
    $first_sec=0;
    delete $pmtime{$f}; delete $pus{$f} ; delete $first_us{$f};  delete $last_us{$f}; 
    delete $prev_1s{$f}; delete $prev_100ms{$f};
    delete $per1s{$f}; delete $per100ms{$f};
    delete $sum1s{$f}; delete $sum100ms{$f};
    $n100ms{$f}=0; $sum100ms{$f}=0; $ss100bps{$f}; $akk100ms{$f}=0;
    $sumbyte{$f}=$sumps{$f}=0;			# sum byte per second
    $maxbps{$f}=0; $max100ms{$f}=0;
    $dup{$f}=$lost{$f}=$late{$f}=0;
    $maxbyte{$f}=0; $minbyte{$f}=2**31;
    $sumgap{$f}=$ssgap{$f}=0;
    $ngap{$f}=0; $maxgap{$f}=0, $mingap{$f}=2**63-5;
    $njitter{$f}=$sumjitter{$f}=$ssjitter{$f}=0;
    $maxjitter{$f}=-2**63-5; $minjitter{$f}=2**63-5; delete $est_jitter{$f}; delete $pest_jitter{$f}; 
    $first_port = 0;  $n_seq=-1;
    delete $tstart{$f}; delete $tend{$f}; 

#MPEG-TS initalizing and resetting (new period!)
    $discontinuities = 0;
    %last_ccs = ();
    @cc = ();
# r_stats {filename/handle}
    @pcr_stats{$f} = ();
#P-stats {filename/handle}
    @pts_stats{$f} = ();
    @mpeg_stats{$f} = ();
#P stuff
    $late_video_packets = 0;
#MPEG-TS END
    %dist1 = ();
    %dist2 = ();
    $num_cc = 0;

}

sub init_file {
    undef $t0;
    my $f=shift;
    undef $tjoined; delete $tjoined{$f}; delete $setuptime{$f}; delete($joined{$file{$f}});
    delete $file{$f};
    @myflows=() if !$opt_sum;
    $flow_pkts{$f}=0;  # new flow on new file

    undef $pus;
    undef $psec;
    delete($pmtime{$f});
    delete($pseq_num{$f});
    delete($pus{$f});
}

sub seq_stat {
    my ($f, $seq_num, $mtime)=@_;

    my $seq_diff= $seq_num - $pseq_num{$f};
    unless ($opt_crude && $seq_diff == 1){ #
	$gap = $us - $pus{$f};		# packet time gap
	$ngap{$f}++;
	$sumgap{$f} +=  $gap;
	$ssgap{$f} += $gap**2;
	$maxgap{$f}=  $gap if $gap > $maxgap{$f};
	$mingap{$f}=  $gap if $gap < $mingap{$f};

	foreach $i (0 .. $#bins){
	    if( $gap < ($bins[$i]*1000)){
		$nbin[$i]++;
		last;
	    }
	}
	printf "Tx=%.6f gap=%d gap=%.1fms\n", $ptx{$f}, $seq_diff-1, $gap/1000, if $opt_log;
    }

    if ($seq_diff != 1) {
#	printf "seq %d,", $rtp->seq_num;

	if ($seq_diff == 0) { # Duplicated
	    $dup{$f} ++;
	} elsif ($seq_diff < 0) { # Out Of Order
	    $late{$f}++;
#		$lost{$f}--;
	} else {			# Lost                   
	    $lost{$f}+= $seq_diff - 1;
	    $nloss{$f}++;
	    $sumloss{$f} += $seq_diff;
	    $ssumloss{$f} += $seq_diff**2;
	    $maxloss{$f} = $seq_diff if  $seq_diff >  $maxloss{$f};
	}
    } elsif ( $mtime <= $pmtime{$f}){
	$ntimeerr{$f}++;
    } else { # the jitter can be measured
#  compute jitter relative to previous packet
	$mgap= ($mtime - $pmtime{$f} ); # media time gap in samples 8khz=us
	$jitter=abs($gap - $mgap);
#ok compute jitter relative to start
#ok	    $mgap= ($mtime - $first_mtime); # media time gap in us
#ok	    $jitter=abs($us - $mtime - $first_jitter);
	$njitter{$f}++;
	$sumjitter{$f} +=  $jitter;
	$ssjitter{$f} += $jitter**2;
	$maxjitter{$f}=  $jitter if $jitter > $maxjitter{$f};
	$minjitter{$f}=  $jitter if $jitter < $minjitter{$f};

	# rfc 3550 
	if (defined($est_jitter{$f})){
	    $est_jitter{$f} = $est_jitter{$f} + ( $jitter - $est_jitter{$f} )/16;
	} else { # inital value
	    $est_jitter{$f} = $jitter;
	}
    }

}


sub pkt_stats {
    my $f = shift; # file or stream name
    my $packet = shift;
    my $us = shift;
    my $dlen = shift;
    my $rtp = new Net::RTP::Packet();

    if (! exists($tstart{$f})){
	$tstart{$f} = $tc;
	$first_us{$f}=$prev1s_us{$f}=$prevt100ms_us{$f}=$us;
	if (! exists($setuptime{$f}) ){
	    if ( $pinterval){
		$setuptime{$f}=0;
	    } else {
		$setuptime{$f}=tv_interval ( $tjoined{$f}, $tstart{$f}) * 10**3;
	    }
	    printf "setup $f : %.3f\n",  $setuptime{$f} if $opt_debug;
	}
    }
    $tend{$f} = $tc;
    $last_us{$f}=$us;
    
    #	next if $dlen < 20;  # typically small background noice packets added

    $npkt{$f}++;
    $sumbyte{$f}+=$dlen; $sum1s{$f}+=$dlen; $sum100ms{$f}+=$dlen;
    $ssbyte{$f}+=$dlen**2;
    $maxbyte{$f}=  $dlen if $dlen > $maxbyte{$f};
    $minbyte{$f}=  $dlen if $dlen < $minbyte{$f};

    $sumps{$f}+=$dlen;		# sum per second

    if ($opt_dump){
	print DUMP $packet;
    }

    if ($opt_mpeg){ 
	mpeg_stats($f, $packet, $us);
    }

    if ($opt_rtp){
	if ( $rtp->decode($packet) ) { #assume RTP
#	    printf "RTP: %10s %5d %5d\n", $rtp->ssrc(), $rtp->timestamp(), $rtp->size();
#	    $dlen = $rtp->size();
	} else {
	    print "??? $f # $npkt{$f} $dlen B\n";
		if ($opt_debug){
		    my $n=0;
		    foreach $w( unpack("C*", $packet)) {
			printf "%0X",$w;
			if ($n > 0){
			    print " " if $n % 4 == 0;
			    print "\n" if $n % 32 == 0;
			}
			$n++;
		    }
		    printf "\n";
		}
	    next;
	}
	# Calculate next number in sequence
	$n_seq = $rtp->seq_num()+1;
	if ($n_seq > 65535) {
	    $n_seq=0;
	}
	$mtime{$f}=$rtp->timestamp() * &time_unit($rtp->payload_type);  # mediatime in us
	$seq_num{$f}=$rtp->seq_num;
    }

    if (defined($pus{$f})) {
	$gap = $us - $pus{$f};		# packet time gap
	$ngap{$f}++;
	$sumgap{$f} +=  $gap;
	$ssgap{$f} += $gap**2;
	$maxgap{$f}=  $gap if $gap > $maxgap{$f};
	$mingap{$f}=  $gap if $gap < $mingap{$f};

	if ($opt_rtp){

	    &seq_stat($f, $seq_num{$f}, $mtime{$f});

	}
	if ($adapt == 1 ) { # adaption field coming
	    my ($length, $aflags, $rest2)=unpack "C1C1C*", $rest;
	    my ($dif, $raif, $espif, $pcrf, $opcrf, $splicef, $tpdf, $afextf)=split(//, $aflags); 
	    if ($pcrf == 1){
		;
	    }
	    printf "%5d %6d %5d %5d %5d %5d %5d\n",$rtp->seq_num,$us,$rtp->timestamp(), $gap,  $mgap, $jitter, $first_jitter if $opt_debug;
	    
	} elsif ($opt_crude){
	    &seq_stat($f, $seq_num{$f}, $mtime{$f});

	} elsif ($opt_mpeg){ # incomplete code under exploration
	    # meg transport stream layout http://en.wikipedia.org/wiki/Transport_stream
	    my ($flag, $b2, $b3, $rest)= unpack "H2n1C1C*", $packet;
#	my ($ts_err, $payload, $ts_pri)=split( //, unpack "B3", $b2); # does not work
	    my ($ts_err, $payload, $ts_pri) = ( $b2 & 0x8000 , $b2 & 0x4000,  $b2 & 0x2000);
	    my $pid= $b2 & 0x1F00; 
#	my ($scramble1, $scramble2, $adapt, $data) = split( //, unpack "B4", $b3); # not working
	    my ($scramble1, $scramble2, $adapt, $data) = ( $b3 & 0x80,  $b3 & 0x40, $b3 & 0x20, $b3 & 0x10);
	    my $contin= $b3 & 0x0F;
#	printf "%5d %5s %5d %5X %5d %5d %s\n", $npkt{$f}, $flag, $pid, $b2, $data, $contin, unpack("H*",$packet)
#	     if $payload > 0 ; # $contin == 13;  #$ts_pri > 0 ;  # $b2 > 256;

	    if ($payload && $pid > 0){  # PES follows http://en.wikipedia.org/wiki/Packetized_Elementary_Stream
		my ($ts_prefix, $code, $id, $pes_lth, $pes_rest )= unpack("H8H6C1C2C*", $packet);
#	    printf "%5d %5s %5d %5X %5d %5d %s\n", $npkt{$f}, $flag, $pid, $code, $data, $id, unpack("H*",$packet);

		if ($code == 0x000001){
		    if ( $id >= 0xE0 ){ # video PES-header
			my ($ts_prefix, $pes_head, $pes_opt, $head_lth)=unpack("H8H6C2C1", $packet);
			if ($pes_opt & 0x0080){ # PTS ind
			    printf "%5d %5s %5d %5X %5d %5d %s\n", $npkt{$f}, $flag, $pid, $id, $pes_lth, $head_lth, unpack("H*",$packet) if $opt_debug;
			}
		    }
		    
		}
	    }
	    if ($adapt == 1 ) { # adaption field coming
		my ($length, $aflags, $rest2)=unpack "C1C1C*", $rest;
		my ($dif, $raif, $espif, $pcrf, $opcrf, $splicef, $tpdf, $afextf)=split(//, $aflags); 
		if ($pcrf == 1){
		    ;
		}
	    }
	}

	&burst($f);

    }
    if ($opt_period){
	$interval=int($us/10**6 / $opt_period);
	if ($interval > $pinterval){
	    &mpeg_debug() if $opt_debug;
	    &display_stats() ;
	    $pinterval=$interval;
	    $pinterval{$f}=$interval;
	    $endstream = $opt_nperiod && $interval >= $opt_nperiod;
	}
    }

    $pus=$us;
    $psec=$sec;
    $pmtime{$f}=$mtime{$f} if $opt_rtp || $opt_crude;
    $pseq_num{$f}= $seq_num{$f} if $opt_rtp ||  $opt_crude;
    $pus{$f}=$us;
} # of pkt_stats


sub burst {
    my $f=shift;
    my $span=$last_us{$f}-$first_us{$f};
    my $per1s=int($span/10**6);
    #printf "per1s: %f\n", $per1s;
    my $per100ms=int($span/10**5);
    #printf "per1s_f%f\n", $per1s{$f};
    
    if ( $per1s > $per1s{$f} ){ # period expired
	my $per=($last_us{$f}-$prev_1s{$f})/10**6;
	#		printf "per: %f\n", $per;
	if ($per>0){
	    my $bps = $sum1s{$f} * 8 / $per;
	    $maxbps{$f} =  $bps if $bps > $maxbps{$f};
	    $per1s{$f} =  $per1s;
	    $prev_1s{$f} = $last_us{$f};
	    $sum1s{$f}=0;
	}
    }
    if ( $per100ms > $per100ms{$f} ){ # period expired
	my $per=($last_us{$f}-$prev_100ms{$f})/10**6;
	if ($per>0){
	    my $bps = $sum100ms{$f} * 8 / $per;
	    $max100ms{$f} =  $bps if $bps > $max100ms{$f};
	    $n100ms{$f}++;
	    $ss100bps{$f}+= $bps**2;
	    $akk100ms{$f}+=$bps;
	    $per100ms{$f} =  $per100ms;
	    $prev_100ms{$f} = $last_us{$f};
	    $sum100ms{$f}=0;
#	print "$bps ";
	}
    }

}

sub display_stats{

#    my $f = shift;

    my $media="file";
#    $media = "host, group" if $opt_net;
    $media = "host:port" if $opt_net;

    @myflows=@streams if ($#myflows < 0); # no response

    foreach $f (sort @myflows) { 
	my $id=$f;
	next if $npkt{$f} <= $opt_flow_min;
	my $source=$f; 
	$source=$source{$f} if  ! $opt_v && $source{$f} ne "";
	$source=$opt_id if $opt_id;
	my $jitter_stat= $opt_rtp || $opt_crude;

	my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime($tstart{$f}->[0]);
	my $time=sprintf "%02d:%02d:%02d", $hour, $min, $sec;
	if ($format eq "full"){
	    $time.=sprintf ".%06d", $tstart{$f}->[1];
	} else {
	    my $time=sprintf "%02d:%02d:%02d", $hour, $min, $sec;
	}
	my $date=sprintf "%4d-%02d-%02d", $year+1900, $mon+1, $mday;
	my $span=&tv_interval($tstart{$f}, $tend{$f}); # second

	&burst($f); # in case a second is interrupted
#	printf "dus %d min %.2f\n", $us-$pus, $max100ms{$f};
#	$max100ms{$f}=-1  if  $max100ms{$f} > 2**61;
#	$maxbps{$f}=-1  if ! $maxbps{$f} > 0;
	
	if ($opt_v) {

	    if ($date ne $pdate) {
		if (! $opt_nohead){
		    print "\nDate $date            packet  size(byte)           packet distance(ms) ";
		    print "       time    sequence         loss               packet jitter(ms)        " if $jitter_stat;

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

		    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 %5s ", "numb", "avg", "sdv", "min", "max", "rfc", "MOS" if $jitter_stat;


		    printf " %5s ", "ency", "" if $opt_mpeg;
		    printf " %4s %5s %5s %5s ", "avg", "sdv", "min", "max" if $opt_mpeg;
		    printf " %4s %5s %5s %5s ", "avg", "sdv", "min", "max" if $opt_mpeg;

		    printf "  %5s %5s %5s %5s  %-15s\n", 
		    "avg",  "1s", "100ms", "sdv", $media;
		}

	    }

	    if ( $npkt{$f} >= 0 && $ngap{$f}>= 0) {
		if (! $jitter_stat) { $minjitter{$f}=0; $maxjitter{$f}=0;}
		if ($ngap{$f} < 1 ){
		    $mingap{$f}=0; $maxgap{$f}=0;
		    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);
		my $n=$npkt{$f}; $n=1 if $npkt{$f} <1; # divide by zero protection
		printf " %4s %4s %4d %4d  ", &lesbar($sumbyte{$f}/$n,4), 
		&lesbar(sdv($npkt{$f}, $sumbyte{$f}, $ssbyte{$f})),	
		$minbyte{$f}, $maxbyte{$f};

		printf " %6s %4s %4s %4s %4s ", 
		&lesbar($ngap{$f},6), &lesbars( avg( $sumgap{$f}, $ngap{$f}) / 1000 ),
		&lesbars(sdv($ngap{$f}, $sumgap{$f},$ssgap{$f})/1000),
		&lesbars($mingap{$f}/1000),  &lesbars($maxgap{$f}/1000);

#				printf (" %5s %.1e",
#					&lesbar($mpeg_stats{$f}{'num'},3),
#					(&largest_deviation(1/16, $mpeg_stats{$f}{'num'}, \@cc),6)) if $opt_mpeg;
		printf (" %4s",
			&lesbar(compare_cc_dists(\%dist1,\%dist2)),
			"") if $opt_mpeg;


		printf " %5s %5s %5s %5s", 
		&lesbar( (10**3) * avg($pcr_stats{$f}{'pcr_jitter_sum'}, $pcr_stats{$f}{'num'})),
		&lesbar( (10**3) * sdv($pcr_stats{$f}{'num'}, 
				       $pcr_stats{$f}{'pcr_jitter_sum'}, 
				       $pcr_stats{$f}{'pcr_jitter_square'})),
		&lesbar( (10**3) * $pcr_stats{$f}{'min_jit'}),
		&lesbar( (10**3) * $pcr_stats{$f}{'max_jit'})
		    if $opt_mpeg;

		printf " %5s %5s %5s %5s ", 
		lesbar( (10**3) * avg($pts_stats{$f}{'dwell_sum'}, $pts_stats{$f}{'num'})),
		lesbar( (10**3) * sdv($pts_stats{$f}{'num'}, 
				      $pts_stats{$f}{'dwell_sum'}, 
				      $pts_stats{$f}{'dwell_squared'})),
		lesbar( (10**3) * $pts_stats{$f}{'min_dwell'}),
		lesbar( (10**3) * $pts_stats{$f}{'max_dwell'}),
		if $opt_mpeg;

		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}))
		    if $jitter_stat;

		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, $mos
		    if $jitter_stat;

		printf " %5s %5s %5s %5s  %s\n",
		&lesbar( avg( $sumbyte{$f}*8, $span ) ), 
		&lesbar($maxbps{$f}), 
		&lesbar($max100ms{$f}), 
		&lesbar(sdv($n100ms{$f},
			    $akk100ms{$f}, 
			    $ss100bps{$f})), 
		$source;


	    } else { # not reached
		printf "%-8s %6s %6s %8s %4s %4d %4d %4d %6d %4s %4s %4s %4s %5d %5d %5d %5d %s\n", $time, 0, 0, &lesbar($npkt{$f},4), 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, $source if $f ne "" && $npkt{$f}>=0;
	    }

	} else {  #not -v
	    if ($date ne $pdate) {
		if (! $opt_nohead){
		    if($opt_mpeg){
			print "\nDate $date       packet      flu-    ";
		    } else {
			print "\nDate $date       packet      size(byte)   ";
		    }
		    if ($jitter_stat){
			print "jitter(ms)                " ;
		    } elsif ($opt_mpeg) {
			print "pcrjit(ms)  dwell (ms)    " ;
		    } else  {
			print "gap(ms)                   " ;
		    }
		    print "thrust(bps)       source\n";                      

		    if($opt_mpeg) {
			printf "%-8s %8s %4s %4s %4s %4s ", 
                        #"time","numb", "dup", "late", "lost", "numb", "dscnt";
			"time","numb", "dup", "late", "lost", "ency";
		    } else {
			printf "%-8s %8s %4s %4s %4s %4s %4s", 
			"time","numb", "dup", "late", "lost", "avg", "sdv";
		    }	


		    if($opt_mpeg) {
			printf "   %5s %5s %5s %5s",
			"avg", "sdv", "avg", "sdv";

		    } else {
			printf "  %5s %5s %5s %5s",
			"avg", "sdv", "min", "max" ; # if $jitter_stat;
		    }	

		    printf "    %4s %5s %5s   %s\n",
		    "avg", "1s", "100ms", $media;
		}

	    }


	    if ( $npkt{$f}>= 0 && $ngap{$f}>= 0) { 
		my @jitterstat=(0,0,0,0);
		my $n=$npkt{$f}; $n=1 if $npkt{$f} <1; # divide by zero protection

		if ($jitter_stat){
		    @jitterstat= (&lesbars($sumjitter{$f}/$njitter{$f}/1000),
				  &lesbars(sdv($njitter{$f}, $sumjitter{$f}, $ssjitter{$f})/1000),
				  &lesbars($minjitter{$f}/1000), 
				  &lesbars($maxjitter{$f}/1000));
		} elsif ($opt_mpeg) {


		    @jitterstat = (

				   &lesbar( (10**3) * avg($pcr_stats{$f}{'pcr_jitter_sum'}, $pcr_stats{$f}{'num'})),
				   &lesbar( (10**3) * sdv($pcr_stats{$f}{'num'}, 
							  $pcr_stats{$f}{'pcr_jitter_sum'}, 
							  $pcr_stats{$f}{'pcr_jitter_square'})),
				   &lesbar( (10**3) * avg($pts_stats{$f}{'dwell_sum'}, $pts_stats{$f}{'num'})),
				   &lesbar( (10**3) * sdv($pts_stats{$f}{'num'}, 
							  $pts_stats{$f}{'dwell_sum'}, 
							  $pts_stats{$f}{'dwell_squared'})),
				   );	

		} else {
		    @jitterstat=( &lesbars( avg($sumgap{$f}, $ngap{$f}) / 1000),
				  &lesbars(sdv($ngap{$f}, $sumgap{$f},$ssgap{$f})/1000),
				  &lesbars($mingap{$f}/1000),  &lesbars($maxgap{$f}/1000) );
		}
		printf "%8s %8d %4d %4d %4s",
		$time, $npkt{$f}, $dup{$f}+$ntimeerr{$f}, $late{$f}, &lesbar($lost{$f},4);

		if($opt_mpeg){

		    printf (" %4s  ",  
			    &lesbar(compare_cc_dists(\%dist1, \%dist2))) if $opt_mpeg;

		} else {
		    printf "%4d %4.1f", $sumbyte{$f}/$n, sdv($n, $sumbyte{$f}, $ssbyte{$f});
		}

		printf   "  %5s %5s %5s %5s",@jitterstat;
		printf "   %5s %5s %5s  %s\n", 
		&lesbar( avg($sumbyte{$f}*8, $span)), 
		&lesbar($maxbps{$f}), 
		&lesbar($max100ms{$f}),
		$source;

		#printf "in-print: \n\n%f\n\n", $maxbps{$f};
	    } else {
		printf "%-8s %8s %68s     %s\n",$time, lesbar($npkt{$f},4),"", $source if $f ne "" && $npkt{$f}>0;
	    }

	}
	$opt_nohead=1;
	$pdate=$date;
	$flow_pkts{$f}=0;
	&init_stats($f);

    }
}

sub display_bins{
    print "\n\nBins    ";
    foreach $bin (@bins){
	printf "%5d ", $bin;
    }
    print "\nGaps    ";
    foreach $i (0 .. $#bins){
	printf "%5d ", $nbin[$i];
    }
    print "\n";
}

sub avg {  
    my ($sum, $n)=@_;
    return 0 if $n <= 0;
    return $sum/$n;
}
sub sdv {
    my ($n, $sum, $sumsq)=@_;
    if ($n > 1) {
	return sqrt( abs($sumsq/$n -  ($sum/$n)**2));
    } else {
	return  0;
    }
}


sub lesbar{
    my ($n, $figs)=@_;
    $figs = 3 if ! $figs;
    return sprintf ("%.1f", $n) if $format ne "pretty";

    @dekade=("", "K", "M", "G", "T", "P");
    for ($d=0; $d<=$#dekade; $d++){
	$dekade=10**(3*$d);
	my $nd=$n / $dekade;

	if ( $nd < 10**$figs){ # 3 siffer OK
	    if ( $nd > 100 || ($n < 10000 && $n == int($n))){
		return (sprintf "%d", $nd) . $dekade[$d];
	    } else { 
		return sprintf "%.1f%s", $nd, $dekade[$d];
	    }

	}
    }
    return sprintf "%E", $n;

}

sub lesbars {
    my ($n, $figs)=@_;
    return &lesbar($n, $figs) if $format ne "pretty";
    if ($n > 1000){
	my $num = &lesbar($n /1000, $figs) . "s";
	$num=~s/ //; 
	return $num;
    } else {
	return  &lesbar($n, $figs);
    }

}

# compute timeunits in mediatime from rfc 3551
# limited to static payload-types
sub time_unit{
    my $payload_type=shift;
    $sample_freq{0}=8000;
    $sample_freq{3}=8000;
    $sample_freq{4}=8000;
    $sample_freq{5}=8000;
    $sample_freq{6}=16000;
    $sample_freq{7}=8000;
    $sample_freq{8}=8000;
    $sample_freq{9}=8000;
    $sample_freq{10}=44100;
    $sample_freq{11}=44100;
    $sample_freq{12}=8000;
    $sample_freq{13}=8000;
    $sample_freq{14}=90000;
    $sample_freq{15}=8000;
    $sample_freq{16}=11025;
    $sample_freq{17}=22050;
    $sample_freq{18}=8000;
    $sample_freq{726}=8000;
    $sample_freq{726}=8000;
    $sample_freq{726}=8000;
    $sample_freq{726}=8000;
    $sample_freq{729}=8000;
    $sample_freq{729}=8000;
    $sample_freq{25}=90000;
    $sample_freq{26}=90000;
    $sample_freq{28}=90000;
    $sample_freq{31}=90000;
    $sample_freq{32}=90000;
    $sample_freq{33}=90000;
    $sample_freq{34}=90000;

    return 10**6/($sample_freq{$payload_type}||8000);
}

sub list_flows {
    my $n=0;
    printf "\n\n%4s   %-50s %10s    %s\n", "No", "    Flow key", "Packets", "File" if  ! $opt_nohead;;
    foreach $f ( @flows){
	printf "%4d   %-50s %10d    %s\n", $n+++1, $f, $flow_pkts{$f}, $file{$f};
    }    
    $opt_nohead=1;
}


sub mpeg_stats
{
    my ($flow_id, $packet, $us) = @_;
#Assume the packet-size is 188 bytes. Might be 204 if error-correcting codes 
#are in use.
    my $l = length ($packet);
#my $n = $l / 188;
    my $n = $CHECK_NUM_TSP;

#One UDP-packet contains $n number of MPEG-TS packets, check every one:
    for(my $j = 0; $j < $n; ++$j) {
	$mpeg_stats{$flow_id}{'num'} += 1;
	my $tspack = mpegtspacket->new (substr($packet, $j*188, 188));


#Update continuity counters
	cc_update_stats($flow_id, $tspack, $us);


	if ($tspack->is_valid() != 1) { next; }

#PCR calcs			
	if ($tspack->pcr_flag() == 1) {
	    pcr_update_stats($flow_id, $tspack);
	}


#Dwell calculations
	if ($tspack->payload_unit_start() == 1) {
	    dwell_update_stats($flow_id, $tspack);
	} else {
	    next;
	}

    }

#Set the time for the last UDP packet
    $mpeg_stats{$flow_id}{'last_udp_time'} = time();
}



#This function uses globals:
#%pcr_stats()
sub pcr_update_stats
{
    my ($f, $tspack) = @_;

# We might have wrapped around, which creates weird results.
# 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
	$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'};

	$pcr_stats{$f}{'pcr_jitter_sum'} += $sample;
	$pcr_stats{$f}{'pcr_jitter_square'} += ($sample**2);
	++$pcr_stats{$f}{'num'};
#Max/Min-checks for PCR jitter

	if ($sample > $pcr_stats{$f}{'max_jit'}) {
	    $pcr_stats{$f}{'max_jit'} = $sample;
	} 
	if ($sample < $pcr_stats{$f}{'min_jit'}) {
	    $pcr_stats{$f}{'min_jit'} = $sample;
	}

#This calculates the PCR arrival difference mean and std 	
#$sample = $d_pcr;
	$tm = $pcr_stats{$f}{'running_dpcr_mean'};
	$ts = $pcr_stats{$f}{'running_dpcr_std'};
	$pcr_stats{$f}{'running_dpcr_mean'} = $tm + (($sample - $tm)/$k);
	$tm2 = $pcr_stats{$f}{'running_dpcr_mean'};	
	$pcr_stats{$f}{'running_dpcr_std'} = $ts + (($sample - $tm)*($sample - $tm2));
#		print $pcr_stats{$f}{'running_dpcr_mean'}."\n";

#Max/Min-checks for delta-PCR
	if ($sample > $pcr_stats{$f}{'max_dpcr'}) {
	    $pcr_stats{$f}{'max_dpcr'} = $sample;
	} 
	if ($sample < $pcr_stats{$f}{'min_jit'}) {
	    $pcr_stats{$f}{'min_dpcr'} = $sample;
	}



	$pcr_stats{$f}{'last_pcr'} = $tspack->pcr_s();
	$pcr_stats{$f}{'last_pcrb'} = $tspack->pcrb_s();

	$packets_between_pcr = 0;

#Got a new PCR-value, so synch the clocks
    } else {
#This is the first time we see a PCR, intialize hash.
#Time between this pcr and the last measured.
#$d_pcr = (($tspack->pcr() - $pcr_stats{$f}{'last_pcr'}))/(27*10**6);
#Time between our local clocks
#$d_t = time() - $pcr_stats{$f}{'real_last_time'};


	$pcr_stats{$f}{'running_jit_mean'} = 0; #($d_pcr - $d_t);
	$pcr_stats{$f}{'running_jit_std'} = 0;
	$pcr_stats{$f}{'running_dpcr_mean'} = 0; #first time
	$pcr_stats{$f}{'running_dpcr_std'} = 0;
	$pcr_stats{$f}{'real_last_time'} = time();
	$pcr_stats{$f}{'last_pcr'} = $tspack->pcr_s();
	$pcr_stats{$f}{'last_pcrb'} = $tspack->pcrb_s();
	$pcr_stats{$f}{'pcr_jitter_sum'} = 0;
	$pcr_stats{$f}{'num'} = 0;
	$pcr_stats{$f}{'max_jit'} = 0;
	$pcr_stats{$f}{'min_jit'} = (2**32)-1;
	$pcr_stats{$f}{'max_dpcr'} = 0;
	$pcr_stats{$f}{'min_dpcr'} = (2**32)-1;
	$pcr_stats{$f}{'initialized'} = true;
	$pcr_stats{$f}{'missing_pcrs'} = 0;
    }

}

#uses :
#last_ccs
#
sub cc_update_stats
{
    my ($f, $tspack, $us) = @_;
    my $secs = $us*10**(-6);
    my $prior = $last_ccs{$pid};
    $last_ccs{$tspack->PID()} = $tspack->continuity_counter();
    check_continuity($tspack->PID(), 
		     $prior,
		     $tspack->{'discontinuity'}, 	
		     $tspack->payload_unit_start() );

#printf("%d %d\n", $tspack->PID(), $tspack->continuity_counter());

#update continuity overview
    $last_ccs{$tspack->PID()} = $tspack->continuity_counter();

#update the per-counter continuity.
    $cc[$tspack->continuity_counter()]++;

#distribution updates
    if(!$opt_period || ($secs % $opt_period < ($opt_period/2))) {
	$dist1{$tspack->PID()} += 1;
    } elsif($secs % $opt_period <= ($opt_period)) {
	$dist2{$tspack->PID()} += 1;
    } 
}

sub compare_cc_dists
{
    my ($a, $b) = @_;
#	print Dumper(%$a); 
#	print "\n";
#	print Dumper(%$b); 
    my $metric = 0;
    foreach my $pid (keys %$a) {
	my $diff = $a->{$pid} - $b->{$pid};
	my $sum = $a->{$pid} + $b->{$pid};
#	print "sum: $sum diff: $diff\n";
#       $metric += $diff / $sum;
	$metric = exp($diff/1000);
    }

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



sub dwell_update_stats
{
    my ($f, $tspack) = @_;
    my $pespack = mpegpespacket->new ($tspack->payload());	
    next if ($pespack->is_valid != 1);

#pts_stuff
    if($pespack->{'PTS_flag'} == 1) {
	++$pts_stats{$f}{'num'};

	if($pts_stats{$f}{'initialized'} == 1) {
#check if frame is late
	    if($pespack->pts_s() < $pcr_stats{$f}{'last_pcr'}) {
		$pts_stats{'late_frames'}++;
	    }

# 4 - Dwell-time calculations
	    if($pcr_stats{$f}{'last_pcr'} > 0 and $pespack->pts_s() > 0) {   #we want to compare against valid PCR
		my $dwell = ((abs($pcr_stats{$f}{'last_pcr'} - $pespack->pts_s())) + $PTS_MAX) % $PTS_MAX;
		$pts_stats{$f}{'dwell_sum'} += $dwell;
		$pts_stats{$f}{'dwell_squared'} += $dwell**2;
		$pts_stats{$f}{'last_dwell'} = $dwell;

#min/max calcoolatins
		if ($dwell > $pts_stats{$f}{'max_dwell'}) {
		    $pts_stats{$f}{'max_dwell'} = $dwell;
		} 
		if ($dwell < $pts_stats{$f}{'min_dwell'}) {
		    $pts_stats{$f}{'min_dwell'} = $dwell;
		}
	    }

	} else {
	    my $dwell = ((abs($pcr_stats{$f}{'last_pcr'} - $pespack->pts_s())) + $PTS_MAX) % $PTS_MAX;
	    $presentation_time = $pespack->pts_s();
	    $pts_stats{$f}{'mean_dwell'} = 0;
	    $pts_stats{$f}{'dwell_squared'} = 0;
	    $pts_stats{$f}{'std_dwell'} = $dwell;
	    $pts_stats{$f}{'max_dwell'} = 0;
	    $pts_stats{$f}{'min_dwell'} = (2**32)-1;
	    $pts_stats{$f}{'last_dwell'} = $dwell;
	    $pts_stats{$f}{'late_frames'} = 0;
	    $pts_stats{$f}{'num'} = 0;
	    $pts_stats{$f}{'initialized'} = 1;
	}
    }

}


sub check_continuity
{
    my ($pid, $cc, $discont, $payload) = @_;

    if (not exists $last_ccs{$pid}) {
	return;
    }	
    if ($last_ccs{$pid} == $cc and $payload == 0) {
	return;	
    }
    if (($last_ccs{$pid} + 1)%16 == $cc) {
	return;
    }
    if ($last_ccs{$pid} == $cc) {
	return;
    }
    $discontinuities++;
}




sub mpeg_debug
{
    print "Debug: ";
    print Dumper(%pcr_stats);
    print Dumper(%pts_stats);
    print Dumper(%mpeg_stats);
    print Dumper(@cc);
}

sub largest_deviation
{
    my ($expected, $num, $values) = @_;
#	print "\n\nhnnng: $num\n\n";
    my $max = 100000000;
    foreach $id (@$values) {
	my $val = abs(($id/$num) - $expected) if $num != 0;
#		printf("\nValue:%s, Expected: %s\n", $id/$num, $expected);
	if( $val < $max) {
	    $max = $val;
	}
    }
    return $max;
}




#########################################################################
#
# 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;

################
#
# Net::MPEG-TS::Packet: MPEG-TS packet object (perl only)
#
# UNINETT 2009 (www.uninett.no)
# Dabbler in code: Odd Rune S. Lykkebø <lykkebo@pvv.ntnu.no>
#
# References:
# [1] - ISO 13818-1

use strict;
use Carp; #The mother hen is watchful.
use Bit::Vector;

use vars qw/$VERSION/;
$VERSION="0.05";


sub new 
{
    my $class = shift;
    my ($packet) = @_;

    my $self = {
	valid => 1,
	size => undef,
	sync_byte => 0x47,
	transport_error => 0,
	transport_priority => 0,
	payload_unit_start => 0,
	PID => undef,
	transport_scrambling => 0,
	adaption_field_control => undef,
	continuity_counter => 0,

#Adaption fields
	adaption_field_length => 0,
# flags
	discontinuity => 0,
	random_access => 0,
	elementary_stream_priority => 0,
	PCR_flag => 0,   #ok
	OPCR_flag => 0, #ok
	splicing_point => 0,
	transport_private_data => 0,
	adaption_field_extension_flag => 0,

# optional adaption fields existence based on flags
	program_clock_reference_base => 0, #ok
	program_clock_reference_extension => 0, 
	PCR => 0,
	original_program_clock_reference_base => 0, #ok
	original_program_clock_reference_extension => 0, 
	splice_countdown => 0,
	transport_private_data_length => 0,
	private_data => 0,
	adaption_field_extension_length => 0,
# mu?
	ltw_flag => 0,
	piecewise_rate_flag => 0,
	seamless_splice_flag => 0,
	ltw_valid_flag => 0,
	ltw_offset => 0,
	piecewise_rate => 0,
	splice_type => 0,
	DTS => 0,

	payload => undef,

    };
    bless ($self, $class);

    my $res = undef;
    if (defined ($packet)) {
	$self->decode ($packet);
    } else {
#Ok, want to create a packet -- that's... easy ?
#so, stuff stuff inside and hope for merry xmas.
#mpeg-ts doesn't really have any complicated fields,
#should be okeydokey
    }

    return $self;
}

sub is_valid
{
    my $self = shift;
    return $self->{'valid'};
}

sub adaption_field_control
{
    my $self = shift;
    my ($afc) = @_;

    if (defined $afc and (
($afc != 1) or
($afc != 2) or 
($afc != 3))) {
	carp "Bad value to adaption field control!";
    }
$self->{'adaption_field_control'} =	$afc if (defined $afc);

return $self->{'adaption_field_control'};
}

sub sync_byte
{
    my $self = shift;
    my ($sync_byte) = @_;
    $self->{'sync_byte'} = $sync_byte if (defined $sync_byte);
    return $self->{'sync_byte'};
}

sub continuity_counter
{
    my $self = shift;
    my ($continuity_counter) = @_;
    $self->{'continuity_counter'} = $continuity_counter if (defined $continuity_counter);

    return $self->{'continuity_counter'};
}

sub compare_cc_dists
{
    my (%a, %b) = @_;
    my $metric = 0;
    foreach my $pid (keys %a) {
	my $diff = $a{$pid} - $b{$pid};
	my $sum = $a{$pid} + $b{$pid};
	$metric += $diff * $sum;
	$metric /= 10;
    }
}


sub program_clock_reference_base
{
    my $self = shift;
    my ($program_clock_reference_base) = @_;
    $self->{'program_clock_reference_base'} = $program_clock_reference_base if (defined $program_clock_reference_base);
    return $self->{'program_clock_reference_base'};
}


sub transport_error
{
    my $self = shift;
    my ($transport_error) = @_;
    $self->{'transport_error'} = $transport_error if (defined $transport_error);
    return $self->{'transport_error'};
}

sub payload_unit_start
{
    my $self = shift;
    my ($payload_unit_start) = @_;
    $self->{'payload_unit_start'} = $payload_unit_start if (defined $payload_unit_start);
    return $self->{'payload_unit_start'};
}

sub PID
{
    my $self = shift;
    my ($PID) = @_;
    $self->{'PID'} = $PID if (defined $PID);
    return $self->{'PID'};
}

sub pcr_flag
{
    my $self = shift;
    return $self->{'PCR_flag'};
}

sub opcr_flag
{
    my $self = shift;
    return $self->{'OPCR_flag'};
}


sub pcrb
{
    my $self = shift;
    my ($pcrb) = @_;
    if (defined $pcrb) {
	$self->{'program_clock_reference_base'} = $pcrb;
	$self->{'PCR_flag'} = 1;
    }
    return $self->{'program_clock_reference_base'};
}

sub pcrb_s
{
    my $self = shift;
    return $self->{'program_clock_reference_base'}/(90000);
}

sub pcr
{
    my $self = shift;
    return $self->{'program_clock_reference_base'}*300 + $self->{'program_clock_reference_extension'};
}

sub pcr_s
{
    my $self = shift;
    my $base = Math::BigFloat->new($self->{'program_clock_reference_base'})->bdiv(90000);
    return $base;
}


#TODO: Reed-Solomon fault tolerance stuff.
sub payload
{
    my $self = shift;
    my ($payload) = @_;
    if (defined $payload) {
	$self->{'payload'} = $payload;
    }
    return $self->{'payload'};
}

sub pp_header 
{
    my $self = shift;
    print "== MPEG TS HEADER ==\n";
    printf ("PID: %d\n".
	    "Transport priority: %d\n".
	    "Payload unit start: %d\n".
	    "Transport error: %d\n".
	    "Scrambling: %s\n".
	    "Adaption field: %s\n".
	    "Continuation counter: %s\n",
	    $self->{'PID'},
	    $self->{'transport_priority'},
	    $self->{'payload_unit_start'},
	    $self->{'transport_error'},
	    $self->{'transport_scrambling'},
	    $self->{'adaption_field_control'},
	    $self->{'continuity_counter'});
    print ":: == MPEG TS HEADER END == ::\n";
}

#pretty-print the adaption field
sub pp_adaption_field
{
    my $self = shift;
    if ($self->{'adaption_field_length'} < 1) {
	carp "No adaption field to print.\n";
	return;
    }
    print "== Adaption field ==\n";
    printf ("Length(bytes): %d", $self->{'adaption_field_length'});
    printf ("Discontinuity: %d".
	    "Random Access: %d".
	    "Priority: %d".
	    "PCR-flag: %d ".
	    "OPCR-flag: %d ".
	    "Splicing point: %d ".
	    "Transport Private data: %d".
	    "Adaption field extension: %d ",
	    $self->{'discontinuity'},
	    $self->{'random_access'}, 
	    $self->{'elementary_stream_priority'}, 
	    $self->{'PCR_flag'}, 
	    $self->{'OPCR_flag'}, 
	    $self->{'splicing_point_flag'}, 
	    $self->{'transport_private_data'},
	    $self->{'adaption_field_extension'});

    print ("= Adaption field extensions =\n");
    if ($self->{'PCR_flag'} == 1) {
	printf ("PCR: %d", 
		$self->{'program_clock_reference_base'});
    }

    if ($self->{'OPCR_flag'} == 1) {
	printf ("OPCR: 0x%x", 
		$self->{'original_program_clock_reference_base'});
    }
    print "\n";

}

#The eternal mess of 7 demons. 
sub decode 
{
    my $self = shift;
    my ($packet) = @_;

#Extract the length of the packet. Should be 188 or 204 bytes.
    $self->{'size'} = length ($packet);
    carp "Bad packet size ($self->{'size'})\n" if ($self->{'size'} > (204*8));

#Decode the standard 4-byte header (ISO/IEC 13818-1)
    my ($sync, 
	$flags_pid, 
	$scradaptcont, 
	$adapt_and_payload) = unpack ("C n C a*", $packet);
    if ($sync != 0x47) {
	carp "Wops: bad sync byte ($sync)!\n";
	$self->{'valid'} = 0;
    }

    $self->{'sync_byte'} = $sync;
#	$flags_pid = hex $flags_pid;
    $self->{'transport_error'} = ($flags_pid & 0x8000) >> 15;
    $self->{'payload_unit_start'} = ($flags_pid >> 14) & 0x0001;
    $self->{'transport_priority'} = ($flags_pid & 0x2000) >> 13;
    $self->{'PID'} = ($flags_pid & 0x1FFF);
    $self->{'transport_scrambling'} = ($scradaptcont & 0xC0) >> 6;
    $self->{'adaption_field_control'} = ($scradaptcont & 0x30) >> 4;
    $self->{'continuity_counter'} = ($scradaptcont & 0x0F);

    if ($self->{'adaption_field_control'} == 2 or
	($self->{'adaption_field_control'} == 3)) {
#decode adaption field
	$self->_decode_adaption_field ($adapt_and_payload);
    }

#only payload
    if ($self->{'adaption_field_control'} == 1){
	$self->{'payload'} = $adapt_and_payload;
	my ($pes_prefix, $stream_id, $pes_length, $baz) = unpack "H6 H2 H4 a*", $adapt_and_payload;
    }	

    elsif ($self->{'adaption_field_control'} == 3) {
	my ($adapt_length) = unpack ("C", $adapt_and_payload);
	my $l = $adapt_length + 1;
	$self->{'payload'} = unpack ("x$l a*", $adapt_and_payload);
    }
}

sub _decode_adaption_field
{
    my $self = shift;
    my ($bin) = @_;

#Check for presence of adaption field
    my ($length, 
	$flags, 
	$adaption_opts_and_payload) = unpack ("C B8 a*", $bin);

#Sanity checks, concerning the length of the adaption fields, page 21 in [1]
    my $adaptc = $self->{'adaption_field_control'};
    if ($adaptc == 2 and $length != 183) {
	carp "Adaption field of illegal size ($adaptc) ($length)\n";
	return;
    }

    if ($adaptc == 3 and $length > 182) { 
	carp "Adaption field of illegal size ($adaptc) ($length)\n";
	return;
    }

    return if ($length < 1);

    $self->{'adaption_field_length'} = $length;


    ($self->{'discontinuity'} ,
     $self->{'random_access'} ,
     $self->{'elementary_stream_priority'} ,
     $self->{'PCR_flag'} ,
     $self->{'OPCR_flag'} ,
     $self->{'splicing_point_flag'} ,
     $self->{'transport_private_data'} ,
     $self->{'adaption_field_extension'} ) = split(//, $flags);

#Based on flags, we fill in the adaption field structures.
    my $field_position = 0;
    my ($foo, $bar); #because perl doesn't allow me to skip bits
    if ($self->{'PCR_flag'} == 1) {
	
	my ($bs) = unpack ("B48", $adaption_opts_and_payload);
	my $base = substr ($bs, 0, 33);
	my $extension = substr ($bs, 39, 9);

	$self->{'program_clock_reference_base'} = Math::BigInt->new("0b".$base);
	$self->{'program_clock_reference_extension'} = oct ("0b".$extension);
	$field_position += 6;
    }


    if ($self->{'OPCR_flag'} == 1) {
	($self->{'original_program_clock_reference_base'}, $foo,
	 $self->{'original_program_clock_rererence_extension'},
	 $bar) = unpack ("x$field_position B33 B6 B9 a*", 
			 $adaption_opts_and_payload);
	$field_position += 6;	
    }

    return $self;

}

sub bin2dec 
{
    return unpack("N", pack("B32", substr("0" x 32 . shift, -32)));
}


return (1);












package mpegpespacket;

################
#
# Net::MPEG-TS::PES-Packet: 
#
# UNINETT 2009 (www.uninett.no)
# Dabbler in code: Odd Rune S. Lykkebø <lykkebo@pvv.ntnu.no>
#
# References:
# [1] - ISO 13818-1


use strict;
use Carp;
use Bit::Vector;

use vars qw/$VERSION/;
$VERSION="0.05";


use constant PES_PREFIX => 0x000001;

use constant PROGRAM_STREAM_MAP => 0xBC;
use constant PRIVATE_STREAM_1 => 0xBD;
use constant PADDING_STREAM => 0xBE;
use constant PRIVATE_STREAM_2 => 0xBF;
use constant ECM_STREAM => 0xF0;
use constant EMM_STREAM => 0xF1;
use constant DSMCC_STREAM => 0xF2;
use constant ITU_222_A => 0xF4; 
use constant ITU_222_B => 0xF5; 
use constant ITU_222_C => 0xF6; 
use constant ITU_222_D => 0xF7; 
use constant ITU_222_E => 0xF8; 
use constant ANCILLARY_STREAM => 0xF9;
use constant ISO_14496_PACKETIZED_STREAM => 0xFA;
use constant ISO_13396_FLEXMUX_STREAM => 0xFB;
use constant PROGRAM_STREAM_DIRECTORY => 0xFF;

# From [1]:
#      1011 1100       1 program_stream_map
#      1011 1101       2 private_stream_1
#      1011 1110         padding_stream
#      1011 1111       3 private_stream_2
#      110x xxxx         ISO/IEC 13818-3 or ISO/IEC 11172-3 or ISO/IEC 13818-7 or ISO/IEC
#                        14496-3 audio stream number x xxxx
#      1110 xxxx         ITU-T Rec. H.262 | ISO/IEC 13818-2 or ISO/IEC 11172-2 or ISO/IEC
#                        14496-2 video stream number xxxx
#      1111 0000       3 ECM_stream
#      1111 0001       3 EMM_stream
#      1111 0010       5 ITU-T Rec. H.222.0 | ISO/IEC 13818-1 Annex A or ISO/IEC 13818-
#                        6_DSMCC_stream
#      1111 0011       2 ISO/IEC_13522_stream
#      1111 0100       6 ITU-T Rec. H.222.1 type A
#      1111 0101       6 ITU-T Rec. H.222.1 type B
#      1111 0110       6 ITU-T Rec. H.222.1 type C
#      1111 0111       6 ITU-T Rec. H.222.1 type D
#      1111 1000       6 ITU-T Rec. H.222.1 type E
#      1111 1001       7 ancillary_stream
#      1111 1010         ISO/IEC14496-1_SL-packetized_stream
#      1111 1011         ISO/IEC14496-1_FlexMux_stream
#1111 1100 ... 1111 1110   reserved data stream
#      1111 1111       4 program_stream_directory

sub new
{
    my $class = shift;
    my ($packet) = @_;

    my $self = {
	valid => 0,
	packet_start_code_prefix => undef,
	stream_id => undef,
	packet_length => 0,
	scrambling_control => 0,
	priority => 0,
	data_alignment_indicator => 0,
	copyright => 0,
	original_or_copy => 1,
	PTS_flag => 0,
	DTS_flag => 0,
	ESCR_flag => 0,
	ES_rate_flag => 0,
	DSM_trick_mode_flag => 0,
	additional_copy_info_flag => 0,
	CRC_flag => 0,
	extension_flag => 0,
	header_data_length => 0,
	marker_bit => 0,
	PTS => undef,
	DTS => undef,
	ESCR_base => undef,
	ESCR_extension => undef,
	ES_rate => undef,
	trick_mode_control => undef,
	slow_motion => undef,
	freeze_frame => undef,
	fast_reverse => undef,
	slow_reverse => undef,
	field_id => undef,
	rep_cntrl => undef,
	additional_copy_info => undef,
	previous_PES_packet_CRC => undef,
	PES_private_data_flag => 0,
	pack_header_field_flag => 0,
	program_packet_sequence_counter_flag => 0,
	P_STD_buffer_flag => 0,
	PES_extension_flag_2 => 0,
	PES_private_data => undef,
	pack_field_length => undef,
	program_packet_sequence_counter => undef,
	MPEG1_MPEG2_identifier => undef,
	original_stuff_length => undef,
	P_STD_buffer_scale => undef,
	P_STD_buffer_size => undef,
	PES_extension_field_length => undef,
	stuffing_byte => undef,
	PES_packet_data_byte => undef,
	padding_byte => 0xFF,   #wtf

	payload => undef, 
    }; 
    bless ($self, $class);

    if (defined ($packet)) {
	$self->decode ($packet);
    } else {
#TODO: Encode packet :-) 
    }

    return $self;
}

sub ppme
{
    my $self = shift;
    if ($self->{'packet_start_code_prefix'} != PES_PREFIX) {
	print "NOT PES PACKET HEADER:";
	print $self->{'packet_start_code_prefix'};
	print "\n";
	return;
    }
    printf(	"- PES -\n".
		"Prefix: 0x%x\n".
		"Length: %s\n".
		"Stream-id: 0x%x (%s)\n".
		"Optional sections:\n".
		"pes_scrambling_control: %s\n".
		"pes_priority: %s\n".
		"data_aligmnent_indicator: %s\n".
		"copyright: %s\n".
		"original_or_copy: %s\n".
		"pts_flag: %s\n".
		"dts_flag: %s\n".
		"escr_flag: %s\n".
		"es_rate_flag: %s\n".
		"dsm_trick_mode_flag: %s\n".
		"additional_copy_info_flag: %s\n".
		"crc_flag: %s\n".
		"extension_flag: %s\n".
		"header_data_length: %s\n".
		"ES-rate (bytes/s): %d\n".
		"PTS: %f\n".
		"- END PES -\n",
		$self->{'packet_start_code_prefix'}, 
		$self->{'packet_length'},
		$self->{'stream_id'},
		$self->_stream_id_type ($self->{'stream_id'}),
		$self->{'PES_scrambling_control'} ,
		$self->{'PES_priority'} ,
		$self->{'data_aligmnent_indicator'} ,
		$self->{'copyright'} ,
		$self->{'original_or_copy'} ,
		$self->{'PTS_flag'} ,
		$self->{'DTS_flag'} ,
		$self->{'ESCR_flag'} ,
		$self->{'ES_rate_flag'} ,
		$self->{'DSM_trick_mode_flag'} ,
		$self->{'additional_copy_info_flag'} ,
		$self->{'CRC_flag'} ,
		$self->{'extension_flag'} ,
		$self->{'header_data_length'},
		50 * $self->{'ES_rate'},
		(300*$self->{'PTS'})/(27*(10**6))
		);
}


sub is_valid
{
    my $self = shift;
    return $self->{'valid'};
}


sub payload
{
    my $self = shift;
    my ($payload) = @_;
    if (defined $payload) {
	$self->{'payload'} = $payload;
    }
    return $self->{'payload'};
}



sub _stream_id_type
{
    my $self = shift;

    if ($self->{'stream_id'} >= 0xe0 and $self->{'stream_id'} <= 0xef) {
	return "video";
    } elsif ($self->{'stream_id'} >= 0xc0 and $self->{'stream_id'} <= 0xcf) {
	return "audio";
    } elsif ($self->{'stream_id'} == PRIVATE_STREAM_1) { return "private_stream_1"; }	
    elsif ($self->{'stream_id'} == PROGRAM_STREAM_MAP) { return "program_stream_map"; }	
    else {
	return "unk";
    }
}


sub decode 
{
    my $self = shift;
    my ($packet) = @_;

    my $flags_data = undef;	
    my $prefix = undef;
    ($prefix,		
     $self->{'stream_id'},
     $self->{'packet_length'},
     $flags_data) = unpack ("H6 C H4 a*", $packet);

    $prefix = hex $prefix;
    if ($prefix != PES_PREFIX) {
#carp "This is not a PES-packet: $prefix $self->{'packet_start_code_prefix'}";
	$self->{'valid'} = 0;
	return;
    } else {
	$self->{'packet_start_code_prefix'} = $prefix;
	$self->{'valid'} = 1;
    }

    if ($self->{'stream_id'} != PROGRAM_STREAM_MAP and
	$self->{'stream_id'} != PADDING_STREAM and
	$self->{'stream_id'} != PRIVATE_STREAM_2 and 
	$self->{'stream_id'} != ECM_STREAM and
	$self->{'stream_id'} != EMM_STREAM and
	$self->{'stream_id'} != PROGRAM_STREAM_DIRECTORY and
	$self->{'stream_id'} != DSMCC_STREAM and
	$self->{'stream_id'} != ITU_222_E) {

# The code below does NOT do what you expect it to. There are subtle issues with un
#packing these things. do not be fooled. It does not work.
#
#			my ($flags, $data) = unpack ("a3 a*", $flags_data);
#			my ($weird,$data) = undef;
#			($weird, 
#			$self->{'PES_scrambling_control'} ,
#			$self->{'PES_priority'} ,
#			$self->{'data_aligmnent_indicator'} ,
#			$self->{'copyright'} ,
#			$self->{'original_or_copy'} ,
#			$self->{'PTS_DTS_flags'} ,
#			$self->{'ESCR_flag'} ,
#			$self->{'ES_rate_flag'} ,
#			$self->{'DSM_trick_mode_flag'} ,
#			$self->{'additional_copy_info_flag'} ,
#			$self->{'CRC_flag'} ,
#			$self->{'extension_flag'} ,
#			$self->{'header_data_length'},
#			$data) = unpack ("B2 B2 B B B B B2 B B B B B B C a*", $flags_data);

#This, however -- does work.
	my @header = unpack("C3", $flags_data);
	my ($data) = unpack("x3 a*", $flags_data);

	my $control_bits = ($header[0] >> 6) & 0x03;
	$self->{'PES_scrambling_control'} = ($header[0] >> 4) & 0x03;
	$self->{'PTS_flag'} = ($header[1] >> 7) & 0x01;
	$self->{'DTS_flag'} = ($header[1] >> 6) & 0x01;
	$self->{'ES_rate_flag'} = ($header[1] >> 4) & 0x01;  #last bit of first nibble in second byte

#			printf("%x\n", $header[1]);
	carp "Weird and confusing (i'm scared mom), no '10($control_bits)' where it should've been" if ($control_bits != 0x2); 

#pull out the payload while we're at it.
	my $header_length = ($header[2]);
	($self->{'payload'}) = unpack("x$header_length a*", $data);

	my $data_pos = 0;
	if ($self->{'PTS_flag'} == 1 and $self->{'DTS_flag'} == 0){
	    $self->_decode_pts($data);
	    $data_pos += 5;
	} elsif ($self->{'PTS_flag'} == 1 and $self->{'DTS_flag'} == 1){
	    $self->_decode_pts_dts($data);
	    $data_pos += 10;
	}

	if ($self->{'ECSR_flag'} eq '1') {
	    my ($ecsr) = unpack ("x$data_pos a*", $data);
	    $self->_decode_ecsr ($ecsr);
	    $data_pos += 6;
	}

	if ($self->{'ES_rate_flag'} eq '1') {
	    my ($es) = unpack ("x$data_pos a*", $data);
	    $self->_decode_es_rate ($es);
	    $data_pos += 3;
	}

	if ($self->{'DSM_trick_mode_flag'} == '1') {
	    my ($trick) = unpack ("x$data_pos a*", $data);
	    $self->_decode_dsm_trick_mode ($trick);
	    $data_pos += 1;	
	}

	if ($self->{'additional_copy_info_flag'} == '1') {
	    my ($addcop) = unpack ("x$data_pos a*", $data);
	    $self->_decode_additional_copy_info ($addcop);
	    $data_pos += 1;	
	}

	if ($self->{'CRC_flag'} == '1') {
	    my ($crc) = unpack ("x$data_pos a*", $data);
	    $self->_decode_CRC ($crc);
	    $data_pos += 2;	
	}

	if ($self->{'extension_flag'} == '1') {
	    my ($ext) = unpack ("x$data_pos a*", $data);
	    $self->_decode_pes_extension ($ext);
#end.
	}
#stuffing bytes and stuff goes here 
    } elsif (
	     $self->{'stream_id'} == PROGRAM_STREAM_MAP or
	     $self->{'stream_id'} == PRIVATE_STREAM_2 or 
	     $self->{'stream_id'} == ECM_STREAM or
	     $self->{'stream_id'} == EMM_STREAM or
	     $self->{'stream_id'} == PROGRAM_STREAM_DIRECTORY or
	     $self->{'stream_id'} == DSMCC_STREAM or
	     $self->{'stream_id'} == ITU_222_E) { 

	$self->{'data'} = unpack ("x3 a*", $flags_data);		
    } elsif ($self->{'stream_id'} == PADDING_STREAM) {
#this isn't very interesting.
    }
}
sub pts_s
{
    my $self = shift;
    my $x = Math::BigFloat->new($self->{'PTS'})/(90000);
    return $x;
}

#This function is BROKEN since bit-unpacking does not behave like expected
sub _decode_es_rate
{
    my $self = shift;
    my ($data) = @_;
    my ($mark, $mark2, $mark3) = undef;
#	print "i shouldn't be here!\n";
    ($mark,
     $mark2,  
     $mark3) = unpack ("B B22 B", $data);
    if ($mark != '1' or $mark3 != '1'){
#		print "oh really\n"22;
	carp "ES_rate marker bits not 1";
    }
    $self->{'ES_rate'} = bin2dec ($mark2);
    return $self->{'ES_rate'};
}


#Unfinished
sub _decode_pts
{
    my $self = shift;
    my ($data) = @_;

    my (@pts) = unpack ("C5", $data);
    my ($p1) = unpack ("H10", $data);
    my $check = ($pts[0] >> 4) & 0x0F;

#	my ($pts) = unpack ("B40", $data);
#	my $check = substr($pts, 0, 4);
    if ($check != 0b0010) {
	carp "PTS decoding error. Must-be 0b0010 was $check, dropping packet\n";
	return;
    }

#The 9 next nibbles encodes the PTS, 
#but there's marker bits so we need to shift the data a bit.
# 
#1 byte			2 byte			   		  2 byte
#aaaa xxxb  xxxxxxxxxxxxxxxxxxb xxxxxxxxxxxxxxxxxxb

    my ($bits) = unpack ("B40", $data);
    my $p1 = substr($bits, 4, 3);
    my $p2 = substr($bits, 8, 15);
    my $p3 = substr($bits, 24, 15);
#
#	my $pesv = Bit::Vector->new(40);
#	$pesv->from_Bin($bits);
#
#	$pesv->Delete(36,4);
#	$pesv->Delete(32,1);
#	$pesv->Delete(15,1);
#	$pesv->Delete(0,1);

#	my $pts_decoded = $pesv->to_Dec();
    $self->{'PTS'} = Math::BigInt->new("0b".$p1.$p2.$p3);
}
#TODO:
sub _decode_pts_dts {
#todo!
}
sub _decode_ecsr {}
sub _decode_dsm_trick_mode {}
sub _decode_additional_copy_info {}
sub _decode_CRC {}
sub _decode_pes_extension {}


sub bin2dec 
{
    return unpack("N", pack("B32", substr("0" x 32 . shift, -32)));
}



return (1);

package mpegespacket;

################
#
# Net::MPEG-TS::ES-packet
#
# UNINETT 2009 (www.uninett.no)
# Dabbler in code: Odd Rune S. Lykkebø <lykkebo@pvv.ntnu.no>
#
# References:
# [1] - Wikipedia... 


use strict;
use Carp;

use vars qw/$VERSION/;
$VERSION="0.05";


sub new
{
    my $class = shift;
    my ($packet) = @_;

    my $self = {
	valid => 0,
	start_code => 0x000001B3,
	horizontal_size => undef,
	vertical_size => undef,
	aspect_ratio => undef,
	frame_rate => undef,
	bit_rate => undef,
	rest => undef, # :-)
    }; 
    bless ($self, $class);

    if (defined ($packet)) {
	$self->decode ($packet);
    } else {
#TODO: Encode packet :-) 
    }

    return $self;
}

sub ppme
{
    my $self = shift;

    if($self->{'start_code'} != 0x000001b3) {
	print "NO HEADER IN ELEMENTAL STREAM\n";
	return;
    } 
    printf(	"- ES -\n".
		"Start-code: 0x%x\n".
		"Horizontal: %d\n".
		"Vertical: %d\n".
		"Aspect: %d\n".
		"Frame-rate(fps): %d\n".
		"Byte-rate(actual byte/s): %d\n".
		"- END ES -\n",

		$self->{'start_code'},
		$self->{'horizontal_size'},
		$self->{'vertical_size'},
		$self->{'aspect_ratio'},
		$self->_frame_rate_code($self->{'frame_rate'}),
		(400/8)*$self->{'bit_rate'} );
}

sub payload
{
    my $self = shift;
    my ($payload) = @_;
    if (defined $payload) {
	$self->{'payload'} = $payload;
    }
    return $self->{'payload'};
}


sub is_valid
{
    my $self = shift;
    return $self->{'valid'};
}


sub decode 
{
    my $self = shift;
    my ($packet) = @_;

    my ($hor,$ver,$asp,$fr,$br);
    ($self->{'start_code'},
     $hor,
     $ver,
     $asp,
     $fr,
     $br) = unpack ("N H3 H3 H H N", $packet);

    $self->{'horizontal_size'} = hex $hor;
    $self->{'vertical_size'} = hex $ver;
    $self->{'aspect_ratio'} = hex $asp;
    $self->{'frame_rate'} = hex $fr;
    $self->{'bit_rate'} = ($br >> 14) & 0x0003FFFF;

    if($self->{'start_code'} != 0x000001B3) {
	carp "Not a ES-packet header ($self->{'start_code'})";
	$self->{'valid'} = 0;
    } else {
	$self->{'valid'} = 1;
    }
}


sub _frame_rate_code
{
    my $code = @_;
    return 30 if $code == 5;
    return 24 if $code == 2;
}
sub bin2dec 
{
    return unpack("N", pack("B32", substr("0" x 32 . shift, -32)));
}



return (1);