receive_USRP_data.m

% MATLAB script that sends (transmits) the simulated radar plus noise waveforms
% Step 1: load a single simulated waveform Matlab workspace generated by the NIST CBRS band radar waveform generator (see link below)
% Step 2: Configure the USRP radio to receive waveform sample data
% Step 3: Enter loop to receive a finite number of large USRP IQ sample buffers and store them in an array
% Step 4: Find and extract the first complete transmitted dataset in the received waveform vector
% Step 5: Plot magnitude of time domain IQ samples for basic data quality checking
%
% NIST CBRS radar waveform generator sourcecode: https://github.com/usnistgov/SimulatedRadarWaveformGenerator
% This script was tested with MATLAB 2020a on a home build PC running Ubuntu Linux 18.04 (i7-920 Quad Core, 2.67 GHz, 24GB RAM, 1TB HDD)
% Tested with a USRP-N210 with Daughter Board: 'SBXv3 RX'
% Note: because of the relatively high processing requirements for this script, I've found that it works best by 
% first closing all other MATLAB plot windows before running the script from the start. Otherwise, I get overflows that corrupt the data samples.
% Also, disable any other applications or background processes that could reduce MATLAB performance.
%
% Author Info: Alex Lackpour, alackpour@gmail.com
% Date: April 5th, 2021
% Version: 1.0, public release

clear; clc

receivingUSRPIPAddress = '192.168.11.2';

%% Step 1: load a single simulated waveform Matlab workspace generated by the NIST CBRS band radar waveform generator
% waveform signal labels are taken from this workspace file

fprintf('Step 1: Loading original (source) Matlab workspace that is being sent to this receiving node.\n');

originalDataPath = '../Group6_data/Group6_orig/Group6/';        % Directory path to the original source data that is being sent
experimentalDataPath = '../Group6_data/Group6_exp/Group6/';     % Directory path to the experimental data that is being received and stored for training and testing
cd(originalDataPath)

% User needs to select a particular MATLAB workspace containing the waveform that is being sent
[fileName, dataRootFolder] = uigetfile('group6*.mat');

load(fileName)

%% Step 2: Configure the USRP radio to receive data
fprintf('Step 2: Configuring USRP to receive data\n');
% USRP N200/N210/USRP2 master clock rate is 100 MHz
USRPclockRate_MHz = 100e6;
actualSampleRate_MHz = 10e6; % NOTE: on my Titan Ubuntu 18 Workstation, 20MSps had many overflows, 12.5 MSps worked for a few seconds until it overflowed, I don't get any overslows out to 20 seconds at 10 MSps
finalSampleRate_MHz = 10e6;

decimationFactor = USRPclockRate_MHz / actualSampleRate_MHz;
Fs_Hz = USRPclockRate_MHz/decimationFactor;

rxSDR = comm.SDRuReceiver('IPAddress', receivingUSRPIPAddress, 'DecimationFactor', decimationFactor);

% Experimentally found that requesting the largest buffer size that MATLAB USRP library supports is best for receiving data without overflows
rxSDR.SamplesPerFrame = 375000;

% Use double precision for the highest sample dynamic range for waveform accuracy
% Also tried using single precision for sample transport but did not appear to reduce occurance of overflow
rxSDR.TransportDataType = 'int16';         
rxSDR.OutputDataType = 'double';

% Notes on characteristics of USRP Rx Daughter Board: 'SBXv3 RX'
% Minimum Frequency: 380 MHz, Maximum Frequency, 4420 MHz
% Minimum Gain: 0, Maximum Gain: 38

rxSDR.Gain = 10;
rxSDR.CenterFrequency = 1500e6;
% NOTE: It was important to include a LO offset on the transmitter,
% setting an LO offset on the receiver didn't seem to have much effect on the floor of the signal during the start sync vector
rxSDR.LocalOscillatorOffset = 10e6;

%% Step 3: Enter loop to receive a finite number of large USRP IQ sample buffers and store them in an array
fprintf('Step 3: Entering USRP receive loop\n');
% Note: The duration of each waveform signal and number of signals being sent in
% a single dataset batch file needs to match the configuration at the sending node
SigDuration_sec = 0.08;
NSigs = 200;

NSamplesPerSig = SigDuration_sec*finalSampleRate_MHz;

fprintf('Releasing the Receiver Object\n');
release(rxSDR);

