read_trace.py

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

import struct


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_file(filename):
    d = open(filename).read()

    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])

        print idx, hex(idx), nam, ":", val
        header[nam] = val
        idx += nbytes

    resu = []

    for i in range(header["Number of elements"]):
        resu.append(decode_float(d[idx: idx+4]))
        idx += 4

    #print "resu = ", resu
    #return
    import pylab
    import numpy
    resu = numpy.array(resu, dtype=float)
    print "max = ", max(resu)
    #xr = numpy.linspace(0, header['Delta X-axis'], len(resu))
    sf = header['Start freq value']
    xr = numpy.linspace(sf, sf+20000, len(resu))
    mn = min(resu[resu>0])
    resu[resu==0] = mn
    pylab.plot(xr, 10*numpy.log10(resu))
    pylab.show()

    # tt=0
    # for i, (nam, dtype, nbytes) in enumerate(HEADER):
    #     if dtype == str:
    #       nb = ord(struct.unpack('c', d[idx])[0])
    #       val = d[idx+1:idx+1+nb]
    #     else:
    #         v = struct.unpack('>d', d[idx: idx+(nbytes*4)])[0]
    #         if isinstance(dtype, dict):
    #             val = dtype.get(int(v), "N/A")
    #         else:
    #             val = dtype(v)
    #     print idx, nam, ":", val
    #     idx += nbytes*4
        
                             
if __name__ == "__main__":
    import sys
    decode_file(sys.argv[1])