From 38001c28d981fed7853972af541d48c7830dfa5d Mon Sep 17 00:00:00 2001 From: Andrew Elder Date: Fri, 19 Jul 2013 15:27:27 -0400 Subject: [PATCH] test: add python scripts for extracting and analyzing AS timestamps in 1722 AVTP packets --- test/README.rst | 60 ++++++++++++++++ test/astime_fitline.py | 185 +++++++++++++++++++++++++++++++++++++++++++++++++ test/avtp_astimes.py | 165 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 410 insertions(+) create mode 100644 test/README.rst create mode 100644 test/astime_fitline.py create mode 100644 test/avtp_astimes.py diff --git a/test/README.rst b/test/README.rst new file mode 100644 index 0000000..9f06862 --- /dev/null +++ b/test/README.rst @@ -0,0 +1,60 @@ +Test tools +========== + +This directory contains AVB related tools + +AVTP Timestamps +--------------- + +This section covers tools and scripts related to AVTP timestamp analysis. + +Extracting timestamps +..................... + +The python script *avtp_astimes.py* extracts AVTP AS presentation times from +AVTP packets in a Wireshark capture file. The input capture file should +be stored in libpcap format. If there is more than one active AVB stream +in the packet capture, the python script can be run twice to extract multiple +timestamp sequences. + +The below operation +:: + $python avtp_astimes.py -c 60000 -m 91:e0:f0:00:15:20 capture.libpcap a.csv + $python avtp_astimes.py -c 60000 -m 91:e0:f0:00:35:80 capture.libpcap b.csv + +extracts 60,000 AS presentation times from capture.libpcap and stores the times +in a.csv and b.csv. Recall that for an AVB stream running with a 48 kHz sample +rate, AS timestamps are inserted every 8 samples, making 6,000 AS timestamps +per second. The count of 60000 timestamps would therefore correspond to 10 seconds +of audio. The -m option specifies the destination multicast MAC that the AVB stream +is being sent to. The first three octets of the MAC address are 91:e0:f0 as reserved +by IEEE to 1722 AVTP transport. The destination MAC address is determined by manually +reading the packet capture and determining what destination MAC addresses are of interest. + +Note that avtp_astimes.py does not check for missing packets and does not use the +packet capture timestamp at all. Additionally avtp_astimes.py "unwraps" the AS times +so that the output times are monotonically increasing. + +Tested using python 2.6.5. + +Line fitting +............ + +The python script *astime_fitline.py* fits a line to a .csv file that contains +AVTP presentation times. Least mean square error is used to perform the line fit. +The script can optionally fit two lines to two separate .csv sequences +and output the difference in slope and y intercept between the two lines. + +The below operation +:: + $python avtp_fitline.py -c 60000 a.csv b.csv + +Fits lines to the a.csv AS presentation sequence and the b.csv presentation time +sequence. The first 60,000 AS presentation times are used for both lines. + +Use the -p option to plot the distrubition of AS time deviations from the fitted line. +The FFT of the deviations is also plotted if the count option is greater than or equal +to 32768. + +Tested using python 2.6.5. + diff --git a/test/astime_fitline.py b/test/astime_fitline.py new file mode 100644 index 0000000..2858f7d --- /dev/null +++ b/test/astime_fitline.py @@ -0,0 +1,185 @@ +#!/usr/bin/env python +# +############################################################################## +# +# Copyright (c) 2013, AudioScience, Inc +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# 3. Neither the name of the Intel Corporation nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +############################################################################## +# +# Fits a line to a sequence of AVTP AS timestamps stored +# in a csv format ascii file. +# + + + +usage = """astime_fitline.py -c 100 ts_file.txt +""" + +from sys import exit, platform +from time import sleep +from struct import unpack +import matplotlib.pyplot as plt +import scipy +import scipy.fftpack +import numpy as np +import csv + + +def fit_line(count, plots, file1): + # see http://en.wikipedia.org/wiki/Simple_linear_regression + SX = 0.0 + SY = 0.0 + SXX = 0.0 + SXY = 0.0 + n = 0.0 + res = {} + with open(file1, 'rU') as csvfile: + vals = list(csv.reader(csvfile)) + for s in vals[0:count]: + x = float(s[0]) + y = float(s[1]) + SX = SX + x + SXX = SXX + (x * x) + SY = SY + y + SXY = SXY + (x * y) + n = n + 1 + m = (n * SXY - SX * SY)/(n * SXX - SX * SX) + b = SY/n - m * SX/n + x_vals = np.array([float(vals[i][0]) for i in range(count)]) + y_line_array = x_vals * m + b + y_array = np.array([float(vals[i][1]) for i in range(count)]) + line_deviation = y_line_array - y_array + max_error = np.amax(abs(line_deviation)) + rnd_max_error = np.int_(max_error) + histogram_bins = range((rnd_max_error + 1) * 2) - (rnd_max_error + 1); + hist, bins = np.histogram(line_deviation, histogram_bins) + res['histogram'] = hist + res['bins'] = bins + res['file'] = file1 + res['line_slope'] = m + res['line_y_intercept'] = b + res['line_max deviataion'] = max_error + res['astimestamps'] = y_array + if plots == True: + if True: # Distribution bar graph, switch to True to run + figh = plt.figure() + center = (bins[:-1]+bins[1:])/2 + plt.bar(center, hist, align = 'center', width = 0.7) + plt.xlabel('deviation from fitted line in ns') + plt.ylabel('count') + plt.title(file1 + ' AS timestamp distribution\nof deviations from fitted line') + res['hist_fig'] = figh + + if False: # Fitted line plot, switch to True to run + figl = plt.figure() + plt.plot(range(count), y_array, 'ro', range(count), y_line_array, 'b-') + plt.xlabel('AS timestamp count') + plt.ylabel('AS timestamp') + plt.title(file1 + ' fitted line') + res['line_fig'] = figl + + if True: # Deviation from fitted line plot, switch to True to run + figl = plt.figure() + plt.plot(range(count), line_deviation, 'ro') + plt.xlabel('AS timestamp count') + plt.ylabel('AS timestamp deviation') + plt.title(file1 + ' deviation from fitted line') + res['line_err'] = figl + + if count > 32768: + if False: # Frequency transform of deviation from fitted line, switch to True to run + figf = plt.figure() + fft_count = 32678 + # fft bin spacing is Fs/N, where Fs = 48000/8 and N = fft_count + fft_bin_spacing = float(48000/8) / float(fft_count) + freq_bins = [float(i) * fft_bin_spacing for i in range(fft_count/2 + 1)] + FFT = abs(np.fft.rfft(line_deviation[range(fft_count)])) + plt.plot(freq_bins, np.log10(FFT)) + plt.xlabel('Hz') + plt.title(file1 + ' deviation periodicity') + res['fft_fig'] = figf + + csvfile.close() + return res + +def print_line(res): + print 'File: ' + res['file'] + print 'Line slope: ' + str(res['line_slope']) + print 'Line y intercept: ' + str(res['line_y_intercept']) + print 'Maximum deviation from fitted line (ns): ' + str(res['line_max deviataion']) + print 'Histgram of deviations (ns): ' + print '\tBins (ns): ' + print res['bins'] + print '\tCounts: ' + print res['histogram'] + +if __name__ == '__main__': + from optparse import OptionParser + + parser = OptionParser(usage = usage) + parser.add_option('-c','--count',type='int',dest='count', + help='Number of sample points to fit. Default=%default', default=100) + parser.add_option('-p', action='store_true', dest='plot', + help='Plot lines, histogram and fft. Default=%default', default=False) + + opts,args = parser.parse_args() + + if len(args) == 0: + parser.print_usage() + exit() + # compute + result1 = fit_line(opts.count, opts.plot, args[0]) + if len(args) == 2: + result2 = fit_line(opts.count, opts.plot, args[1]) + + #plot + if opts.plot: + plt.show() + + # print + print_line(result1) + if len(args) == 2: + print_line(result2) + print "=== The following output compares fitted lines and data between the two input data sets ===" + ppb_freq = abs(result1['line_slope'] - result2['line_slope'])/result1['line_slope'] * 1000000000 + print "PPM freq difference is: %f" % ppb_freq + ASTime_diff = abs(result1['line_y_intercept'] - result2['line_y_intercept']) + time_diff_in_samples = ASTime_diff * 48000 / 1000000000 + print "AStime offset in is %f ns, or %f samples" % (ASTime_diff, time_diff_in_samples) + floor_time_diff = np.int(time_diff_in_samples + 0.5) + adj_time_diff = float(floor_time_diff) * 1000000000 / 48000 + print "AStime offset rounded to nearest sample is %d samples, or %f ns" % (floor_time_diff, adj_time_diff) + if result1['line_y_intercept'] < result2['line_y_intercept']: + adj_time_diff = adj_time_diff * -1 + time_astimestamps_diff = np.amax(abs((result1['astimestamps'] - result2['astimestamps']) - adj_time_diff)) + print "Maximum instaneous offset per AS timestamp is %f ns (%f of sample period)" % \ + (time_astimestamps_diff, time_astimestamps_diff / (1000000000 / 48000)) + print "\tfor reference 0.25 of a sample at 48 kHz is %f ns" % ( 10000000000 / 48000 / 4) + diff --git a/test/avtp_astimes.py b/test/avtp_astimes.py new file mode 100644 index 0000000..633a96d --- /dev/null +++ b/test/avtp_astimes.py @@ -0,0 +1,165 @@ +#!/usr/bin/env python +# +############################################################################## +# +# Copyright (c) 2013, AudioScience, Inc +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, +# this list of conditions and the following disclaimer. +# +# 2. Redistributions in binary form must reproduce the above copyright +# notice, this list of conditions and the following disclaimer in the +# documentation and/or other materials provided with the distribution. +# +# 3. Neither the name of the Intel Corporation nor the names of its +# contributors may be used to endorse or promote products derived from +# this software without specific prior written permission. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +# POSSIBILITY OF SUCH DAMAGE. +# +############################################################################## +# +# Extracts AVTP AS timestamps from a libpcap captre file and stores them in +# a csv format ascii file. +# + +import math +import scapy.all as s + +usage = """avtp_astimes.py [options] input.libpcap output.csv + +This script extracts AVTP packet timestamps and writes them +in a text file for later processing. Timestamps are modified +so as to be always incrementing (ie the 4s wrap is removed). + +""" + + +wraps = 0 +prev_pkt_ts = 0 +ts_count = 0 +seq = {} + +class AVTP(s.Packet): + name = "AVTP" + fields_desc=[ + s.XByteField("controlData", None), + s.XByteField("flags", None), + s.XByteField("sequence", None), + s.XByteField("timestampUncertain", None), + s.IntField("streamID0", None), + s.IntField("streamID1", None), + s.IntField("ptpTimestamp", None), + s.IntField("gateway", None), + s.ShortField("pktDataLength", None), + s.XByteField("pkt1394format", None), + s.XByteField("pkt1394tcode", None), + s.XByteField("sourceId", None), + s.XByteField("dataBlockSize", None), + s.XByteField("packing", None), + s.XByteField("DBC", None), + s.XByteField("formatId", None), + s.XByteField("SYT", None), + s.IntField("ptpUpper", None), + s.IntField("ptpLower", None), + ] +s.bind_layers(s.Ether, AVTP, type=0x22f0) + +def pkt_avtp(pkt, fout, pkt_count): + global wraps, prev_pkt_ts, ts_count, ts_accum, seq + + n = 0 + avtp = pkt[AVTP] + if avtp.controlData == 0x0: + if seq['init'] == False: + seq['last'] = avtp.sequence + seq['init'] = True + else: + if avtp.sequence < seq['last']: + if avtp.sequence != 0: + print "Sequence wrap error at packet number %d" % (pkt_count) + else: + if avtp.sequence - seq['last'] != 1: + print "Sequence error at packet number %d" % (pkt_count) + seq['last'] = avtp.sequence + + if avtp.flags == 0x81: + if ts_count == 0: + ts_accum = avtp.ptpTimestamp + else: + if avtp.ptpTimestamp > prev_pkt_ts: + ts_delta = avtp.ptpTimestamp - prev_pkt_ts + else: + ts_delta = avtp.ptpTimestamp + 0x100000000 - prev_pkt_ts + ts_accum = ts_accum + ts_delta + fout.write("%d, %d\n" % (ts_count, ts_accum)) + n = 1 + prev_pkt_ts = avtp.ptpTimestamp + ts_count = ts_count + 1 + return n + +def main(count, input_cap, output_txt, mac_filter): + mac_counts = {} + pkt_count = 0 + seq['last'] = 0 + seq['init'] = False + capture = s.rdpcap(input_cap) + foutput = open(output_txt, 'wt') + for pkt in capture: + n = 0; + pkt_count += 1 + try: + if pkt[s.Ether].dst in mac_counts: + mac_counts[pkt[s.Ether].dst] = mac_counts[pkt[s.Ether].dst] + 1 + else: + mac_counts[pkt[s.Ether].dst] = 1 + + # Deal with AVTP packets + if pkt[s.Ether].type == 0x22f0: + # look for the requested MAC + if pkt[s.Ether].dst == mac_filter: + n = pkt_avtp(pkt, foutput, pkt_count) + except IndexError: + print "Unknown ethernet type" + count = count - n + if count == 0: + break + foutput.close(); + if count != 0: + print "Could not find the specified MAC, or MAC count" + print "Mac counts" + print mac_counts + print "Complete" + +if __name__ == "__main__": + from optparse import OptionParser + + parser = OptionParser(usage = usage) + parser.add_option('-m','--mac',type='string',dest='dst_mac', + help='Destination MAC address of the AVTP stream', default=None) + parser.add_option('-c','--count',type='int',dest='count', + help='Number of 802.1AS timestamps to extract. Default=%default', default=100) + + opts,args = parser.parse_args() + + print "Search for %d timestamped AVTP packets sent to dst MAC %s" % (opts.count, opts.dst_mac) + + if len(args) == 0: + parser.print_usage() + exit() + + main(opts.count, args[0], args[1], opts.dst_mac) -- 2.7.4