diff --git a/Matlab/AnalogAmplitudeCalibrationFormat.m b/Matlab/AnalogAmplitudeCalibrationFormat.m
new file mode 100644
index 0000000..d402bb0
--- /dev/null
+++ b/Matlab/AnalogAmplitudeCalibrationFormat.m
@@ -0,0 +1,53 @@
+classdef AnalogAmplitudeCalibrationFormat
+    
+    properties
+        CourseSubarray1Data
+        CourseSubarray2Data
+        FineSubarray1Data
+        FineSubarray2Data
+    end
+    
+    methods
+        function s = toStruct(this)
+            s.CourseSubarray1Data = this.CourseSubarray1Data;
+            s.CourseSubarray2Data = this.CourseSubarray2Data;
+            s.FineSubarray1Data = this.FineSubarray1Data;
+            s.FineSubarray2Data = this.FineSubarray2Data;
+        end
+
+        function [cal1,cal2] = getCourseCalValues(this)
+            cal1 = this.getCourseCal1();
+            cal2 = this.getCourseCal2();
+        end
+
+        function [cal1,cal2] = getFineCalValues(this)
+            cal1 = this.getFineCal1();
+            cal2 = this.getFineCal2();
+        end
+    end
+
+    methods (Access = private)
+        function cal1 = getCourseCal1(this)
+            cal1 = this.getCal(this.CourseSubarray1Data);
+        end
+
+        function cal1 = getFineCal1(this)
+            cal1 = this.getCal(this.FineSubarray1Data);
+        end
+
+        function cal2 = getCourseCal2(this)
+            cal2 = this.getCal(this.CourseSubarray2Data);
+        end
+
+        function cal2 = getFineCal2(this)
+            cal2 = this.getCal(this.FineSubarray2Data);
+        end
+
+        function cal = getCal(~,data)
+            % get amplitudes
+            amp = mean(max(abs(fft(data))));
+            scaling = min(amp) ./ amp;
+            cal = reshape(scaling,[1 4]);
+        end
+    end
+end
\ No newline at end of file
diff --git a/Matlab/AnalogPhaseCalibrationFormat.m b/Matlab/AnalogPhaseCalibrationFormat.m
new file mode 100644
index 0000000..c22f5cd
--- /dev/null
+++ b/Matlab/AnalogPhaseCalibrationFormat.m
@@ -0,0 +1,45 @@
+classdef AnalogPhaseCalibrationFormat
+    
+    properties
+        PhaseSetting
+        Subarray1Measurements
+        Subarray2Measurements
+    end
+    
+    methods
+        function s = toStruct(this)
+            s.PhaseSetting = this.PhaseSetting;
+            s.Subarray1Measurements = this.Subarray1Measurements;
+            s.Subarray2Measurements = this.Subarray2Measurements;
+        end
+
+        function [cal1,cal2] = getCalValues(this)
+            cal1 = this.sub1Cal();
+            cal2 = this.sub2Cal();
+        end
+    end
+
+    methods (Access = private)
+        function cal1 = sub1Cal(this)
+            cal1 = this.subCal(this.Subarray1Measurements);
+        end
+
+        function cal2 = sub2Cal(this)
+            cal2 = this.subCal(this.Subarray2Measurements);
+        end
+
+        function cal = subCal(this,data)
+            % Get signal amplitude
+            datafft = abs(fft(data));
+            amp = max(datafft,[],1);
+
+            % Get location of min amplitude
+            [~,phaseIdx] = min(amp,[],2);
+            
+            % Get calibration phase
+            calPhase = wrapTo360(this.PhaseSetting(phaseIdx)-180);
+            cal = [0;calPhase.'];
+        end
+    end
+end
+
diff --git a/Matlab/AntennaInteractor.m b/Matlab/AntennaInteractor.m
new file mode 100644
index 0000000..59d445a
--- /dev/null
+++ b/Matlab/AntennaInteractor.m
@@ -0,0 +1,167 @@
+classdef AntennaInteractor < handle
+
+    properties
+        ArrayControl
+        PlutoControl
+        Model
+        Fc
+        SubSteer
+        ArraySteer
+        LastPhase
+        LastGain
+        AnalogWeights
+        DigitalWeights
+    end
+
+    methods
+        function this = AntennaInteractor(fc,calValues)
+            [rx,bf,model] = setupAntenna(fc);
+            this.ArrayControl = bf;
+            this.PlutoControl = rx;
+            this.Model = model;
+            this.Fc = fc;
+            this.SubSteer = phased.SteeringVector("SensorArray",this.Model.Subarray,'NumPhaseShifterBits',7);
+            this.ArraySteer = phased.SteeringVector("SensorArray",this.Model);
+            this.AnalogWeights = calValues.AnalogWeights;
+            this.DigitalWeights = calValues.DigitalWeights;
+        end
+
+        function updateCalibration(this,calValues)
+            this.updateAnalogWeights(calValues.AnalogWeights);
+            this.updateDigitalWeights(calValues.DigitalWeights);
+        end
+
+        function updateAnalogWeights(this,analogWeights)
+            this.AnalogWeights = analogWeights;
+        end
+
+        function updateDigitalWeights(this,digitalWeights)
+            this.DigitalWeights = digitalWeights;
+        end
+
+        function [patternData,rxdata] = capturePattern(this,steerangles)
+            patternData = zeros(1024,numel(steerangles));
+            for ii = 1 : numel(steerangles)
+                [analogweights,digitalweights] = this.getAllWeights(steerangles(ii));
+                rxdata = this.steerAnalog(analogweights);
+                patternData(:,ii) = rxdata * conj(digitalweights);
+            end
+        end
+
+        function [sumdiffampdelta,sumdiffphasedelta,sumpatterndata,diffpatternData] = captureMonopulsePattern(this,steerangles)
+            sumpatterndata = zeros(1024,numel(steerangles));
+            diffpatternData = zeros(1024,numel(steerangles));
+            for ii = 1 : numel(steerangles)
+                % capture data
+                [analogweights,digitalweights] = this.getAllWeights(steerangles(ii));
+                rxdata = this.steerAnalog(analogweights);
+
+                % sum
+                sumpatterndata(:,ii) = rxdata * conj(digitalweights);
+
+                % diff
+                diffdigitalweights = digitalweights .* [1;-1];
+                diffpatternData(:,ii) = rxdata * conj(diffdigitalweights);
+            end
+
+            % calulate sum and diff amplitude and phase deltas
+            sumamp = mag2db(helperGetAmplitude(sumpatterndata));
+            diffamp = mag2db(helperGetAmplitude(diffpatternData));
+            sumdiffampdelta = sumamp-diffamp;
+
+            sumphase = helperGetPhase(sumpatterndata);
+            diffphase = helperGetPhase(diffpatternData);
+            sumdiffphasedelta = sign(wrapTo180(sumphase-diffphase));
+        end
+
+        function patternData = capturePatternWithNull(this,steerangles,nullangle)
+            patternData = zeros(1024,numel(steerangles));
+            for ii = 1 : numel(steerangles)
+                [analogweights,digitalweights] = this.getAllWeightsNull(steerangles(ii),nullangle);
+                rxdata = this.steerAnalog(analogweights);
+                patternData(:,ii) = rxdata * conj(digitalweights);
+            end
+        end
+
+        function rxdata = steerAnalog(this,analogWeights)
+            % Set analog phase shifter
+            phases = this.getPhaseCodes(analogWeights);
+            if ~isequal(this.LastPhase,phases)
+                this.ArrayControl.RxPhase(:) = phases;
+                this.LastPhase = phases;
+            end
+        
+            % Set analog gain
+            gainCode = this.getGainCodes(analogWeights);
+            if ~isequal(this.LastGain,gainCode)
+                this.ArrayControl.RxGain(:) = gainCode;
+                this.LastGain = gainCode;
+            end
+        
+            % receive data
+            this.ArrayControl.LatchRxSettings();
+            this.PlutoControl();
+            rxdata = this.PlutoControl();
+        end
+
+        function phases = getPhaseCodes(this,analogWeights)
+            sub1weights = analogWeights(:,1);
+            sub2weights = analogWeights(:,2);
+            sub1phase = this.getPhase(sub1weights);
+            sub2phase = this.getPhase(sub2weights);
+            phases = [sub1phase',sub2phase'];
+            phases = phases - phases(1);
+            phases = wrapTo360(phases);
+        end
+
+        function codes = getGainCodes(~,analogWeights)
+            codes = helperGainCodes(analogWeights);
+        end
+    end
+
+    methods (Access = private)
+        function [analogweights,digitalweights] = getAllWeights(this,steerangle)
+            defaultAnalogWeights = this.AnalogWeights;
+            defaultDigitalWeights = this.DigitalWeights;
+
+            % get analog weights
+            analogweights = this.getWeights(steerangle,defaultAnalogWeights,this.SubSteer);
+
+            % get digital weights, flip to ensure they are being applied to
+            % the correct channel.
+            flippedDigitalWeights = [defaultDigitalWeights(2);defaultDigitalWeights(1)];
+            digitalweights = this.getWeights(steerangle,flippedDigitalWeights,this.ArraySteer);
+            digitalweights = [digitalweights(2);digitalweights(1)];
+        end
+
+        function [analogweights,digitalweights] = getAllWeightsNull(this,steerangle,nullangle)
+            defaultAnalogWeights = this.AnalogWeights;
+            defaultDigitalWeights = this.DigitalWeights;
+            
+            % get analogweights
+            analogweights = this.getWeightsNull(steerangle,nullangle,defaultAnalogWeights,this.SubSteer);
+
+            % get digital weights
+            flippedDigitalWeights = [defaultDigitalWeights(2);defaultDigitalWeights(1)];
+            digitalweights = this.getWeightsNull(steerangle,nullangle,flippedDigitalWeights,this.ArraySteer);
+            digitalweights = [digitalweights(2);digitalweights(1)];
+        end
+
+        function weights = getWeights(this,steerangle,defaultweights,sv)
+            initialweights = sv(this.Fc,steerangle);
+            weights = initialweights .* defaultweights;
+        end
+
+        function weights = getWeightsNull(this,steerangle,nullangle,defaultweights,sv)
+            steerweight = sv(this.Fc,steerangle);
+            nullweight = sv(this.Fc,nullangle);
+            rn = nullweight'*steerweight/(nullweight'*nullweight);
+            steernullweight = steerweight-nullweight*rn;
+            weights = steernullweight .* defaultweights;
+        end
+
+        function phase = getPhase(~,weights)
+            phase = rad2deg(angle(weights));
+        end
+    end
+end
\ No newline at end of file
diff --git a/Matlab/ArrayFactorDataFormat.m b/Matlab/ArrayFactorDataFormat.m
new file mode 100644
index 0000000..49380f9
--- /dev/null
+++ b/Matlab/ArrayFactorDataFormat.m
@@ -0,0 +1,27 @@
+classdef ArrayFactorDataFormat
+    
+    % Carry array factor data for calibration example
+    
+    properties
+        SteeringAngle
+        UncalibratedPattern
+        AnalogCourseAmplitudeCalPattern
+        AnalogPhaseCalPattern
+        AnalogFineAmplitudeCalPattern
+        DigitalAmplitudeCalPattern
+        FullCalibration
+    end
+
+    methods
+        function s = toStruct(this)
+            s.SteeringAngle = this.SteeringAngle;
+            s.UncalibratedPattern = this.UncalibratedPattern;
+            s.AnalogCourseAmplitudeCalPattern = this.AnalogCourseAmplitudeCalPattern;
+            s.AnalogPhaseCalPattern = this.AnalogPhaseCalPattern;
+            s.AnalogFineAmplitudeCalPattern = this.AnalogFineAmplitudeCalPattern;
+            s.DigitalAmplitudeCalPattern = this.DigitalAmplitudeCalPattern;
+            s.FullCalibration = this.FullCalibration;
+        end
+    end
+end
+
diff --git a/Matlab/CalibrationDataFormat.m b/Matlab/CalibrationDataFormat.m
new file mode 100644
index 0000000..097ae56
--- /dev/null
+++ b/Matlab/CalibrationDataFormat.m
@@ -0,0 +1,25 @@
+classdef CalibrationDataFormat
+    
+    % Store calibration data
+    
+    properties
+        ExampleData
+        CalibrationWeights = CalibrationWeightFormat();
+        AntennaPattern = ArrayFactorDataFormat();
+        AnalogPhaseCalibration = AnalogPhaseCalibrationFormat();
+        AnalogAmplitudeCalibration = AnalogAmplitudeCalibrationFormat();
+        DigitalCalibration = DigitalCalibrationFormat();
+    end
+
+    methods
+        function s = toStruct(this)
+            s.ExampleData = this.ExampleData;
+            s.CalibrationValues = this.CalibrationWeights.toStruct();
+            s.AntennaPattern = this.AntennaPattern.toStruct();
+            s.AnalogPhaseCalibration = this.AnalogPhaseCalibration.toStruct();
+            s.AnalogAmplitudeCalibration = this.AnalogAmplitudeCalibration.toStruct();
+            s.DigitalCalibration = this.DigitalCalibration.toStruct();
+        end
+    end
+end
+
diff --git a/Matlab/CalibrationValueFormat.m b/Matlab/CalibrationValueFormat.m
new file mode 100644
index 0000000..2b696af
--- /dev/null
+++ b/Matlab/CalibrationValueFormat.m
@@ -0,0 +1,17 @@
+classdef CalibrationValueFormat
+    
+    % Store the antenna calibration values
+    
+    properties
+        AnalogWeights
+        DigitalWeights
+    end
+    
+    methods
+        function s = toStruct(this)
+            s.AnalogWeights = this.AnalogWeights;
+            s.DigitalWeights = this.DigitalWeights;
+        end
+    end
+end
+
diff --git a/Matlab/CalibrationWeightFormat.m b/Matlab/CalibrationWeightFormat.m
new file mode 100644
index 0000000..5a9d133
--- /dev/null
+++ b/Matlab/CalibrationWeightFormat.m
@@ -0,0 +1,25 @@
+classdef CalibrationWeightFormat
+    
+    % Hold calibration weights at each step of the calibration process.
+    
+    properties
+        UncalibratedWeights = CalibrationValueFormat()
+        AnalogCourseAmplitudeWeights = CalibrationValueFormat()
+        AnalogPhaseWeights = CalibrationValueFormat()
+        AnalogFineAmplitudeWeights = CalibrationValueFormat()
+        DigitalAmplitudeWeights = CalibrationValueFormat()
+        FinalCalibrationWeights = CalibrationValueFormat()
+    end
+    
+    methods
+        function s = toStruct(this)
+            s.UncalibratedWeights = this.UncalibratedWeights.toStruct();
+            s.AnalogCourseAmplitudeWeights = this.AnalogCourseAmplitudeWeights.toStruct();
+            s.AnalogPhaseWeights = this.AnalogPhaseWeights.toStruct();
+            s.AnalogFineAmplitudeWeights = this.AnalogFineAmplitudeWeights.toStruct();
+            s.DigitalAmplitudeWeights = this.DigitalAmplitudeWeights.toStruct();
+            s.FinalCalibrationWeights = this.FinalCalibrationWeights.toStruct();
+        end
+    end
+end
+
diff --git a/Matlab/DigitalCalibrationFormat.m b/Matlab/DigitalCalibrationFormat.m
new file mode 100644
index 0000000..81cfc34
--- /dev/null
+++ b/Matlab/DigitalCalibrationFormat.m
@@ -0,0 +1,28 @@
+classdef DigitalCalibrationFormat
+    
+    properties
+        MeasuredData
+    end
+
+    methods        
+        function s = toStruct(this)
+            s.AmplitudeCalData = this.MeasuredData;
+        end
+
+        function gain = channelGain(this)
+            fReceive = max(abs(fft(this.MeasuredData)));
+            gain = max(fReceive) ./ fReceive;
+        end
+
+        function channelWeights = getChannelPhaseOffset(this)
+            phase = deg2rad(0:360);
+            st_vec = [ones(size(phase)); exp(1i*phase)];
+            sig = this.MeasuredData*conj(st_vec);
+            sigfft = fft(sig);
+            af = max(abs(sigfft));
+            [~, ind_phase] = max(af);
+            channelWeights = st_vec(:,ind_phase);
+        end
+    end
+end
+
diff --git a/Matlab/GainProfile.mat b/Matlab/GainProfile.mat
new file mode 100644
index 0000000..8120a41
Binary files /dev/null and b/Matlab/GainProfile.mat differ
diff --git a/Matlab/Phaser_steeringAngle_rev4.m b/Matlab/Phaser_steeringAngle_rev4.m
new file mode 100644
index 0000000..7ee3c80
--- /dev/null
+++ b/Matlab/Phaser_steeringAngle_rev4.m
@@ -0,0 +1,141 @@
+clear; close all;
+warning('off','MATLAB:system:ObsoleteSystemObjectMixin');
+
+% Key Parameters
+signal_freq = 10.145e9;  % this is the HB100 frequency
+%signal_freq = findTxFrequency();
+plutoURI = 'ip:192.168.2.1';
+phaserURI = 'ip:phaser.local';
+run_calibration = true;
+use_calibration = true;
+
+%% Run Calibration
+if run_calibration == true
+    CalibrationData = calibrationRoutine(signal_freq);
+    finalcalweights = CalibrationData.CalibrationWeights.FinalCalibrationWeights;
+    save("calibration_data.mat", "finalcalweights");
+else
+    if isfile("calibration_data.mat") == true
+        load("calibration_data.mat", "finalcalweights");
+    else
+        use_calibration = false;
+    end
+end
+
+%% Load the measured gain profiles
+load('GainProfile.mat','subArray1_NormalizedGainProfile','subArray2_NormalizedGainProfile','gaincode'); 
+
+% Setup the pluto
+rx = setupPluto(plutoURI);
+
+% Setup the phaser
+bf = setupPhaser(rx,phaserURI,signal_freq);
+bf.RxPowerDown(:) = 0;
+bf.RxGain(:) = 127;
+
+% Create the model of the phaser, which consists of 2 4-element subarrays    
+nElements = 4;
+subarrayModel = phased.ULA('NumElements',nElements,'ElementSpacing',bf.ElementSpacing);
+nSubarrays = 2;
+arrayModel = phased.ReplicatedSubarray("Subarray",subarrayModel,"GridSize",[1,nSubarrays],'SubarraySteering','Custom');
+
+% Create the steering vector for the subarray and the main array
+subarraySteer = phased.SteeringVector("SensorArray",subarrayModel,'NumPhaseShifterBits',7);
+arraySteer = phased.SteeringVector("SensorArray",arrayModel);
+
+% Set the analog and digital calibration weights. To see the uncalibrated
+% array factor, all weights are set to 1.
+if use_calibration == true
+    % These will need to be set based on the calibration routine
+    analogWeights = finalcalweights.AnalogWeights;
+    digitalWeights = finalcalweights.DigitalWeights;
+else
+    % Weights all equal to 1 is an uncalibrated antenna
+    analogWeights = ones(4,2);
+    digitalWeights = ones(2,1);
+end
+
+%% Sweep the steering angle and capture data
+steeringAngle = -90 : 90;
+ArrayFactor = zeros(size(steeringAngle));
+for ii = 1 : numel(steeringAngle)
+    currentangle = steeringAngle(ii);
+
+    % Get the steering weights for the subarray elements
+    subarraysteerweights = subarraySteer(signal_freq,currentangle);
+    adjustedsubarrayweights = subarraysteerweights .* analogWeights;
+
+    % Get element phase from steering weights
+    sub1weights = adjustedsubarrayweights(:,1);
+    sub2weights = adjustedsubarrayweights(:,2);
+    sub1phase = rad2deg(angle(sub1weights));
+    sub2phase = rad2deg(angle(sub2weights));
+    phases = [sub1phase',sub2phase'];
+    phases = phases - phases(1);
+    phases = wrapTo360(phases);
+
+    % Get the element gain from steering weights
+    sub1gain = mag2db(abs(sub1weights));
+    sub2gain = mag2db(abs(sub2weights));
+    calibGainCode = zeros(1,8);
+    for nch = 1 : 4
+        % Set gain codes for array 1
+        xp = sub1gain(nch);
+        if xp == -Inf
+            calibGainCode(nch) = 0;
+        else
+            calibGainCode(nch) = round(interp1(subArray1_NormalizedGainProfile(:,nch),gaincode,xp));
+        end
+    
+        % Set gain codes for array 2
+        xp = sub2gain(nch);
+        if xp == -Inf
+            calibGainCode(nch+4) = 0;
+        else
+            calibGainCode(nch+4) = round(interp1(subArray2_NormalizedGainProfile(:,nch),gaincode,xp));
+        end
+        
+        % Make sure the values of the gain code is always <= 127
+        calibGainCode(calibGainCode>127 | isnan(calibGainCode)) = 127;
+    end
+
+    % Set element phase and gain
+    bf.RxPhase(:) = phases;
+    bf.RxGain(:) = calibGainCode;
+
+    % Collect data
+    bf.LatchRxSettings();
+    receivedSig_HW = rx();
+
+    % Get the steering weights for the digital channels, we need to flip
+    % them before applying the steering vector because the receive data is
+    % out of order
+    flipdigitalweights = [digitalWeights(2);digitalWeights(1)];
+    arraysteerweights = arraySteer(signal_freq,currentangle);
+    adjustedarrayweights = arraysteerweights .* flipdigitalweights;
+    adjustedarrayweights = [adjustedarrayweights(2);adjustedarrayweights(1)];
+
+    % Apply the digital steering weights
+    receivedSig_HW_sum = receivedSig_HW * conj(adjustedarrayweights);
+
+    % Get the array factor
+    receivedFFT = fft(receivedSig_HW_sum);
+    ArrayFactor(ii) = (max(abs(receivedFFT)));
+end
+
+%% Compare the measured array factor and model
+[~,ind] = max(ArrayFactor);
+EmitterAz = steeringAngle(ind);
+figure(101)
+subarrayWeights = conj(subarraySteer(signal_freq,EmitterAz));
+arrayWeights = arraySteer(signal_freq,EmitterAz);
+pattern(arrayModel,signal_freq,-90:90,0,'CoordinateSystem', ...
+    'Rectangular','Type','powerdb','ElementWeights',[subarrayWeights,subarrayWeights],'Weights',arrayWeights)
+hold on;
+
+% Plot the measured data and the model
+plot(steeringAngle,mag2db(ArrayFactor./max(abs(ArrayFactor))))
+
+%%
+bf.release();
+rx.release();
\ No newline at end of file
diff --git a/Matlab/calibrationRoutine.m b/Matlab/calibrationRoutine.m
new file mode 100644
index 0000000..4340ebf
--- /dev/null
+++ b/Matlab/calibrationRoutine.m
@@ -0,0 +1,352 @@
+function CalibrationData = calibrationRoutine(fc_hb100)
+
+% Setup:
+%
+% Place the HB100 in front of the Phaser - 0 degree azimuth. It should be
+% sufficiently far from the antenna so that the wavefront is approximately
+% planar.
+%
+% Notes:
+% 
+% The amplitude and phase of all 8 analog channels as well as the two
+% digital channels must be calibrated so that the amplitudes in each
+% channel are the same and the phase shifts align as expected during
+% steering.
+%
+% The order of the calibration steps is described below.
+%
+% 1. Course Analog Amplitude Calibration: The first calibration step is an
+% initial amplitude calibration so that the amplitude in each of the 4
+% elements in each of the 2 subarrays are the same. This is done by
+% independently measuring the signal amplitude in each channel when the gain
+% is set to maximum.The gain is set so that all of the channel signal
+% amplitudes match the lowest signal amplitude channel.
+%
+% 2. Analog Phase Calibration: The second calibration step is to account for
+% the different phase shift in each analog channel. For each subarray, 
+% select a reference channel. Then measure the phase of other channels
+% within the subarray with respect to the reference channel. The search is
+% brute force, by sweeping all the possible phases on the chip.
+%
+% 3. Fine Analog Amplitude Calibration: Changing the default phase shift in
+% each analog channel in step (2) causes the amplitude of the signal in
+% each channel to change. Therefore, the amplitude calibration is performed
+% again.
+%
+% 4. Digital Amplitude Calibration: Once the analog channels have been
+% calibrated, the difference in amplitude between the two digital channels
+% is measured. Whenever the digital channel signals are combined, an
+% amplitude adjustment is applied based on this different in amplitude.
+%
+% 5. Digital Phase Calibration: The signal phase into each of the two
+% digital channels will not necessarily align. The phase difference is
+% measured by shifting the phase of channel 2 until the amplitude of the
+% combined signal is maximized. This phase shift is applied when the two
+% digital signals are combined.
+
+%% Setup
+
+% Setup antenna and model
+[rx,bf,~] = setupAntenna(fc_hb100);
+
+% Setup calibration data storage object
+CalibrationData = CalibrationDataFormat();
+
+% Setup the inital analog and digital antenna weights. These are the values
+% that get adjusted during calibration. The amplitude and phase of the weights
+% represent the gain and phase adjustments that will be made on each channel.
+% Initially there is no phase shift or amplitude adjustment.
+analogWeights = ones(4,2);
+digitalWeights = ones(2,1);
+
+% Store uncalibrated weights for later analysis
+CalibrationData.CalibrationWeights.UncalibratedWeights.AnalogWeights = analogWeights;
+CalibrationData.CalibrationWeights.UncalibratedWeights.DigitalWeights = digitalWeights;
+
+%% Course Amplitude Calibration
+% The amplitude in each analog channel is measured when the gain is
+% set to the maximum level. Gain in each channel is adjusted to account
+% for differences.
+
+% Set the gain of all 8 channels to maximum.
+bf.RxGain(:) = 127;
+
+% Setup data variables
+nCapture = 20;
+rx1data = zeros(1024,nCapture,4);
+rx2data = zeros(1024,nCapture,4);
+
+for nCh = 1:4
+    % Turn off all the channels.
+    bf.RxPowerDown(:) = 1;
+    
+    % Turn on channels under test (subarray 1 and subarray 2)
+    bf.RxPowerDown(nCh) = 0; 
+    bf.RxPowerDown(nCh+4) = 0;
+    bf.LatchRxSettings();
+
+    % Capture multiple snapshots
+    for ncapture = 1 : nCapture        
+        rx();
+        receivedSig = rx();
+        rx1data(:,ncapture,nCh) = receivedSig(:,2);
+        rx2data(:,ncapture,nCh) = receivedSig(:,1);
+    end
+end
+
+% Turn all channels back on
+bf.RxPowerDown(:) = 0;
+
+% Calculate the average amplitude for each element in each subarray
+avgamp1 = mean(max(abs(fft(rx1data))));
+avgamp2 = mean(max(abs(fft(rx2data))));
+
+% Reshape into a vector
+avgamp1 = reshape(avgamp1,[1 4]);
+avgamp2 = reshape(avgamp2,[1 4]);
+
+% Scale each amplitude value to the minimum measured amplitude value
+sub1gain = min(avgamp1) ./ avgamp1;
+sub2gain = min(avgamp2) ./ avgamp2;
+
+% Plot the analog amplitude calibration values
+helperPlotElementGainCalibration(avgamp1,sub1gain,avgamp2,sub2gain);
+
+% Update analog weights with array gain calibration values
+analogWeights = analogWeights .* [sub1gain',sub2gain'];
+
+% Store the calibration data for later analysis
+CalibrationData.AnalogAmplitudeCalibration.CourseSubarray1Data = rx1data;
+CalibrationData.AnalogAmplitudeCalibration.CourseSubarray2Data = rx2data;
+CalibrationData.CalibrationWeights.AnalogCourseAmplitudeWeights.AnalogWeights = analogWeights;
+
+%% Calibrate phase
+% Place the emitter in front of the phased array.
+% Turn on two adjacent channel. On the first channel put 0 phase, and on
+% the adjacent channel sweep the phase shift setting from 0 to 360. Minimum
+% signal strength occurs when the phase offset between channels is 180
+% degrees. The calibration value is set such that the phase shift is
+% actually 180 degrees when shifter is set to 180 degrees.
+
+% Setup phase steps to test and variables to hold data
+PhaseResolutionBits = 7;
+phase_step_size = 360 / (2 ^ PhaseResolutionBits);
+channel2phase = double(0:phase_step_size:360);
+rawdata1 = zeros(rx.SamplesPerFrame,numel(channel2phase),3);
+rawdata2 = zeros(rx.SamplesPerFrame,numel(channel2phase),3);
+
+% Set the reference channel and channels under test
+referencechannel = 1;
+testchannels = 2:4;
+
+% set the receiver gain based on the amplitude calibration. The gain codes
+% are set based on data that was previously captured to determine gain code
+% correspondance to actual gain.
+bf.RxGain = helperGainCodes(analogWeights);
+
+% Loop through each channel under test
+for nChannel = testchannels
+    % Set all phase shifts to 0
+    bf.RxPhase(:) = 0;
+
+    % Power down all elements
+    bf.RxPowerDown(:) = 1;
+
+    % Turn on element 1 and element under test subarray 1
+    bf.RxPowerDown(referencechannel) = 0;
+    bf.RxPowerDown(nChannel) = 0;
+
+    % Turn on element 1 and element under test subarray 2
+    bf.RxPowerDown(referencechannel+4) = 0;
+    bf.RxPowerDown(nChannel+4) = 0;
+    bf.LatchRxSettings();
+
+    % Loop through each phase offset, set element under test phase shifter,
+    % capture data.
+    for ii = 1 : numel(channel2phase)
+        % Set phase for subarray 1 and 2 element
+        bf.RxPhase(nChannel) = channel2phase(ii);
+        bf.RxPhase(nChannel+4) = channel2phase(ii);
+        bf.LatchRxSettings();
+
+        % Capture and store data
+        rx();
+        receivedSig_HW = rx();
+        rawdata1(:,ii,nChannel-1) = receivedSig_HW(:,2);
+        rawdata2(:,ii,nChannel-1) = receivedSig_HW(:,1);
+    end
+end
+
+% Get signal amplitude for both subarrays
+amp1 = max(abs(fft(rawdata1)),[],1);
+amp2 = max(abs(fft(rawdata2)),[],1);
+
+% Reshape into a 2d array which is Number Phases x Number Elements
+sub1amplitudes = reshape(amp1,[numel(channel2phase) 3]);
+sub2amplitudes = reshape(amp2,[numel(channel2phase) 3]);
+
+% Get location of min amplitude
+[~,phase1Idx] = min(sub1amplitudes,[],1);
+[~,phase2Idx] = min(sub2amplitudes,[],1);
+
+% Get calibration phase by substracting 180 from the minimum phase location
+cal1phase = [0;wrapTo360(channel2phase(phase1Idx)-180)'];
+cal2phase = [0;wrapTo360(channel2phase(phase2Idx)-180)'];
+
+% Plot the analog phase shift information
+helperPlotAnalogPhaseCalibration(channel2phase,sub1amplitudes,sub2amplitudes);
+
+% Convert calibration phase to a phase shift for the element weights, and
+% update analog weights
+phaseSteer = exp(deg2rad([cal1phase,cal2phase])*1i);
+analogWeights = analogWeights .* phaseSteer;
+
+% Store recorded data for later analysis
+CalibrationData.AnalogPhaseCalibration.PhaseSetting = channel2phase;
+CalibrationData.AnalogPhaseCalibration.Subarray1Measurements = rawdata1;
+CalibrationData.AnalogPhaseCalibration.Subarray2Measurements = rawdata2;
+CalibrationData.CalibrationWeights.AnalogPhaseWeights.AnalogWeights = analogWeights;
+
+%% Fine Grain Calibration Data
+
+% The amplitude in each analog channel is measured when the gain is
+% set to the maximum level with the new default phase shifts applied.
+% Gain in each channel is adjusted to account for differences. This is the
+% exact same process that was used for the analog course amplitude
+% calibration, although now the phase shifter values are set to the default
+% determined in the previous section.
+
+% Set the gain of all 8 channels to maximum.
+bf.RxGain(:) = 127;
+
+% Set the phase shifts based on the new analog weights
+phaseshifts = wrapTo360(rad2deg(angle(analogWeights)));
+bf.RxPhase(:) = [phaseshifts(:,1)',phaseshifts(:,2)'];
+
+% Setup data variables
+nCapture = 20;
+rx1data = zeros(1024,nCapture,4);
+rx2data = zeros(1024,nCapture,4);
+
+for nCh = 1:4
+    % Turn off all the channels.
+    bf.RxPowerDown(:) = 1;
+    
+    % Turn on channels under test (subarray 1 and subarray 2)
+    bf.RxPowerDown(nCh) = 0; 
+    bf.RxPowerDown(nCh+4) = 0;
+    bf.LatchRxSettings();
+
+    % Capture multiple snapshots
+    for ncapture = 1 : nCapture
+        rx();
+        receivedSig = rx();
+        rx1data(:,ncapture,nCh) = receivedSig(:,2);
+        rx2data(:,ncapture,nCh) = receivedSig(:,1);
+    end
+end
+
+% Turn all channels back on
+bf.RxPowerDown(:) = 0;
+
+% Calculate the average amplitude for each element in each subarray
+avgamp1 = mean(max(abs(fft(rx1data))));
+avgamp2 = mean(max(abs(fft(rx2data))));
+
+% Reshape into a vector
+avgamp1 = reshape(avgamp1,[1 4]);
+avgamp2 = reshape(avgamp2,[1 4]);
+
+% Scale each amplitude value to the minimum measured amplitude value
+sub1gain = min(avgamp1) ./ avgamp1;
+sub2gain = min(avgamp2) ./ avgamp2;
+
+% Normalize the weights back to amplitude 1
+analogWeights = analogWeights ./ abs(analogWeights);
+
+% Update analog weights with array gain calibration values
+analogWeights = analogWeights .* [sub1gain',sub2gain'];
+
+% Store the calibration data for later analysis
+CalibrationData.AnalogAmplitudeCalibration.FineSubarray1Data = rx1data;
+CalibrationData.AnalogAmplitudeCalibration.FineSubarray2Data = rx2data;
+CalibrationData.CalibrationWeights.AnalogFineAmplitudeWeights.AnalogWeights = analogWeights;
+
+%% Calibrate the two digital channels on Pluto
+% After the analog channels have been calibrated, collect data with
+% analog calibration value settings and calibrate the digital channels for
+% amplitude and phase.
+
+% Collect data with analog weight settings
+analogcaldata = helperSteerAnalog(bf,rx,analogWeights);
+
+% Store calibration data for later analysis
+CalibrationData.DigitalCalibration.MeasuredData = analogcaldata;
+
+% Calculate the channel calibration gain by normalizing the channel
+% amplitudes to the max channel amplitude
+channelamplitude = max(abs(fft(analogcaldata)));
+gain = max(channelamplitude) ./ channelamplitude;
+digitalWeights = digitalWeights .* gain';
+
+% Save the weights for later analysis
+CalibrationData.CalibrationWeights.DigitalAmplitudeWeights.AnalogWeights = analogWeights;
+CalibrationData.CalibrationWeights.DigitalAmplitudeWeights.DigitalWeights = digitalWeights;
+
+% Calculate the channel calibration phase by steering channel 2 from 0 to
+% 360 and finding the maximum combined channel amplitude
+channel2phase = deg2rad(0:360);
+st_vec = [ones(size(channel2phase)); exp(1i*channel2phase)];
+sig = analogcaldata*conj(st_vec);
+sigfft = fft(sig);
+combinedamplitude = max(abs(sigfft));
+[~, phaseIdx] = max(combinedamplitude);
+
+% Plot the digital phase offset pattern and calibration value
+ax = axes(figure); hold(ax,"on");
+title(ax,"Digital Phase Calibration"); ylabel(ax,"dB"); xlabel(ax,"Channel 2 Phase Offset (deg)");
+plot(ax,channel2phase,combinedamplitude,"DisplayName","Combined Channel Power");
+scatter(ax,channel2phase(phaseIdx),combinedamplitude(phaseIdx),"DisplayName","Selected Phase Offset");
+legend('location','southeast');
+
+% Set the new digital weights
+defaultsteervec = st_vec(:,phaseIdx);
+digitalWeights = digitalWeights .* defaultsteervec;
+
+% Save the final calibration weights
+CalibrationData.CalibrationWeights.FinalCalibrationWeights.AnalogWeights = analogWeights;
+CalibrationData.CalibrationWeights.FinalCalibrationWeights.DigitalWeights = digitalWeights;
+
+%% Release Hardware
+bf.release(); delete(bf);
+rx.release(); delete(rx);
+
+%% Plot final data: uncalibrated data vs. calibrated data vs. simulated data
+
+% Use a custom class called AntennaInteractor to help steer the antenna
+initialcalvalues = CalibrationData.CalibrationWeights.UncalibratedWeights;
+antennaInteractor = AntennaInteractor(fc_hb100,initialcalvalues);
+steerangles = -90:0.5:90;
+
+% Capture the pattern with the initial calibration values
+[initialpattern,~] = antennaInteractor.capturePattern(steerangles);
+initialamp = mag2db(helperGetAmplitude(initialpattern));
+
+% Capture the pattern with the final calibration values
+finalcalvalues = CalibrationData.CalibrationWeights.FinalCalibrationWeights;
+antennaInteractor.updateCalibration(finalcalvalues);
+[finalpattern,~] = antennaInteractor.capturePattern(steerangles);
+finalamp = mag2db(helperGetAmplitude(finalpattern));
+
+% Simulate the expected antenna pattern
+rxpos = [0;0;0]; % phaser at 0
+txpos = [0;10;0]; % transmitter 10 meters away at 0 boresight
+simpattern = mag2db(helperSimulateAntennaSteering(fc_hb100,rxpos,txpos,steerangles));
+
+% Plot data
+patternax = axes(figure); hold(patternax,"on"); legend(patternax);
+xlabel(patternax,"Steer Angle (deg)"); ylabel(patternax,"Amplitude (dB)"); title(patternax,"Calibration Effect");
+plot(patternax,steerangles,initialamp - max(initialamp),'DisplayName','Uncalibrated Array Factor');
+plot(patternax,steerangles,finalamp - max(finalamp),'DisplayName','Calibrated Array Factor');
+plot(patternax,steerangles,simpattern - max(simpattern),'DisplayName','Simulated Array Factor');
+end
diff --git a/Matlab/calibration_data.mat b/Matlab/calibration_data.mat
new file mode 100644
index 0000000..49eec3a
Binary files /dev/null and b/Matlab/calibration_data.mat differ
diff --git a/Matlab/cleanupAntenna.m b/Matlab/cleanupAntenna.m
new file mode 100644
index 0000000..7038674
--- /dev/null
+++ b/Matlab/cleanupAntenna.m
@@ -0,0 +1,4 @@
+function cleanupAntenna(rx,bf)
+    rx.release();
+    bf.release();
+end
\ No newline at end of file
diff --git a/Matlab/finalcalweights.mat b/Matlab/finalcalweights.mat
new file mode 100644
index 0000000..1b0e5e5
Binary files /dev/null and b/Matlab/finalcalweights.mat differ
diff --git a/Matlab/findTxFrequency.m b/Matlab/findTxFrequency.m
new file mode 100644
index 0000000..6fc2cd2
--- /dev/null
+++ b/Matlab/findTxFrequency.m
@@ -0,0 +1,58 @@
+function txFrequency = findTxFrequency()
+% Setup:
+%
+% Place the HB100 in front of the Phaser - 0 degree azimuth.
+%
+% Notes:
+% 
+% The frequency of the phaser is swept between 10 and 10.5 GHz. The peak
+% frequency is measured. This is assumed to be the frequency of the HB100.
+%
+
+% Setup the frequencies to scan.
+f_start = 10.e9;
+f_stop = 10.5e9;
+f_step = 5e6; 
+fvec = f_start : f_step : f_stop;
+
+% Setup the antenna, setting the frequency to f_start.
+[rx,bf,~] = setupAntenna(f_start);
+
+% Setup variables for capturing scan amplitude and frequency.
+full_ampl = [];
+full_freqs = [];
+N = rx.SamplesPerFrame;
+
+% Loop through the frequency vector and save receive data at each center
+% frequency.
+for centerfrequency = fvec
+    % The LO is set to 4x the Frequency setting.
+    bf.Frequency = (centerfrequency + rx.CenterFrequency)/4;
+    rx();
+
+    % Capture the data, sum the channels, and convert to the frequency domain
+    data = rx();
+    data = sum(data,2);
+    famplitude = mag2db(1/N * fftshift(abs(fft(data))));
+    full_ampl = [full_ampl, famplitude];
+
+    % Get the frequency span for this data sample
+    df = rx.SamplingRate/N;
+    band = (-rx.SamplingRate/2 : df : rx.SamplingRate/2 - df);
+    freqspan = centerfrequency-band;
+    full_freqs = [full_freqs, freqspan'];
+end
+
+% Get the max amplitude for each measurement
+[maxamplitudes,maxidxs] = max(full_ampl);
+
+% Get the max amplitude in the whole dataset
+[~,maxframeidx] = max(maxamplitudes);
+maxdatapointidx = maxidxs(maxframeidx);
+
+txFrequency = full_freqs(maxdatapointidx,maxframeidx);
+
+ax = axes(figure);
+plot(ax,full_freqs/1e9,full_ampl); xlabel('Frequency (GHz)'); ylabel('Amplitude (dB)');
+title(ax,['HB100 Frequency = ' num2str(txFrequency/1e9) 'GHz']);
+end
diff --git a/Matlab/helperCalculateAmplitude.m b/Matlab/helperCalculateAmplitude.m
new file mode 100644
index 0000000..a3c7e26
--- /dev/null
+++ b/Matlab/helperCalculateAmplitude.m
@@ -0,0 +1,13 @@
+function amplitude = helperCalculateAmplitude(data,maxvalue)
+    % Get the signal amplitude
+
+    % Scale data
+    datascaled = data / maxvalue;
+    [nsamples,~] = size(data);
+
+    % Convert the signal to the frequency domain
+    fexampledata = mag2db(abs(fft(datascaled)) / nsamples);
+
+    % Amplitude is the largest frequency value
+    amplitude = max(fexampledata);
+end
diff --git a/Matlab/helperGainCodes.m b/Matlab/helperGainCodes.m
new file mode 100644
index 0000000..df15f21
--- /dev/null
+++ b/Matlab/helperGainCodes.m
@@ -0,0 +1,36 @@
+function calibGainCode = helperGainCodes(analogWeights)
+    % get gain for each subarray
+    sub1weights = analogWeights(:,1);
+    sub2weights = analogWeights(:,2);
+    sub1gain = getGain(sub1weights);
+    sub2gain = getGain(sub2weights);
+
+    % Loading the measured gain profiles
+    load('GainProfile.mat','subArray1_NormalizedGainProfile','subArray2_NormalizedGainProfile','gaincode'); 
+    
+    calibGainCode = zeros(1,8);
+    for nch = 1 : 4
+    
+        xp = sub1gain(nch);
+        if xp == -Inf
+            calibGainCode(nch) = 0;
+        else
+            calibGainCode(nch) = round(interp1(subArray1_NormalizedGainProfile(:,nch),gaincode,xp));
+        end
+    
+        xp = sub2gain(nch);
+        if xp == -Inf
+            calibGainCode(nch+4) = 0;
+        else
+            calibGainCode(nch+4) = round(interp1(subArray2_NormalizedGainProfile(:,nch),gaincode,xp));
+        end
+    end
+    % Make sure the values of the gain code is always <= 127
+    calibGainCode(calibGainCode>127 | isnan(calibGainCode)) = 127;
+
+    function gain = getGain(weights)
+        amp = abs(weights);
+        gain = mag2db(amp);
+    end
+end
+
diff --git a/Matlab/helperGetAmplitude.m b/Matlab/helperGetAmplitude.m
new file mode 100644
index 0000000..489903a
--- /dev/null
+++ b/Matlab/helperGetAmplitude.m
@@ -0,0 +1,3 @@
+function [amp] = helperGetAmplitude(signal)
+    amp = max(abs(fft(signal)));
+end
\ No newline at end of file
diff --git a/Matlab/helperGetPhase.m b/Matlab/helperGetPhase.m
new file mode 100644
index 0000000..7fade84
--- /dev/null
+++ b/Matlab/helperGetPhase.m
@@ -0,0 +1,8 @@
+function phase = helperGetPhase(signal)
+    fsig = fft(signal);
+    [~,ampidx] = max(fsig);
+    phase = zeros(1,numel(ampidx));
+    for i = 1:numel(ampidx)
+        phase(i) = rad2deg(angle(fsig(ampidx(i),i)));
+    end
+end
\ No newline at end of file
diff --git a/Matlab/helperPlotAnalogPhaseCalibration.m b/Matlab/helperPlotAnalogPhaseCalibration.m
new file mode 100644
index 0000000..dfb3ecf
--- /dev/null
+++ b/Matlab/helperPlotAnalogPhaseCalibration.m
@@ -0,0 +1,19 @@
+function helperPlotAnalogPhaseCalibration(phasesetting,sub1amplitudes,sub2amplitudes)
+    % Visualize the element-wise phase calibration data for the entire
+    % array.
+    figure; tiledlayout(1,2); nexttile();
+    helperPlotPhaseSubarrayData("Subarray 1",phasesetting,mag2db(sub1amplitudes)); nexttile();
+    helperPlotPhaseSubarrayData("Subarray 2",phasesetting,mag2db(sub2amplitudes));
+end
+
+function helperPlotPhaseSubarrayData(name, phasesetting, phaseOffsetAmplDb)
+    % Visualize the element-wise phase calibration data for a subarray.
+    lines = plot(phasesetting, phaseOffsetAmplDb);
+    lines(1).DisplayName = "Element 2";
+    lines(2).DisplayName = "Element 3";
+    lines(3).DisplayName = "Element 4";
+    title([name,' - Phase Calibration Data'])
+    ylabel("Amplitude (dB) Element X + Element 1")
+    xlabel("Test Element Phase Shift (deg)")
+    legend("Location","southoutside")
+end
\ No newline at end of file
diff --git a/Matlab/helperPlotElementGainCalibration.m b/Matlab/helperPlotElementGainCalibration.m
new file mode 100644
index 0000000..8075db7
--- /dev/null
+++ b/Matlab/helperPlotElementGainCalibration.m
@@ -0,0 +1,41 @@
+function [array1gaincal,array2gaincal] = helperPlotElementGainCalibration(sub1meanamp,sub1gaincal,sub2meanamp,sub2gaincal)
+    % Calculate and visualize the element-wise amplitude calibration for
+    % both subarrays.
+
+    sub1meanamp = mag2db(sub1meanamp);
+    sub1gaincal = mag2db(sub1gaincal);
+    sub2meanamp = mag2db(sub2meanamp);
+    sub2gaincal = mag2db(sub2gaincal);
+    
+    % Setup figure
+    figure; tiledlayout(1,2); a = nexttile();
+
+    % Calculate and plot the gain calibration for subarray 1
+    array1gaincal = helperElementSubarrayGainCalibration(a,'Subarray 1',sub1meanamp, sub1gaincal); a = nexttile();
+    
+    % Calculate and plot the gain calibration for subarray 2
+    array2gaincal = helperElementSubarrayGainCalibration(a, 'Subarray 2', sub2meanamp, sub2gaincal);
+end
+
+function arraygaincal = helperElementSubarrayGainCalibration(ax,name,amplitudes,arraygaincal)
+    % Calculate and visualize the element-wise amplitude calibration for
+    % one subarray.
+
+    hold(ax,"on");
+
+    % Normalize amplitude for each element in the array
+    dbNormAmplitudes = amplitudes - max(amplitudes);
+    
+    % Plot normalized amplitudes and gain adjustments
+    b = bar(ax,[dbNormAmplitudes',arraygaincal'],'stacked');
+    b(1).DisplayName = "Initial Normalized Amplitude";
+    b(2).DisplayName = "Gain Adjustment";
+
+    % Plot a line showing the final amplitude of all elements
+    plot(ax,[0,5],[min(dbNormAmplitudes),min(dbNormAmplitudes)],"DisplayName","Final Element Amplitude","LineWidth",2,"Color","k")
+
+    xlabel('Antenna Element')
+    ylabel('dB');
+    title([name ' - Gain Calibration'])
+    legend('Location','southoutside')
+end
\ No newline at end of file
diff --git a/Matlab/helperPlotGainCalibration.m b/Matlab/helperPlotGainCalibration.m
new file mode 100644
index 0000000..4002d6a
--- /dev/null
+++ b/Matlab/helperPlotGainCalibration.m
@@ -0,0 +1,18 @@
+function gain = helperPlotGainCalibration(name, amplitudes)
+% Normalize amplitude for each element in the array
+dbNormAmplitudes = mag2db(amplitudes./max(amplitudes));
+
+% Calculate gain adjustment
+gain = -mag2db(amplitudes./min(amplitudes));
+
+% Plot normalized amplitudes and gain adjustments
+b = bar([dbNormAmplitudes',gain'],'stacked');
+b(1).DisplayName = "Initial Normalized Amplitude";
+b(2).DisplayName = "Gain Adjustment";
+xlabel('Antenna Element')
+ylabel('dB');
+title([name ' - Gain Calibration'])
+legend('Location','southoutside')
+
+end
+
diff --git a/Matlab/helperSimulateAntennaSteering.m b/Matlab/helperSimulateAntennaSteering.m
new file mode 100644
index 0000000..f6c8b7c
--- /dev/null
+++ b/Matlab/helperSimulateAntennaSteering.m
@@ -0,0 +1,86 @@
+function [rxamp,rxphase] = helperSimulateAntennaSteering(fctransmit,rxpos,txpos,steerangle,analogweight,digitalweight,nullangle)
+% This helper function simulates antenna steering for the ADI workshop.
+% analogweight and digitalweight are optional inputs that default to 1.
+arguments
+    fctransmit (1,1) double
+    rxpos (3,1) double
+    txpos (3,1) double
+    steerangle (1,:) double
+    analogweight (4,2) double = ones(4,2)
+    digitalweight (2,1) double = ones(2,1)
+    nullangle = []
+end
+
+% Setup the transmitter
+txelement = phased.IsotropicAntennaElement;
+radiator = phased.Radiator("Sensor",txelement,"OperatingFrequency",fctransmit);
+
+% Setup the Phased Array Receiver based on the Phaser specs
+frange = [10.0e9 10.5e9];
+sampleRate = 30e6;
+nsamples = 1024;
+element = phased.IsotropicAntennaElement('FrequencyRange',frange);
+hRange = freq2wavelen(frange);
+spacing = hRange(2)/2;
+subarrayElements = 4;
+subarray = phased.ULA('Element',element,'NumElements',subarrayElements,'ElementSpacing',spacing);
+array = phased.ReplicatedSubarray('Subarray',subarray,'GridSize',[1,2],"SubarraySteering","Custom");
+collector = phased.Collector("Sensor",array,"OperatingFrequency",fctransmit,"WeightsInputPort",true);
+
+% CW signal is used
+signal = ones(nsamples,1);
+
+% Set up a channel for radiating signal
+channel = phased.FreeSpace("OperatingFrequency",fctransmit,"SampleRate",sampleRate);
+
+% Setup geometry
+rxvel = [0;0;0];
+txvel = [0;0;0];
+[~,ang] = rangeangle(rxpos,txpos);
+
+% Radiate signal
+sigtx = radiator(signal,ang);
+
+% Propagate signal
+sigtx = channel(sigtx,txpos,rxpos,txvel,rxvel);
+
+% Create steering vector for the two subarray channels
+steervec = phased.SteeringVector("SensorArray",array);
+substeervec = phased.SteeringVector("SensorArray",subarray);
+
+% Receive signal while steering beam
+rxphase = zeros(1,numel(steerangle));
+rxamp = zeros(1,numel(steerangle));
+for steer = steerangle
+    % Create the subarray weights.
+    singlesubweight = substeervec(fctransmit,steer);
+
+    % Create the replicated array weights
+    repweight = steervec(fctransmit,steer);
+
+    % insert a null if a null angle was passed in
+    if ~isempty(nullangle)
+        % Null the subarray
+        nullsubweight = substeervec(fctransmit,nullangle);
+        singlesubweight = getNullSteer(singlesubweight,nullsubweight);
+
+        % null the replicated array
+        nullrepweight = steervec(fctransmit,nullangle);
+        repweight = getNullSteer(repweight,nullrepweight);
+    end
+
+    % Create the full subarray weights
+    subweight = [singlesubweight,singlesubweight] .* analogweight;
+
+    % Receive the signal
+    sigreceive = collector(sigtx,[0;0],repweight,subweight) * digitalweight;
+    rxphase(steer == steerangle) = mean(rad2deg(angle(sigreceive)));
+    rxamp(steer == steerangle) = mean(abs(sigreceive));
+end
+
+end
+
+function nullsteer = getNullSteer(steerweight,nullweight)
+    rn = nullweight'*steerweight/(nullweight'*nullweight);
+    nullsteer = steerweight-nullweight*rn;
+end
\ No newline at end of file
diff --git a/Matlab/helperSimulateDisabledElement.m b/Matlab/helperSimulateDisabledElement.m
new file mode 100644
index 0000000..2ab18d9
--- /dev/null
+++ b/Matlab/helperSimulateDisabledElement.m
@@ -0,0 +1,9 @@
+function simpattern = helperSimulateDisabledElement(fc,steerangles,iEl)
+    % Simulate what happens when we disable an element on each of the antenna
+    % subarrays.
+    disabletaper = ones(4,2);
+    disabletaper(iEl,:) = 0;
+    rxpos = [0;0;0];
+    txpos = [0;10;0];
+    simpattern = helperSimulateAntennaSteering(fc,rxpos,txpos,steerangles,disabletaper);
+end
\ No newline at end of file
diff --git a/Matlab/helperSimulateMonopulse.m b/Matlab/helperSimulateMonopulse.m
new file mode 100644
index 0000000..65e9514
--- /dev/null
+++ b/Matlab/helperSimulateMonopulse.m
@@ -0,0 +1,15 @@
+function [sumamp,diffamp,phasedelta] = helperSimulateMonopulse(fc_hb100,steerangles)
+    % Simulate the monopulse pattern for the ADI phaser.
+
+    % Get the sum pattern
+    rxpos = [0;0;0];
+    txpos = [0;10;0];
+    [sumamp,sumphase] = helperSimulateAntennaSteering(fc_hb100,rxpos,txpos,steerangles);
+    
+    % The diff pattern is generated by setting the digital weights to
+    % [1;-1]
+    analogweights = ones(4,2);
+    digitalweights = [1;-1];
+    [diffamp,diffphase] = helperSimulateAntennaSteering(fc_hb100,rxpos,txpos,steerangles,analogweights,digitalweights);
+    phasedelta = sign(wrapTo180(sumphase-diffphase));
+end
\ No newline at end of file
diff --git a/Matlab/helperSimulateNull.m b/Matlab/helperSimulateNull.m
new file mode 100644
index 0000000..a09806e
--- /dev/null
+++ b/Matlab/helperSimulateNull.m
@@ -0,0 +1,8 @@
+function pattern = helperSimulateNull(fc_hb100,steerangles,nullangle)
+    % Simulate the pattern with a null at the specified angle
+    analogweights = ones(4,2);
+    digitalweights = ones(2,1);
+    rxpos = [0;0;0];
+    txpos = [0;10;0];
+    [pattern,~] = helperSimulateAntennaSteering(fc_hb100,rxpos,txpos,steerangles,analogweights,digitalweights,nullangle);
+end
\ No newline at end of file
diff --git a/Matlab/helperSourceLocation.m b/Matlab/helperSourceLocation.m
new file mode 100644
index 0000000..1ba136f
--- /dev/null
+++ b/Matlab/helperSourceLocation.m
@@ -0,0 +1,21 @@
+function sourceaz = helperSourceLocation(CalibrationData)
+
+[rx,bf,fc_hb100,phaseCal,channelWeights,steeringVec,~] = setupAntenna(CalibrationData);
+
+% Capture the initial pattern data, find the location of the interferer -
+% assumed to be higher power
+steerangles = -90:0.5:90;
+patternData = helperCapturePattern(steerangles,steeringVec,fc_hb100,bf,rx,phaseCal,channelWeights);
+amp = helperGetAmplitude(patternData);
+smoothamp = smooth(amp,20);
+[~,azidx] = max(smoothamp);
+sourceaz = steerangles(azidx);
+
+ax = axes(figure); hold(ax,"on");
+plot(ax,steerangles,amp,"DisplayName","Scan Pattern")
+scatter(ax,sourceaz,amp(azidx),"DisplayName","Source Location");
+legend();
+
+cleanupAntenna(rx,bf);
+
+end
\ No newline at end of file
diff --git a/Matlab/helperSteerAnalog.m b/Matlab/helperSteerAnalog.m
new file mode 100644
index 0000000..43dcb46
--- /dev/null
+++ b/Matlab/helperSteerAnalog.m
@@ -0,0 +1,27 @@
+function rxdata = helperSteerAnalog(bf,rx,analogWeights)
+    sub1weights = analogWeights(:,1);
+    sub2weights = analogWeights(:,2);
+
+    % Set analog phase shifter
+    sub1phase = getPhase(sub1weights);
+    sub2phase = getPhase(sub2weights);
+    phases = [sub1phase',sub2phase'];
+    if ~isequal(bf.RxPhase,phases)
+        bf.RxPhase(:) = phases;
+    end
+
+    % Set analog gain
+    gainCode = helperGainCodes(analogWeights);
+    if ~isequal(bf.RxGain,gainCode)
+        bf.RxGain(:) = gainCode;
+    end
+
+    % receive data
+    bf.LatchRxSettings();
+    rx();
+    rxdata = rx();
+end
+
+function phase = getPhase(weights)
+    phase = wrapTo360(rad2deg(angle(weights)));
+end
\ No newline at end of file
diff --git a/Matlab/setupAntenna.m b/Matlab/setupAntenna.m
new file mode 100644
index 0000000..85119dc
--- /dev/null
+++ b/Matlab/setupAntenna.m
@@ -0,0 +1,19 @@
+function [rx,bf,phaserModel] = setupAntenna(fc_hb100)
+    % Setup the Pluto, Phaser, and Phaser Model.
+    
+    % Setup the pluto
+    plutoURI = 'ip:192.168.2.1';
+    rx = setupPluto(plutoURI);
+
+    % Setup the phaser
+    phaserURI = 'ip:phaser.local';
+    bf = setupPhaser(rx,phaserURI,fc_hb100);
+    bf.RxPowerDown(:) = 0;
+    bf.RxGain(:) = 127;
+    
+    % Create the model of the phaser    
+    nElements = 4;
+    nSubarrays = 2;
+    subModel = phased.ULA('NumElements',nElements,'ElementSpacing',bf.ElementSpacing);
+    phaserModel = phased.ReplicatedSubarray("Subarray",subModel,"GridSize",[1,nSubarrays]);
+end
\ No newline at end of file
diff --git a/Matlab/setupPhaser.m b/Matlab/setupPhaser.m
new file mode 100644
index 0000000..370ba2b
--- /dev/null
+++ b/Matlab/setupPhaser.m
@@ -0,0 +1,48 @@
+function bf = setupPhaser(rx,phaserURI,fc_hb100)
+%% Configure phaser
+bf = adi.Phaser;
+bf.uri = phaserURI;
+bf.SkipInit = true; % Bypass writing all initial attributes to speed things up
+bf();
+bf.ElementSpacing = 0.014;
+% Put device in Rx mode
+bf.TxRxSwitchControl = {'spi','spi'};
+bf.Mode(:) = {'Disabled'};
+bf.BeamMemEnable(:) = false;
+bf.BiasMemEnable(:) = false;
+bf.PolState(:) = false;
+bf.PolSwitchEnable(:) = false;
+bf.TRSwitchEnable(:) = true;
+bf.ExternalTRPolarity(:) = true;
+
+bf.RxVGAEnable(:) = true;
+bf.RxVMEnable(:) = true;
+bf.RxLNABiasCurrent(:) = 8;
+bf.RxVGABiasCurrentVM(:) = 22;
+
+% Self bias LNAs
+bf.LNABiasOutEnable(:) = false;
+
+% Fire them up
+bf.RxPowerDown(:) = false;
+bf.Mode(:) = {'Rx'};
+
+%% Set up PLL
+bf.Frequency = (fc_hb100 + rx.CenterFrequency) / 4;
+BW = 500e6 / 4; num_steps = 1000;
+bf.FrequencyDeviationRange = BW; % frequency deviation range in H1.  This is the total freq deviation of the complete freq ramp
+bf.FrequencyDeviationStep = int16(BW / num_steps);  % frequency deviation step in Hz.  This is fDEV, in Hz.  Can be positive or negative
+bf.RampMode = "disabled";
+bf.DelayStartWord = 4095;
+bf.DelayClockSource = "PFD";
+bf.DelayStartEnable = false;  % delay start
+bf.RampDelayEnable = false;  % delay between ramps.
+bf.TriggerDelayEnable = false;  % triangle delay
+bf.SingleFullTriangleEnable = false;  % full triangle enable/disable -- this is used with the single_ramp_burst mode
+bf.TriggerEnable = false;  % start a ramp with TXdata
+%% Flatten phaser phase/gain
+bf.RxGain(:) = 127;
+bf.RxAttn(:) = 0;
+bf.RxPhase(:) = 0;
+bf.RxLNAEnable(:) = true;
+bf.LatchRxSettings();
diff --git a/Matlab/setupPluto.m b/Matlab/setupPluto.m
new file mode 100644
index 0000000..d90a108
--- /dev/null
+++ b/Matlab/setupPluto.m
@@ -0,0 +1,19 @@
+function rx = setupPluto(plutoURI)
+rx = adi.AD9361.Rx('uri', plutoURI);
+rx.EnabledChannels = [1,2];
+rx.SamplesPerFrame = 1024;
+rx.CenterFrequency = 2e9;
+rx.kernelBuffersCount = 2; % Minimize delay in receive data
+rx.GainControlModeChannel0 = 'manual';
+rx.GainControlModeChannel1 = 'manual';
+rx.GainChannel0 = 6;
+rx.GainChannel1 = 6;
+rx.SamplingRate = 30e6;
+rx();
+
+tx = adi.AD9361.Tx('uri', plutoURI);
+tx.EnabledChannels = [1,2];
+tx.CenterFrequency = rx.CenterFrequency+2e6;
+tx.AttenuationChannel0 = -89;
+tx.AttenuationChannel1 = -89;
+end
\ No newline at end of file