% Note: On the PC that this script was developed, a 26 second receive capture often correctly captured a 16 second waveform dataset (200 radar waveform samples, 80ms per saveform)
% This amount of time did not lead to any receive overflows, which is critical for obtaining good time domain samples without dropouts (missing sections)
% However, if I had a PC with more memory and processing speed, I probably would have increased the capture time to 32 seconds and guaranteed that the receive capture buffer contained the full 16 seconds of transmitted data
totalRxTime_seconds = 26;    

TotalRxPackets = round(totalRxTime_seconds*Fs_Hz/rxSDR.SamplesPerFrame);

OVER = zeros(1, TotalRxPackets);
savedData = zeros(rxSDR.SamplesPerFrame, TotalRxPackets);

skipCountDisplay = 50;
skipCounts = 0;
fprintf('Starting the reception of samples\n');
for ind = 1:TotalRxPackets
    
    [savedData(:, ind), ~, OVER(ind)] = step(rxSDR); % Receive buffer of IQ samples from USRP radio
    
    if mod(ind, skipCountDisplay) == 0
        bufInds = skipCounts*skipCountDisplay + (1:skipCountDisplay);
        bufNumOverFlows = sum(OVER(bufInds));
        fprintf('Loop iteration %d of %d. Number of receiver overflows in this loop iteration = %d\n', ind, TotalRxPackets, bufNumOverFlows);
        skipCounts = skipCounts + 1;
    end
    
    % Experimentally found that waiting 6/10ths of a fraction of the amount of time needed to read
    % in 375k samples avoids overflows and return-to-1 sprurious sample reads. If I decreased the pause time, 
    % then there would be many IQ data samples with a value of 1. If I increased the pause time, then many overflows would occur
    pause(rxSDR.SamplesPerFrame*6/10/Fs_Hz)
    
end

release(rxSDR)

% Convert the matrix of received data into a data vector for post-processing
totalPlotPackets = size(savedData,2);
dataComplexVector = reshape(savedData, 1, prod(size(savedData)));

clear savedData

if actualSampleRate_MHz > finalSampleRate_MHz
    dataComplexVector = resample(dataComplexVector, finalSampleRate_MHz*2, actualSampleRate_MHz*2);
end

dataMagLogVector = single(20*log10(abs(dataComplexVector)));    % convert to single precision FLOAT to reduce RAM
sampInds = single(1:numel(dataMagLogVector));                   % convert to single precision INT to reduce RAM

%% Step 4: Find and extract the first complete transmitted dataset in the received waveform vector
fprintf('Step 4: Searching for a complete batch of the sent radar waveforms');
threshold_dB = -40;
aboveInds = single(find(dataMagLogVector > threshold_dB));
belowInds = single(find(dataMagLogVector < threshold_dB));
lengthAboveInds = single(diff(aboveInds));

% Following operation attempts to detect the presence of an empty signal slot sent at the beginning of each vector of radar signals
% This is done as a post processing synchronization of the sent waveform labels and the asynchronously received data samples
syncStartThresholdedInds = find(lengthAboveInds > NSamplesPerSig*9/10);
syncStartInds = sampInds(aboveInds(syncStartThresholdedInds));
syncStartInds = double(syncStartInds);

multipleSyncStartIndsDetected = numel(syncStartInds) > 1;

