diff options
| -rwxr-xr-x | tools/gr-usrptest/apps/uhd_phase_alignment.py | 545 |
1 files changed, 545 insertions, 0 deletions
diff --git a/tools/gr-usrptest/apps/uhd_phase_alignment.py b/tools/gr-usrptest/apps/uhd_phase_alignment.py new file mode 100755 index 000000000..ff36c4e7c --- /dev/null +++ b/tools/gr-usrptest/apps/uhd_phase_alignment.py @@ -0,0 +1,545 @@ +#!/usr/bin/env python +# +# Copyright 2018 Ettus Research, a National Instruments Company +# +# SPDX-License-Identifier: GPL-3.0-or-later +# +""" +UHD Phase Alignment: Phase alignment test using the UHD Python API. +""" + + +import argparse +from builtins import input +from datetime import datetime, timedelta +import itertools as itt +import sys +import time +import logging +import numpy as np +import numpy.random as npr +import uhd + + +CLOCK_TIMEOUT = 1000 # 1000mS timeout for external clock locking +INIT_DELAY = 0.05 # 50mS initial delay before transmit +CMD_DELAY = 0.05 # set a 50mS delay in commands +NUM_RETRIES = 10 # Number of retries on a given trial before giving up +# TODO: Add support for TX phase alignment + + +def parse_args(): + """Parse the command line arguments""" + description = """UHD Phase Alignment (Python API) + + Currently only supports RX phase alignment + + Example usage: + - Setup: 2x X310's (one with dboard in slot A, one in slot B) + + uhd_phase_alignment.py --args addr0=ADDR0,addr1=ADDR1 --rate 5e6 --gain 30 + --start-freq 1e9 --stop-freq 2e9 --freq-bands 3 + --clock-source external --time-source external --sync pps + --subdev "A:0" "A:0" --runs 3 --duration 1.0 + + Note: when specifying --subdev, put each mboard's subdev in "" + """ + # TODO: Add gain steps! + parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter, + description=description) + # Standard device args + parser.add_argument("--args", default="", type=str, + help="UHD device address args (requires 2 MBoards)") + parser.add_argument("--rate", type=float, default=5e6, + help="specify to perform a rate test (sps)") + parser.add_argument("--gain", type=float, default=10., + help="specify a gain setting for the device") + parser.add_argument("--channels", default=[0, 1], nargs="+", type=int, + help="which channel(s) to use " + "(specify 0 1 or 0 1 2 3)") + parser.add_argument("--duration", default=0.25, type=float, + help="duration for each capture in seconds") + parser.add_argument("--runs", default=10, type=int, + help="Number of times to retune and measure phase alignment") + # Test configuration + parser.add_argument("--start-freq", type=float, required=True, + help="specify a minimum frequency") + parser.add_argument("--stop-freq", type=float, required=True, + help="specify a maximum frequency") + parser.add_argument("--freq-bands", type=float, required=True, + help="specify the number of frequency bands to test") + parser.add_argument("--start-power", type=float, default=-30, + help="specify a starting output power for the siggen (dBm)") + parser.add_argument("--power-step", type=float, default=0, + help="specify the increase in siggen output power at each step") + parser.add_argument("--tone-offset", type=float, default=1e6, + help="Frequency offset of the input signal (ie. the " + "difference between the device's center frequency " + "and the test tone)") + parser.add_argument("--drift-threshold", type=float, default=2., + help="Maximum frequency drift (deg) while testing a given frequency") + parser.add_argument("--stddev-threshold", type=float, default=2., + help="Maximum frequency deviation (deg) over a single receive call") + # Device configuration + parser.add_argument("--clock-source", type=str, + help="clock reference (internal, external, mimo, gpsdo)") + parser.add_argument("--time-source", type=str, + help="PPS source (internal, external, mimo, gpsdo)") + parser.add_argument("--sync", type=str, default="default", + #choices=["default", "pps", "mimo"], + help="Method to synchronize devices)") + parser.add_argument("--subdev", type=str, nargs="+", + help="Subdevice(s) of UHD device where appropriate. Use " + "a space-separated list to set different boards to " + "different specs.") + # Extra, advanced arguments + parser.add_argument("--plot", default=False, action="store_true", + help="Plot results") + parser.add_argument("--save", default=False, action="store_true", + help="Save each set of samples") + parser.add_argument("--easy-tune", type=bool, default=True, + help="Round the target frequency to the nearest MHz") + args = parser.parse_args() + + # Do some sanity checking + if args.tone_offset >= (args.rate / 2): + logger.warning("Tone offset may be outside the received bandwidth!") + + return args + + +class LogFormatter(logging.Formatter): + """Log formatter which prints the timestamp with fractional seconds""" + @staticmethod + def pp_now(): + """Returns a formatted string containing the time of day""" + now = datetime.now() + return "{:%H:%M}:{:05.2f}".format(now, now.second + now.microsecond / 1e6) + + def formatTime(self, record, datefmt=None): + converter = self.converter(record.created) + if datefmt: + formatted_date = converter.strftime(datefmt) + else: + formatted_date = LogFormatter.pp_now() + return formatted_date + + +def setup_ref(usrp, ref, num_mboards): + """Setup the reference clock""" + if ref == "mimo": + if num_mboards != 2: + logger.error("ref = \"mimo\" implies 2 motherboards; " + "your system has %d boards", num_mboards) + return False + usrp.set_clock_source("mimo", 1) + else: + usrp.set_clock_source(ref) + + # Lock onto clock signals for all mboards + if ref != "internal": + logger.debug("Now confirming lock on clock signals...") + end_time = datetime.now() + timedelta(milliseconds=CLOCK_TIMEOUT) + for i in range(num_mboards): + if ref == "mimo" and i == 0: + continue + is_locked = usrp.get_mboard_sensor("ref_locked", i) + while (not is_locked) and (datetime.now() < end_time): + time.sleep(1e-3) + is_locked = usrp.get_mboard_sensor("ref_locked", i) + if not is_locked: + logger.error("Unable to confirm clock signal locked on board %d", i) + return False + return True + + +def setup_pps(usrp, pps, num_mboards): + """Setup the PPS source""" + if pps == "mimo": + if num_mboards != 2: + logger.error("ref = \"mimo\" implies 2 motherboards; " + "your system has %d boards", num_mboards) + return False + # make mboard 1 a slave over the MIMO Cable + usrp.set_time_source("mimo", 1) + else: + usrp.set_time_source(pps) + return True + + +def setup_usrp(args): + """Create, configure, and return the device + + The USRP object that is returned will be synchronized and ready to receive. + """ + usrp = uhd.usrp.MultiUSRP(args.args) + + # Always select the subdevice first, the channel mapping affects the other settings + if args.subdev: + assert len(args.subdev) == usrp.get_num_mboards(),\ + "Please specify a subdevice spec for each mboard" + for mb_idx in range(usrp.get_num_mboards()): + usrp.set_rx_subdev_spec(uhd.usrp.SubdevSpec(args.subdev[mb_idx]), mb_idx) + + else: + logger.warning("No RX subdev specs set! Please ensure that the correct " + "connections are being used.") + + logger.info("Using Device: %s", usrp.get_pp_string()) + + # Set the reference clock + if args.clock_source and not setup_ref(usrp, args.clock_source, usrp.get_num_mboards()): + # If we wanted to set a reference clock and it failed, return + return None + + # Set the PPS source + if args.time_source and not setup_pps(usrp, args.time_source, usrp.get_num_mboards()): + # If we wanted to set a PPS source and it failed, return + return None + # At this point, we can assume our device has valid and locked clock and PPS + + # Determine channel settings + # TODO: Add support for >2 channels! (TwinRX) + if len(args.channels) != 2: + logger.error("Must select 2 channels! (%s selected)", args.channels) + return None + logger.info("Selected %s RX channels", args.channels if args.channels else "no") + # Set the sample rate + for chan in args.channels: + usrp.set_rx_rate(args.rate, chan) + + # Actually synchronize devices + # We already know we have >=2 channels, so don't worry about that + if args.sync in ['default', "pps"]: + logger.info("Setting device timestamp to 0...") + usrp.set_time_unknown_pps(uhd.types.TimeSpec(0.0)) + elif args.sync == 'mimo': + # For MIMO, we want to set the time on the master and let it propogate + # through the MIMO cable + usrp.set_time_now(uhd.types.TimeSpec(0.0), 0) + time.sleep(1) + logger.info("Current device timestamp: %.8f", + usrp.get_time_now().get_real_secs()) + else: + # This should never happen- argparse choices should handle this + logger.error("Invalid sync option for given configuration: %s", args.sync) + return None + + return usrp + + +def get_band_limits(start_freq, stop_freq, freq_bands): + """Return an array of length `freq_bands + 1`. + Each element marks the start of a frequency band (Hz). + Bands are equal sized (not log or anything fancy). + The last element is the stop frequency. + ex. get_band_limits(10., 100., 2) => [10., 55., 100.] + """ + return np.linspace(start_freq, stop_freq, freq_bands+1, endpoint=True) + + +def window(seq, width=2): + """Returns a sliding window (of `width` elements) over data from the iterable. + s -> (s0,s1,...s[n-1]), (s1,s2,...,sn), ... + Itertools example found at https://docs.python.org/release/2.3.5/lib/itertools-example.html + """ + seq_iter = iter(seq) + result = tuple(itt.islice(seq_iter, width)) + if len(result) == width: + yield result + for elem in seq_iter: + result = result[1:] + (elem,) + yield result + + +def generate_time_spec(usrp, time_delta=0.05): + """Return a TimeSpec for now + `time_delta`""" + return usrp.get_time_now() + uhd.types.TimeSpec(time_delta) + + +def tune_siggen(freq, power_lvl): + """Tune the signal generator to output the correct tone""" + # TODO: support actual RTS equipment, or any automated way + input("Please tune the signal generator to {:.3f} MHz and {:.1f} dBm, " + "then press Enter".format(freq / 1e6, power_lvl)) + + +def tune_usrp(usrp, freq, channels, delay=CMD_DELAY): + """Synchronously set the device's frequency""" + usrp.set_command_time(generate_time_spec(usrp, time_delta=delay)) + for chan in channels: + usrp.set_rx_freq(uhd.types.TuneRequest(freq), chan) + + +def recv_aligned_num_samps(usrp, streamer, num_samps, freq, channels=(0,)): + """ + RX a finite number of samples from the USRP + :param usrp: MultiUSRP object + :param streamer: RX streamer object + :param num_samps: number of samples to RX + :param freq: RX frequency (Hz) + :param channels: list of channels to RX on + :return: numpy array of complex floating-point samples (fc32) + """ + # Allocate a sample buffer + result = np.empty((len(channels), num_samps), dtype=np.complex64) + + # Tune to the desired frequency + tune_usrp(usrp, freq, channels) + + metadata = uhd.types.RXMetadata() + buffer_samps = streamer.get_max_num_samps() * 10 + recv_buffer = np.zeros( + (len(channels), buffer_samps), dtype=np.complex64) + recv_samps = 0 + + stream_cmd = uhd.types.StreamCMD(uhd.types.StreamMode.start_cont) + stream_cmd.stream_now = False + stream_cmd.time_spec = generate_time_spec(usrp) + stream_cmd.num_samps = num_samps + streamer.issue_stream_cmd(stream_cmd) + logger.debug("Sending stream command for T=%.2f", stream_cmd.time_spec.get_real_secs()) + + samps = np.array([], dtype=np.complex64) + while recv_samps < num_samps: + samps = streamer.recv(recv_buffer, metadata) + + if metadata.error_code != uhd.types.RXMetadataErrorCode.none: + # If we get a timeout, retry MAX_TIMEOUTS times + if metadata.error_code == uhd.types.RXMetadataErrorCode.timeout: + logger.error("%s (%d samps recv'd)", metadata.strerror(), recv_samps) + recv_samps = 0 + break + + real_samps = min(num_samps - recv_samps, samps) + result[:, recv_samps:recv_samps + real_samps] = recv_buffer[:, 0:real_samps] + recv_samps += real_samps + + logger.debug("Stopping stream") + stream_cmd = uhd.types.StreamCMD(uhd.types.StreamMode.stop_cont) + streamer.issue_stream_cmd(stream_cmd) + + logger.debug("Flushing stream") + # Flush the remainder of the samples + while samps: + samps = streamer.recv(recv_buffer, metadata) + + if recv_samps < num_samps: + logger.warning("Received too few samples, returning an empty array") + return np.array([], dtype=np.complex64) + return result + + +def plot_samps(samps, alignment): + """ + Show a nice plot of samples and their phase alignment + """ + try: + import pylab as plt + except ImportError: + logger.error("--plot requires pylab.") + return + + plt.tick_params(axis="both", labelsize=20) + # Plot the samples + plt.plot(samps[0][1000:2000].real, 'b') + plt.plot(samps[1][1000:2000].real, 'r') + plt.title("Phase Aligned RX", fontsize=44) + plt.legend(["Device A", "Device B"], fontsize=24) + plt.ylabel("Amplitude (real)", fontsize=35) + plt.xlabel("Time (us)", fontsize=35) + plt.show() + # Plot the alignment + logger.info("plotting alignment") + plt.plot(alignment) + plt.title("Phase Difference between Devices", fontsize=40) + plt.ylabel("Phase Delta (radian)", fontsize=30) + plt.xlabel("Time (us)", fontsize=30) + plt.ylim([-np.pi, np.pi]) + plt.show() + + +def check_results(alignment_stats, drift_thresh, stddev_thresh): + """Print the alignment stats in a nice way + + alignment_stats should be a dictionary of the following form: + {test_freq : [list of runs], ...} + ... the list of runs takes the form: + [{dictionary of run statistics}, ...] + ... the run dictionary has the following keys: + mean, stddev, min, max, test_freq, run_freq + ... whose values are all floats + """ + success = True # Whether or not we've exceeded a threshold + msg = "" + for freq, stats_list in alignment_stats.items(): + # Try to grab the test frequency for the frequency band + try: + test_freq = stats_list[0].get("test_freq") + except (KeyError, IndexError): + test_freq = 0. + logger.error("Failed to find test frequency for test band %.2fMHz", freq) + msg += "=== Frequency band starting at {:.2f}MHz. ===\n".format(freq/1e6) + msg += "Test Frequency: {:.2f}MHz ===\n".format(test_freq/1e6) + + # Allocate a list so we can calulate the drift over a set of runs + mean_list = [] + + for run_dict in stats_list: + run_freq = run_dict.get("run_freq", 0.) + # Convert mean and stddev to degrees + mean_deg = run_dict.get("mean", 0.) * 180 / np.pi + stddev_deg = run_dict.get("stddev", 0.) * 180 / np.pi + if stddev_deg > stddev_thresh: + success = False + + msg += "{:.2f}MHz<-{:.2f}MHz: {:.3f} deg +- {:.3f}\n".format( + test_freq/1e6, run_freq/1e6, mean_deg, stddev_deg + ) + mean_list.append(mean_deg) + + # Report the largest difference in mean values of runs + # FIXME: This won't work around +-180 deg + max_drift = max(mean_list) - min(mean_list) + if max_drift > drift_thresh: + success = False + msg += "--Maximum drift over runs: {:.2f} degrees\n".format(max_drift) + # Print a newline to separate frequency bands + msg += "\n" + + logger.info("Printing statistics!\n%s", msg) + return success + + +def main(): + """RX samples and write to file""" + args = parse_args() + + # Setup a usrp device + usrp = setup_usrp(args) + if usrp is None: + return False + + ### General test description ### + # 1. Split the frequency range of our device into bands. For each of these + # bands, we'll pick a random frequency within the band to be our test + # frequency. + # 2. Again split the frequency range of our device into bands, this time + # using the number of trials we want to run to split the range. Pick a + # random frequency within each run band. Tune to that run frequency, then + # back to our test frequency. + # 3. Receive synchronized samples, and determine the phase alignment. Report + # statistics based on the alignment. + # 4. Once we've iterated through each test frequency, determine whether or + # not the test passed or failed. + + # Determine the frequency bands we need to test + # TODO: allow users to specify test frequencies in args + freq_bands = get_band_limits(args.start_freq, args.stop_freq, args.freq_bands) + # Frequency bands to tune away to + # TODO: make this based on the device's frequency range. This requires + # additional Python API bindings. + run_bands = get_band_limits(args.start_freq, args.stop_freq, args.runs) + + nsamps = int(args.duration * args.rate) + st_args = uhd.usrp.StreamArgs("fc32", "sc16") + st_args.channels = args.channels + streamer = usrp.get_rx_stream(st_args) + + # Make a big dictionary to store all of the reported statistics + # Keys are the starting test frequency of the band + # Values are lists of dictionaries of statistics + all_alignment_stats = {} + # Test phase alignment in each test frequency band + current_power = args.start_power + for freq_start, freq_stop in window(freq_bands): + # Pick a random center frequency between the start and stop frequencies + tune_freq = npr.uniform(freq_start, freq_stop) + if args.easy_tune: + # Round to the nearest MHz + tune_freq = np.round(tune_freq, -6) + # Request the SigGen tune to our test frequency plus some offset away + # the device's LO + tune_siggen(tune_freq + args.tone_offset, current_power) + + # This is where the magic happens! + # Store phase alignment statistics as a list of dictionaries + alignment_stats = [] + for tune_away_start, tune_away_stop in window(run_bands): + # Try to get samples + for i in range(NUM_RETRIES): + # Tune to a random frequency in each of the frequency bands... + tune_away_freq = npr.uniform(tune_away_start, tune_away_stop) + tune_usrp(usrp, tune_away_freq, args.channels) + time.sleep(0.5) + + logger.info("Receiving samples, take %d, (%.2fMHz -> %.2fMHz)", + i, tune_away_freq/1e6, tune_freq/1e6) + + # Then tune back to our desired test frequency, and receive samples + samps = recv_aligned_num_samps(usrp, + streamer, + nsamps, + tune_freq, + args.channels) + if samps.size >= nsamps: + break + else: + streamer = None # Help the garbage collector + time.sleep(1) + streamer = usrp.get_rx_stream(st_args) + + # If we have failed to get good samples, put an empty dict in the stats + else: + logger.error("Failed to receive aligned samples!") + alignment_stats.append({}) + continue + + alignment = np.angle(np.conj(samps[0]) * samps[1])[500:] + + if args.plot: + plot_samps(samps, alignment,) + + if args.save: + # TODO: add frequency data + date_now = datetime.utcnow() + epoch = datetime(1970, 1, 1) + utc_now = int((date_now - epoch).total_seconds()) + np.savez("phaseAligned_{}.npz".format(utc_now), samps) + + # Store the phase alignment stats + alignment_stats.append({ + "mean": np.mean(alignment), + # Subtract the mean before calculating the stddev so we don't + # have rollover errors + "stddev": np.std(alignment - np.mean(alignment)), + "min": alignment.min(), + "max": alignment.max(), + "test_freq": tune_freq, + "run_freq": tune_away_freq + }) + run_means = [run_stats.get("mean", 0.) for run_stats in alignment_stats] + run_stddevs = [run_stats.get("stddev", 0.) for run_stats in alignment_stats] + logger.debug("Test freq %.3fMHz health check: %.1f deg drift, %.2f deg max stddev", + tune_freq/1e6, + max(run_means) - min(run_means), # FIXME: This won't work around +-180 deg + max(run_stddevs) + ) + all_alignment_stats[freq_start] = alignment_stats + # Increment the power level for the next run + current_power += args.power_step + + return check_results(all_alignment_stats, args.drift_threshold, args.stddev_threshold) + + +if __name__ == "__main__": + # Setup the logger with our custom timestamp formatting + global logger + logger = logging.getLogger(__name__) + logger.setLevel(logging.DEBUG) + console = logging.StreamHandler() + logger.addHandler(console) + formatter = LogFormatter(fmt='[%(asctime)s] [%(levelname)s] %(message)s') + console.setFormatter(formatter) + + sys.exit(not main()) |