HP3562A/read_trace.py

# -*- coding: utf-8 -*-

import struct
import numpy


def decode_float(s):
    assert len(s) in [4,8]
    # exponential term 
    e = ord(s[-1])
    if e & 0x80:
        e = e - 256

    # mantissa
    m = [ord(x) for x in s[:-1]]
    M = 0.
    for i in range(len(s)-1):
        #M += m[i]<<(i*8)
        M += float(m[i])/2**((i+1)*8)
    # XXX how do we deal negative numbers?
    #if m[0] & 0x80:
    #    M = M - 2^(len(s))
    return M * 2**(e+1)

def decode_string(s):
    nb = ord(s[0])
    s = s[1:nb+2]
    r = ""
    # XXX why do we need to do this? It's not described in the manual...
    for c in s:
        r += chr(ord(c) & 0x7F)
    return r

EDSP = {0: "No data",
        1: "Frequency response",
        2: "Power spectrum 1",
        3: "Power spectrum 2",
        4: "Coherence",
        5: "Cross spectrum",
        6: "Input time 1",
        7: "Input time 2",
        8: "Input linear spectrum 1",
        9: "Input linear spectrum 2",
        10: "Impulse response",
        11: "Cross correlation",
        12: "Auto correlation 1",
        13: "Auto correlation 2",
        14: "Histogram 1",
        15: "Histogram 2",
        16: "Cumulative density function 1",
        17: "Cumulative density function 2",
        18: "Probability density function 1",
        19: "Probability density function 2",
        20: "Average linear spectrum 1",
        21: "Average linear spectrum 2",
        22: "Average time record 1",
        23: "Average time record 2",
        24: "Synthesis pole-zeros",
        25: "Synthesis pole-residue",
        26: "Synthesis polynomial",
        27: "Synthesis constant",
        28: "Windowed time record 1",
        29: "Windowed time record 2",
        30: "Windowed linear spectrum 1",
        31: "Windowed linear spectrum 2",
        32: "Filtered time record 1",
        33: "Filtered time record 2",
        34: "Filtered linear spectrum 1",
        35: "Filtered linear spectrum 2",
        36: "Time capture buffer",
        37: "Captured linear spectrum",
        38: "Captured time record",
        39: "Throughput time record 1",
        40: "Throughput time record 2",
        41: "Curve fit",
        42: "Weighted function",
        43: "Not used",
        44: "Orbits",
        45: "Demodulation polar",
        46: "Preview demod record 1",
        47: "Preview demod record 2",
        48: "Preview demod linear spectrum 1",
        49: "Preview demod linear spectrum 2",
        }

ECH = {0: "Channel 1",
       1: "Channel 2",
       2: "Channel 1&2",
       3: "No channel",
       }

EOVR = ECH

EDOM = {0: 'Time',
        1: 'Frequency',
        2: 'Voltage (amplitude)',
        }

EVLT = {0: "Peak",
        1: "RMS",
        2: "Volt (indicates peak only)",
        }

EAMP = {0: "Volts",
        1: "Volts squared",
        2: "PSD (V²/Hz)",
        3: "ESD (V²s/Hz)",
        4: "PSD¹² (V/Hz¹²)",
        5: "No unit",
        6: "Unit volts",
        7: "Unit volts²",
        }
        
EXAXIS= {0: "No units",
         1: "Hertz",
         2: "RPM",
         3: "Orders",
         4: "Seconds",
         5: "Revs",
         6: "Degrees",
         7: "dB",
         8: "dBV",
         9: "Volts",
         10: "V Hz¹²",
         11: "Hz/s",
         12: "V/EU",
         13: "Vrms",
         14: "V²/Hz",
         15: "%",
         16: "Points",
         17: "Records",
         18: "Ohms",
         19: "Hertz/octave",
         20: "Pulse/Rev",
         21: "Decades",
         22: "Minutes",
         23: "V²s/Hz",
         24: "Octave",
         25: "Seconds/Decade",
         26: "Seconds/Octave",
         27: "Hz/Point",
         28: "Points/Sweep",
         29: "Points/Decade",
         30: "Points/Octave",
         31: "V/Vrms",
         32: "V²",
         33: "EU referenced to chan 1",
         34: "EU referenced to chan 2",
         35: "EU value",
         }

EMEAS = {0: "Linear resolution",
         1: "Log resolution",
         2: "Swept sine",
         3: "Time capture",
         4: "Linear resolution throughput",
         }

EDEMOD1 = {45: "AM",
           46: "FM",
           47: "PM",
           }

EDEMOD2 = EDEMOD1

EAVG = {0: "No data",
        1: "Not averaged",
        2: "Averaged",}


EWIN = {0: "N/A",
        1: "Hann",
        2: "Flat top",
        3: "Uniforme",
        4: "Exponential",
        5: "Force",
        6: "Force chan 1/expon chan 2",
        7: "Expon chan 1/force chan 2",
        8: "User",
        }