if multipleSyncStartIndsDetected
    
    startInd = syncStartInds(1) + NSamplesPerSig;
    
    % Check that detected radar waveform is approximately the correct length between the two sync vectors
    sampleLengthErrorThreshold = 100; % Reasonable max tolerable difference between the measured and expected length of sequentially transmitted and received radar waveform samples
    measuredLengthInds = syncStartInds(2) - startInd;
    expectedLengthInds = NSigs*NSamplesPerSig;
    measuredDifferenceInds = measuredLengthInds - expectedLengthInds;
    goodWaveformLengthTest = (abs(measuredDifferenceInds) <= sampleLengthErrorThreshold);
    
    if goodWaveformLengthTest
        
        stopInd = startInd + expectedLengthInds - 1;
        
        receivedDataMatrix = reshape(dataComplexVector(startInd:stopInd), NSamplesPerSig, NSigs);
        fprintf('\nGood News: Successfully found two sync vectors in the received data! Difference between measured and expected numbers of samples is %d and that is below threshold %d\n', abs(measuredDifferenceInds), sampleLengthErrorThreshold)
        
        fileNameWOExt = strsplit(fileName, '.');
        fileNameWOExt = fileNameWOExt{1};
        
        finalFileNameUnderScoreInd = max(strfind(fileNameWOExt, '_'))+1;
        
        fileIDString = fileNameWOExt(finalFileNameUnderScoreInd:end);
        waveformSubsetName = ['waveformSubset_' fileIDString];
        radarStateusSubsetName = ['radarStatusSubset_' fileIDString];
        waveformTableSubsetName = ['waveformTableSubset_' fileIDString];
        
        assignin('base', ['waveformSubset_' fileNameWOExt(finalFileNameUnderScoreInd:end)], receivedDataMatrix)
        
        fprintf('Saving the experimental data in the updated MATLAB workspace\n');
        % Move up a couple of directories because we are inside the original dataset directory
        save(['../../' experimentalDataPath fileNameWOExt '.mat'], 'FInfo', radarStateusSubsetName, waveformSubsetName, waveformTableSubsetName')
        
    else
        fprintf('\nPost Processing Error: Found two sync vectors (zeros) but they are not the expected distance appart (%d > %d).\n', abs(measuredDifferenceInds), sampleLengthErrorThreshold);
        fprintf('Please run this script again to get a better capture and/or change its configuration parameters to reduce overflows!\n');
        startInd = syncStartInds(1) + NSamplesPerSig;
        stopInd = syncStartInds(2);
    end
else % Case of what occurs when 
    fprintf('\nPost Processing Error: Cannot find two sync vectors (zeros). Need to run this script again or change its parameters!\n');
end

% Free up some memory prior to plotting the results
clear aboveInds belowInds lengthAboveInds receivedDataMatrix

%% Step 5: Plot magnitude of time domain IQ samples for basic data quality checking

% Flags for enabling/disabling plotting
TIME_DOMAIN_PLOT_FLAG = 1;  % Enable/Disable time domain plotting
PLOT_ALL_DATA_FLAG = 0;     % 1 = plot ALL received IQ time samples, 0 = plot received IQ samples after being detected
SPECTROGRAM_PLOT_FLAG = 0;  % 1 = Plot spectrogram of the detected radar waveform (if radar data not detected, then plot spurious data)

figure(1); clf
stem(OVER); grid on
xlabel('Receive Buffer Index #')
ylabel('Status of Receive Overflow Flag');
title('Status of USRP Receive Overflow Flag per Rx Buffer Index');

if TIME_DOMAIN_PLOT_FLAG % Enable or Disable the Time domain vector plot
    figure(2); clf;
    
    if multipleSyncStartIndsDetected
        fprintf('\nStarting time domain plotting\n');

        if PLOT_ALL_DATA_FLAG   % plot ALL received IQ time samples
            plot(sampInds, dataMagLogVector, 'b', [startInd stopInd], [0 0], 'r*'); grid on
            xlabel('IQ Sample index #')
            ylabel('Unscaled Power IQ samples (dB)');
            title('Plot of All Received IQ Sample Magnitudes')
            
        else    % plot only the data near the start and stop detection inds - showing the zero-valued sync vector
            plotIndRange = 1000;
            dplotIndsRising = [(-plotIndRange:plotIndRange)+startInd];
            dplotIndsFalling = [(-plotIndRange:plotIndRange)+stopInd];
            
            subplot(1,2,1)
            plot(dataMagLogVector(dplotIndsRising), 'b.-'); grid on
            xlabel('IQ Sample index # on beginning of radar waveform signals')
            ylabel('Unscaled Power IQ samples (dB)');
            title('Plot of Beginning of Detected Buffer of Radar waveforms')
            
            subplot(1,2,2)
            plot(dataMagLogVector(dplotIndsFalling), 'b.-'); grid on
            xlabel('IQ Sample index # on ending of radar waveform signals')
            ylabel('Unscaled Power IQ samples (dB)');
            title('Plot of Ending of Detected Buffer of Radar waveforms')
        end
    end
end

if SPECTROGRAM_PLOT_FLAG & multipleSyncStartIndsDetected
    fprintf('\nStarting plotting of spectrogram\n');
    
    figure(3); clf
    
    NFFT = 16384;
    NOLAP = 16;
    F_Hz = linspace(-finalSampleRate_MHz/2, finalSampleRate_MHz/2-1/finalSampleRate_MHz, NFFT);
    winSpecFunc = chebwin(NFFT, 120); %  hamming(NFFT); %
    spectrogram(dataComplexVector(startInd:stopInd), winSpecFunc, NOLAP, F_Hz, finalSampleRate_MHz)
    crange = caxis;
    caxis([crange(2)-60 crange(2)])
    title('Spectrogram of Detected Vector of Sent Radar Signals')
end

drawnow
fprintf('Finished running receiver script!\n');