creotech-sayma-testsuite/plot_sayma_data.py

import numpy as np
import matplotlib.pyplot as plot
from scipy import signal, constants
import argparse
import os
import datetime


def rp_raw_to_numpy(rp_raw):
    # Convert raw buffer strings to numpy arrays
    buff_string = rp_raw.split(',')
    return np.array(list(map(float, buff_string)))


RP_IP_ADDRS = {
    "creotech-1": "192.168.1.104",
    "creotech-2": "192.168.1.105",
    "mlabs-1": "rp-f05cc9",
    "mlabs-2": "rp-f0612e",
}


def main():
    parser = argparse.ArgumentParser(description="Data plotting tool for Sayma DAC/TTL")
    parser.add_argument("dir", help="output directory", type=str)
    parser.add_argument("rps", metavar="NAME:CHANNEL",
        help="input <CHANNEL> of the RedPitaya at <NAME> where data is collected; "
        "CHANNEL must be 1 or 2; "
        "NAME must be any of: " + " ".join(list(RP_IP_ADDRS.keys())),
        type=str, nargs=2)
    parser.add_argument("--log",
        help="path of the log file where the measurement record will be appended",
        type=str)
    parser.add_argument("--noplot",
        help="do not show data plot, which is blocking until the GUI is closed",
        action="store_true")
    args = parser.parse_args()

    # Must only compare 2 data
    name, channel = args.rps[0].split(':')
    y1_filename = 'rp_{}_y{}_raw.bin'.format(name, channel)
    name, channel = args.rps[1].split(':')
    y2_filename = 'rp_{}_y{}_raw.bin'.format(name, channel)
    if y1_filename == y2_filename:
        raise ValueError("Both files are the same.")

    with open(os.path.join(args.dir, y1_filename), 'rb') as f:
        y1_now_iso, y1_raw = [l.decode('utf-8') for l in f.readlines()]
    with open(os.path.join(args.dir, y2_filename), 'rb') as f:
        y2_now_iso, y2_raw = [l.decode('utf-8') for l in f.readlines()]
    if None in [y1_raw, y2_raw]:
        raise IOError("Raw RP string files cannot be opened.")
    if y1_now_iso != y2_now_iso:
        raise ValueError("Timestamps of raw RP files are not identical.")
    now_iso = y1_now_iso.rstrip()
    print("Reading raw RP data collected at {}.".format(now_iso))

    y1 = rp_raw_to_numpy(y1_raw)
    y2 = rp_raw_to_numpy(y2_raw)

    # Define t as an array of timestamps (in seconds) for each sample
    # (Note that RedPitaya's oscilloscope has a sampling rate @ 125MHz)
    t = np.arange(y1.shape[0])/125e6
    # Generate matrix Y by having arrays y1 and y2 as 2 rows
    y = np.c_[y1, y2].T
    # Element-wise multiply an array of cos(2pi*9e6*t) with each row in Y;
    # Then, downsample the array by 10 as Z
    z = signal.decimate(y*np.exp(1j*2*np.pi*9e6*t), q=10, ftype="fir", zero_phase=True)[:, 10:]
    # Downsample Z by 10 again.
    z = signal.decimate(z, q=10, ftype="fir", zero_phase=True)[:, 10:]
    # Element-wise multiply Z[0] with the conjugate of Z[1] to get the phase difference (i.e. angle(z0) - angle(z1)), and use the mean value.
    angle = np.angle(np.mean(z[0]*z[1].conj()))

    # Append the phase difference to the log file
    log = args.log
    if log is not None:
        with open(log, 'a') as f:
            f.write("{}\t{}\n".format(now_iso, angle))
            print("Phase measurement record appended to log: {}".format(log))
    # Print the phase difference
    print(angle)

    if not args.noplot:
        # Normalize y1 and y2 for plotting
        y1 /= abs(y1).max()
        y2 /= abs(y2).max()
        plot.plot(y1)
        plot.plot(y2)
        plot.show()

if __name__ == "__main__":
    main()