HEADER = [ ("Display function", EDSP, 'h', 2),
           ('Number of elements', int, 'h', 2),
           ('Displayed elements', int, 'h', 2),
           ('Number of averages', int, 'h', 2),
           ('Channel selection', ECH, 'h', 2),
           ('Overflow status', EOVR, 'h', 2),
           ('Overlap percentage', int, 'h', 2),
           ('Domain', EDOM, 'h', 2),
           ('Volts peak/rms', EVLT, 'h', 2),
           ('Amplitude units', EAMP, 'h', 2),
           ('X axis units', EXAXIS, 'h', 2),
           ('Auto math label', str, 's', 14),
           ('Trace label', str, 's', 22),
           ('EU label 1', str, 's', 6),
           ('EU label 2', str, 's', 6),
           ('Float/Interger', bool, 'h', 2),
           ('Complex/Real', bool, 'h', 2),
           ('Live/Recalled', bool, 'h', 2),
           ('Math result', bool, 'h', 2),
           ('Real/Complex input', bool, 'h', 2),
           ('Log/Linear data', bool, 'h', 2),
           ('Auto math', bool, 'h', 2),
           ('Real time status', bool, 'h', 2),
           ('Measurement mode', EMEAS, 'h', 2),
           ('Window', EWIN, 'h', 2),
           ('Demod type channel 1', EDEMOD1, 'h', 2),
           ('Demod type channel 2', EDEMOD2, 'h', 2),
           ('Demod active channel 1', bool, 'h', 2),
           ('Demod active channel 2', bool, 'h', 2),
           ('Average status', EAVG, 'h', 2),
           ('Not used', int, 'hh', 4),
           ('Samp freq/2 (real)', decode_float, None, 4),
           ('Samp freq/2 (imag)', decode_float, None, 4),
           ('Not used', decode_float, None, 4),
           ('Delta X-axis', decode_float, None, 4),
           ('Max range', decode_float, None, 4),
           ('Start time value', decode_float, None, 4),
           ('Expon wind const 1', decode_float, None, 4),
           ('Expon wind const 2', decode_float, None, 4),
           ('EU value chan 1', decode_float, None, 4),
           ('EU value chan 2', decode_float, None, 4),
           ('Trig delay chan 1', decode_float, None, 4),
           ('Trig delay chan 2', decode_float, None, 4),
           ('Start freq value', decode_float, None, 8),
           ('Start data value', decode_float, None, 8),
           ]

def decode_trace(data):    
    d = data

    typ = d[:2]
    assert typ == "#A"

    totlen = struct.unpack('>h', d[2:4])[0]
    idx = 4
    tt=0
    header = {}
    for i, (nam, dtype, fmt, nbytes) in enumerate(HEADER):
        if dtype == str:
            val = decode_string(d[idx:])
        else:
            if fmt:
                v = struct.unpack('>'+fmt, d[idx: idx+nbytes])[0]
                if isinstance(dtype, dict):
                    val = dtype.get(int(v), "N/A")
                else:
                    val = dtype(v)
            else:
                val = dtype(d[idx: idx+nbytes])
        header[nam] = val
        idx += nbytes
    resu = []
    for i in range(header["Number of elements"]):
        resu.append(decode_float(d[idx: idx+4]))
        idx += 4
    return header, numpy.array(resu, dtype=float)
    
def format_header(header, head_struct, columns=80):
    todisp = []
    for row in head_struct:
        key = row[0]        
        val = header.get(key, "N/A")
        if isinstance(val, basestring):
            val = repr(val)
        else:
            val = str(val)
        todisp.append((key+":", val))
    maxk = max([len(k) for k, v in todisp])
    maxv = max([len(v) for k, v in todisp])
    fmt = "%%-%ds %%-%ds"%(maxk, maxv)
    w = maxk+maxv+4
    ncols = columns/w
    nrows = len(todisp)/ncols
    print "w=", w
    print "ncols=", ncols
    print "nrows=", nrows
    res = ""
    for i in range(nrows):
        res += "| ".join([fmt%todisp[j*nrows+i] for j in range(ncols)]) + "\n"
    return res
          

if __name__ == "__main__":
    import sys
    import optparse
    opt = optparse.OptionParser("A simple tool for tracing a dumped trace")
    opt.add_option('-f', '--filename', default=None,
                   dest='filename',
                   help='Output filename. If not set, read from stdin')
    opt.add_option('-m', '--mode', default='binary',
                   dest='mode',
                   help='Dumping mode (may be "binary" [default], "ascii" or "ansi")',
                   )
    opt.add_option('-d', '--display-header', default=False,
                   action="store_true",
                   dest="displayheader",
                   help="Display the trace header")
    opt.add_option('-P', '--noplot-trace', default=True,
                   action="store_false",
                   dest="plot",
                   help="Do not display the plot of the trace")
    opt.add_option('-x', '--xmode', default='lin',
                   dest='xmode',
                   help='X coordinate mode (may be "lin" [default] or "log")')
    opt.add_option('-y', '--ymode', default='lin',
                   dest='ymode',
                   help='Y coordinate mode (may be "lin" [default], "log" or "db")')
    
    options, argv = opt.parse_args(sys.argv)


    if options.filename is None:
        print "Can't deal stdin for now..."
        sys.exit(1)
    try:
        header, data = decode_trace(open(options.filename, 'rb').read())
    except Exception, e:
        print "ERROR: can't read %s an interpret it as a HP3562 trace"%options.filename
        print e
        sys.exit(1)

    if options.displayheader:
        print format_header(header, HEADER, 100)
    if options.plot:
        f0 = header['Start freq value']
        dx = header['Delta X-axis']
        n = header['Number of elements']
        x = numpy.linspace(f0, f0+dx*n, len(data)) 
        y = data.copy()

        import pylab
        if options.ymode != "lin":
            minv = min(y[y>0])
            y[y==0] = minv
            y = numpy.log10(y)
        if options.ymode == "db":
            y = y*10
            pylab.ylabel('db')
        pylab.grid()
        pylab.plot(x, y)
        pylab.xlabel('frequency')
        pylab.show()