.gitignore

-# Object files
-*.o
-*.ko
-*.obj
-*.elf
-
-# Precompiled Headers
-*.gch
-*.pch
-
-# Libraries
-*.lib
-*.a
-*.la
-*.lo
-
-# Shared objects (inc. Windows DLLs)
-*.dll
-*.so
-*.so.*
-*.dylib
-
-# Executables
-*.exe
-*.out
-*.app
-*.i*86
-*.x86_64
-*.hex
-
-# Debug files
-*.dSYM/
-*.su
-
-
-# vrpn/build files
-src/vrpn/build*
-src/vrpn/pc_linux64/*
-
-#Exacutables
-./BlueTooth
+test.log
\ No newline at end of file

.gitmodules

+[submodule "groundStation/src/vrpn"]
+	path = groundStation/src/vrpn
+	url = https://github.com/vrpn/vrpn

Makefile

-# Declaration of variables
-GCC=gcc
-GXX=g++
-CFLAGS= -Wall  -std=c99 -g
-CPPFLAGS= -Wall  -std=c++11 -g
-INCLUDES = $(foreach dir, $(INCDIR), -I$(dir))
-
-# Directories
-SRCDIR=src 
-INCDIR=inc src src/vrpn src/vrpn/quat src/vrpn/build 
-OBJDIR=obj
-
-# Final exacutable name
-EXE=BlueTooth
-
-# File names
-CSOURCES := $(foreach dir, $(SRCDIR), $(wildcard $(dir)/*.c )) 
-CPPSOURCES := $(foreach dir, $(SRCDIR), $(wildcard $(dir)/*.cpp ))
-COBJECTS = $(CSOURCES:.c=.o) 
-CPPOBJECTS = $(CPPSOURCES:.cpp=.o)
-OBJECTS = $(CPPOBJECTS) $(COBJECTS) 
-
-LIBS= -lpthread -lbluetooth -lvrpn -lquat -Lsrc/vrpn/build -Lsrc/vrpn/build/quat 
-
-# Default target
-all: logs vrpn/build $(EXE)
-
-# Main target
-$(EXE): $(OBJECTS) 
-	$(GXX) $(CPPFLAGS) $^ -o $@  $(INCLUDES) $(LIBS)
-
-$(COBJECTS) : %.o : %.c
-	$(GCC)  $(CFLAGS) -c $< -o $@ $(INCLUDES) $(LIBS)
-
-$(CPPOBJECTS) : %.o : %.cpp
-	$(GCC)  $(CPPFLAGS) -c $< -o $@ $(INCLUDES) $(LIBS)
-
-vrpn/build:
-	mkdir -p src/vrpn/build
-	cd src/vrpn/build && cmake .. && make
-
-logs:
-	mkdir -p logs
-
-clean_logs:
-	rm -f logs/*
-
-clean:
-	rm -f $(EXE) $(OBJECTS)
-
-debug:
-	@echo $(COBJECTS)
-	@echo $(CPPOBJECTS)

README.md

-# groundStation
\ No newline at end of file

Readme.md

+Repository for 2016-2017 MicroCART project.

controls/MATLAB/IMU/accelerometer_noise_analysis.m

+close all; clear all; clc; 
+
+filename = 'imu.csv';
+data = load_IMU_data_from_csv( filename ) ; 
+
+t = data.t; 
+ax = data.ax; 
+ay = data.ay; 
+az = data.az; 
+gx = data.gx; 
+gy = data.gy; 
+gz = data.gz; 
+
+% analyze the noise (deviation from mean) on each axis and compare the
+% emperical distribution to one generated from a quantized Gaussian
+% distrubtuion with the same variance. 
+quantized_noise_gauss_anal( ax , 'Accelerometer' , 'x-axis' , 'gs' )
+quantized_noise_gauss_anal( ay , 'Accelerometer' , 'y-axis', 'gs' )
+quantized_noise_gauss_anal( az , 'Accelerometer' , 'z-axis', 'gs' )
+
+
+quantized_noise_gauss_anal( gx , 'Gyroscope' , 'x-axis' , 'rad/s')
+quantized_noise_gauss_anal( gy , 'Gyroscope' , 'y-axis' , 'rad/s' )
+quantized_noise_gauss_anal( gz , 'Gyroscope' , 'z-axis' , 'rad/s' )
\ No newline at end of file

controls/MATLAB/IMU/load_IMU_data_from_csv.m

+function data = load_IMU_data_from_csv( filename )
+% load_IMU_data_from_csv load data from csv file into workspace 
+%
+%SYNTAX
+% data = load_IMU_data_from_csv( filename )
+% 
+%DESCRIPTION
+% ...to be written...
+%
+%
+%Inputs:
+% filename - string specifying the file 
+%
+%Options:
+% none 
+%
+%Outputs:
+% data - IMU data loaded into a structure (specify more later)  
+%
+%EXAMPLES
+%
+%NOTES
+% This function currently assumes the csv format that was used as of the 
+% date: 11/15/2016 (m/d/y) 
+%
+%AUTHORS
+% Matt Rich - m87rich@iastate.edu
+%
+%DEVELOPERS
+%
+
+%
+%DEVELOPMENT NOTES
+%
+% dates are in m-d-y format 
+%
+% Initial bare bones function. 
+% - Matt Rich 11-15-2016
+%
+delimiter = ',';
+formatSpec = '%f%f%f%f%f%f%f%[^\n\r]';
+fileID = fopen(filename,'r');
+dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'EmptyValue' ,NaN, 'ReturnOnError', false);
+fclose(fileID);
+data.t = dataArray{:, 1};
+data.ax = dataArray{:, 2};
+data.ay = dataArray{:, 3};
+data.az = dataArray{:, 4};
+data.gx = dataArray{:, 5};
+data.gy = dataArray{:, 6};
+data.gz = dataArray{:, 7};
+
+
+end
+

controls/MATLAB/IMU/microcart_imu_analysis_1.m

+% function microcart_imu_analysis_1
+
+close all;  clc; 
+filename = 'imu.csv';
+delimiter = ',';
+formatSpec = '%f%f%f%f%f%f%f%[^\n\r]';
+fileID = fopen(filename,'r');
+dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'EmptyValue' ,NaN, 'ReturnOnError', false);
+fclose(fileID);
+t = dataArray{:, 1};
+ax = dataArray{:, 2};
+ay = dataArray{:, 3};
+az = dataArray{:, 4};
+gx = dataArray{:, 5};
+gy = dataArray{:, 6};
+gz = dataArray{:, 7};
+clearvars filename delimiter formatSpec fileID dataArray ans;
+
+% magnitude of acceleration in gs 
+a = ( ax.^2 + ay.^2 + az.^2 ).^0.5 ; 
+
+% magnitude of deviation from ideal 
+ea1 = 1 - a ;
+
+% magnitude of deviation from mean 
+ea2 = a - mean(a) ; 
+
+% RMS error 
+ea1rms = sqrt(1/length(ea1)*(ea1')*ea1) ; 
+ea2rms = sqrt(1/length(ea2)*(ea2')*ea2) ; 
+
+subplot(3,1,1); 
+plot(t, a , t , ones(length(t),1)); axis([min(t),max(t),min([a;1])-0.01,max([a;1])+0.01]); grid;
+xlabel('Time (s)'); ylabel('Acceleration (gs)'); 
+legend('Measured','Ideal'); 
+
+
+subplot(3,1,2); 
+plot(t, ea1 ); axis([min(t),max(t),min(ea1)-0.01,max(ea1)+0.01]); grid;
+xlabel('Time (s)'); ylabel('Error From Ideal (gs)'); 
+
+subplot(3,1,3); 
+plot(t, ea2 ); axis([min(t),max(t),min(ea2)-0.01,max(ea2)+0.01]); grid;
+xlabel('Time (s)'); ylabel('Error From Mean (gs)'); 
+
+% 
+figure; 
+X = fft(ea2); 
+N = length(ea2); 
+fs = 1/(t(2)-t(1)); 
+F = ([-N/2:N/2-1]/N)*fs; %frequency axis for plotting FFTs 
+% subplot(2,1,1); 
+plot(F,fftshift(abs(X))); grid; 
+xlabel('frequency (Hz)');
+ylabel(['|FFT(x)|']);
+% subplot(2,1,2); 
+% plot(F,20*log10(fftshift(abs(X)))); grid; 
+% xlabel('frequency (Hz)');
+% ylabel('|FFT(x)| dB');
+suptitle('x = Error From Mean (gs)'); 
+
+
+eax = ax - mean(ax); 
+
+vx = var(ax); 
+sigmax = sqrt(vx);  
+
+figure; 
+
+bins = unique(eax) ; 
+dbins = bins(2)-bins(1); % the quantization levels consistent 
+
+bins_right = bins + dbins/2; 
+bins_left = bins - dbins/2; 
+cgauss_right = normcdf(bins_right,0,sigmax); 
+cgauss_left = normcdf(bins_left,0,sigmax); 
+
+empf = length(eax)*(cgauss_right-cgauss_left); 
+stairs(bins_left,empf,'r','LineWidth',1) ; 
+xlabel('Error From Mean (gs)'); 
+ylabel('Occurences'); 
+
+grid; 
+hold on; 
+h = hist(eax,bins); 
+stem(bins,h); 
+
+legend(['Quantized Gaussian Noise: N(0,',num2str(vx),')'],'Emperical Distribution')
+
+
+% quantized_accel_noise_gauss_anal( ax )
+
+
+
+% end
+
+
+

controls/MATLAB/IMU/quantized_noise_gauss_anal.m

+function varargout = quantized_noise_gauss_anal( a , sensor_name, axis , units )
+% quantized_noise_gauss_anal analysis for quantized sensor data 
+%
+%SYNTAX
+% quantized_noise_gauss_anal(a , sensor_name, axis , units )
+% 
+%DESCRIPTION
+% ...to be written...
+%
+%
+%Inputs:
+%           a - vector of a single axis accelerometer reading 
+% sensor_name - string specifying the axis label 
+%        axis - string specifying the axis label 
+%       units - string specifying the physical units on that axis 
+%
+%Options:
+% none 
+%
+%Outputs:
+% none 
+%
+%EXAMPLES
+%
+%NOTES
+% This function is only really meaningful for a static (not moving) test. 
+%
+%AUTHORS
+% Matt Rich - m87rich@iastate.edu
+%
+%DEVELOPERS
+%
+
+%
+%DEVELOPMENT NOTES
+%
+% dates are in m-d-y format 
+%
+% Initial bare bones function. 
+% - Matt Rich 11-15-2016
+%
+% Change funcdtion name and code to not specify what sensor it is for so it
+% can be used for any sensor producing quantized output
+% - Matt Rich 11-16-2016
+%
+
+mu = mean(a); 
+v = var(a); 
+sigma = sqrt(v); 
+
+ea = a - mu; 
+
+mue = mean(ea); %calculate the mean of the error from mean
+
+
+figure; 
+
+bins = unique(ea) ; 
+dbins = bins(2)-bins(1); % assuming the quantization levels consistent 
+
+bins_right = bins + dbins/2; 
+bins_left = bins - dbins/2; 
+cgauss_right = normcdf(bins_right,mue,sigma); 
+cgauss_left = normcdf(bins_left,mue,sigma); 
+
+emp_dist = length(ea)*(cgauss_right-cgauss_left); 
+stairs(bins_left,emp_dist,'r','LineWidth',1) ; 
+xlabel(['Error From Mean (',units,')']); 
+ylabel('Occurences'); 
+grid; 
+hold on; 
+h = hist(ea,bins); 
+stem(bins,h); 
+legend(['Quantized Gaussian Noise: N(',num2str(mue),',',num2str(v),')'],'Emperical Distribution'); 
+title([sensor_name,' deviation from mean on ',axis]); 
+
+
+end
+

controls/MATLAB/momentOfInertia/getAverageMomentOfInertia.m

+function averageMOI = ...
+    getAverageMomentOfInertia( directoryLocation, massUsed_grams_array )
+    %Given a directory of data labeled "Test_1", "Test_2", etc, will find
+    %the moment of inertia for each data set, using the mass used in grams
+    %at the corresponding index of the massUsed_grams_array. And from
+    %finding all of the moments of inertia, in this way, this function will
+    %return the average.
+    
+    %First find the number of tests we need to run (will be equal to the
+    %length of the massUsed_grams_array)
+    numRuns = length(massUsed_grams_array);
+    
+    %Create a sum variable, to which we will add our calculated MOI values:
+    sumOfMOIs = 0;
+    
+    %For each of the runs, calculate the moment of inertia
+    for i = 1:numRuns
+    
+        %Create the file name for the current run, in string:
+        currFileName = ['Test_', num2str(i), '.csv'];
+        
+        %Create the full file path for the current run:
+        currFullFilePath = [directoryLocation, '\', currFileName];
+        
+        %Get the current mass used in grams
+        massUsed_grams = massUsed_grams_array(i);
+        
+        %Get the moment of inertia for this run
+        curr_MOI = getMomentOfInertia(currFullFilePath, massUsed_grams);
+        
+        %Add this run's moment of inertia value calculated to the sum
+        sumOfMOIs = sumOfMOIs + curr_MOI;
+        
+    end %for i = 1:numRuns
+    
+    %Get the average moment of inertia value, by dividing by the total
+    %number of runs
+    averageMOI = sumOfMOIs / numRuns;
+
+end

controls/MATLAB/momentOfInertia/getMomentOfInertia.m

+function momentOfInertia = getMomentOfInertia( fullFilePath, massUsed_grams )
+%This function will calculate a single moment of inertia value, given a
+%file location to a csv file, along with the mass value used on the mass
+%handle as a torque. (Note, there's no need to have to consider the mass of
+%the handle, which is 5.0 grams, since this function takes care of that)
+
+    %fileLocation = 'C:\Users\Tara\Desktop\Project Documents\Current Project Documents\EE 491\Data\Physics Department Measurements\Test\';
+    %fileName = 'Test_1.csv';
+
+    %Mass used in grams
+    %mass_grams = 200.1;
+    massHandle_grams = 5.0;
+    totalMass_grams = massUsed_grams + massHandle_grams;
+    totalMass_kg = totalMass_grams / 1000.0;
+
+    %Define a cut-off start angle in radians
+    angleStart_cutoff = 0.1;
+
+    %Radius of the knob where the string was pulling from (in meters)
+    radius = 0.04445;   %(4.445 centimeters - or 1.75 inches)
+
+    %Acceleration constant due to gravity (in meters per second squared)
+    g = 9.80665;
+
+    %Define a cut-off "padding" to stop taking data before reaching the end of
+    %the data recorded in the file
+    %For example, if dataEndPadding is 10, we will "cut-off" our data when we
+    %get to 10 data points from the very end of the data run
+    dataEndPadding = 5;
+
+    %Define columns of the table
+    timeCol = 1;
+    angleCol = 2;
+    angularVelCol = 3;
+    %angularAccelCol = 4;
+
+    %Create full file path
+    %fullFilePath = [fileLocation, fileName];
+
+    %Bring in Data
+    dataTable = readtable(fullFilePath);
+
+    %Extract the time, angle (in radians), and angular velocity (in radians per
+    %second) arrays from the data table
+    timeArray = dataTable{:, timeCol};
+    angleArray = dataTable{:, angleCol};
+    angularVelArray = dataTable{:, angularVelCol};
+
+    %Find the first positive angle value
+    i_start = find(angleArray > angleStart_cutoff, 1);
+
+    %Get the start time and start angular position, based on this start index
+    time_start = timeArray(i_start);
+    angle_start = angleArray(i_start);
+
+    %Find the initial angular velocity, which is the corresponding angular 
+    %velocity value at this start index (in radians per second)
+    w_o = angularVelArray(i_start);
+
+    %Determine the end time of our data collection, so first get the length of
+    %the time array
+    timeArray_len = length(timeArray);
+
+    %Get the last index based on the length of the data and the "padding"
+    %amount defined at the beginning
+    i_end = timeArray_len - dataEndPadding;
+
+    %Get the end time and end angular position in radians
+    time_end = timeArray(i_end);
+    angle_end = angleArray(i_end);
+
+    %Get the difference in time and angular position (in radians)
+    angle_diff = angle_end - angle_start;
+    time_diff = time_end - time_start;
+
+    %Using this expression, where alpha is the constant angular acceleration:
+    %   theta_diff = w_o*t_diff + (1/2)*alpha*(t_diff^2)
+    %We can solve for alpha to get
+    %   (theta_diff - w_o*t_diff) = (1/2)*alpha*(t_diff^2)
+    %   2*(theta_diff - w_o*t_diff) / (t_diff^2) = alpha
+
+    %Thus the angular acceleration in radians per second-squared, denoted by 
+    %alpha can be found like this:
+    angularAccel = 2*(angle_diff - w_o*time_diff) / (time_diff * time_diff);
+
+    %This constant angular acceleration should be due to the torque caused by
+    %the gravitational force caused by our falling mass (m*g*r)
+    torque_gravity = totalMass_kg * g * radius;
+
+    %The moment of inertia would be this torque divided by the constant angular
+    %acceleration value (in kilograms meters-squared)
+    momentOfInertia = torque_gravity / angularAccel;
+
+
+end
+

controls/MATLAB/momentOfInertia/getMomentOfInertiaFinalData.m

+%Create the massesUsed_grams_array, which represents the masses we used for
+%each run, which we keep the same for the runs done for each of the
+%testing (yaw, pitch, roll, and their calibration tests as well)
+massesUsed_grams = [200.1, 300.1, 100.1];
+numTrialsPerMass = 5;
+massUsed_grams_array = [ ones(1, numTrialsPerMass) * massesUsed_grams(1), ...
+                         ones(1, numTrialsPerMass) * massesUsed_grams(2), ...
+                         ones(1, numTrialsPerMass) * massesUsed_grams(3)];
+
+%Set the directory locations for the pitch, roll, and yaw data and their
+%calibration runs
+topLevelDir = 'C:\Users\Tara\Desktop\Project Documents\Current Project Documents\EE 491\Data\Physics Department Measurements\Moment of Inertia Data';
+directoryLocation_yawCalib = [topLevelDir , '\Calibration Testing Yaw'];
+directoryLocation_pitchCalib = [topLevelDir, '\Calibration Testing Pitch and Roll'];
+directoryLocation_rollCalib = [topLevelDir, '\Calibration Testing Pitch and Roll'];
+directoryLocation_yawData = [topLevelDir, '\Yaw Data'];
+directoryLocation_pitchData = [topLevelDir, '\Pitch Data'];
+directoryLocation_rollData = [topLevelDir, '\Roll Data'];
+                     
+%Get the average moment of inertia values for the pitch, roll, and yaw data
+%as well as their calibration funs
+averageMOI_yawCalib = getAverageMomentOfInertia( ...
+    directoryLocation_yawCalib, massUsed_grams_array );
+
+averageMOI_pitchCalib = getAverageMomentOfInertia( ...
+    directoryLocation_pitchCalib, massUsed_grams_array );
+
+averageMOI_rollCalib = getAverageMomentOfInertia( ...
+    directoryLocation_rollCalib, massUsed_grams_array );
+
+averageMOI_yawData = getAverageMomentOfInertia( ...
+    directoryLocation_yawData, massUsed_grams_array );
+
+averageMOI_pitchData = getAverageMomentOfInertia( ...
+    directoryLocation_pitchData, massUsed_grams_array );
+
+averageMOI_rollData = getAverageMomentOfInertia( ...
+    directoryLocation_rollData, massUsed_grams_array );
+
+%Now determine the quadcopter's average moment of inertia about the pitch,
+%roll, and yaw axes of rotation
+quadYaw = averageMOI_yawData - averageMOI_yawCalib
+quadPitch = averageMOI_pitchData - averageMOI_pitchCalib
+quadRoll = averageMOI_rollData - averageMOI_rollCalib
\ No newline at end of file

controls/MATLAB/thrustAndDragConstant/calcDragConstant.m

+function [ Kd ] = calcDragConstant( data, Pmin, Pmax, Rm, Kv, Kq, If )
+%CALC_DRAG_CONSTANT 
+%   Calculates the drag constant (Kd) given experimental data. The drag
+%   constant is described in detail in sections 4.2.5, 4.2.6, and 5.5.4.1 
+%   of "Model development, system identification, and control of a 
+%   quadrotor helicopter" by Matt Rich.
+
+%   Input Arguments:
+%   data: experimental data
+%   Pmax: Calculated maximum duty cycle of the function generators PWM wave
+%   that the ESC can handle.
+%   Pmin: Calculated minimum duty cycle of the function generators PWM wave
+%   for initialization of the ESC.
+%   Rm: Motor resistance
+%   Kv: Back-EMF constant of the motor
+%   Kq: Torque constant of the motor
+%   If: No-Load (friction) current of the motor
+
+
+% Convert RPM to angular speed of each rotor.
+rotor_speed_0 = data.(2) * (pi/30);
+rotor_speed_1 = data.(3) * (pi/30);
+rotor_speed_2 = data.(4) * (pi/30);
+rotor_speed_3 = data.(5) * (pi/30);
+
+% Refer to the sections described in the header of this function
+% for a better understanding of what this loop is doing. Note that these
+% equations for each individual rotor. Calculating the total drag constant
+% requires taking into account all four motors. In this function this is
+% done by stacking the data of each motor in a single array and performing
+% a least squares approximation of all of the data.
+%
+%       u: Defined in section 4.2.5
+%       omega: vector of each rotors speed in rad/s
+%       A: column vector of each rotors speed squared
+%       b: Defined in section 5.5.4.1
+%       Kd_vector: Vector containing all experimental Kd values
+%       Vb: Battery voltage
+u1 = ((data.(1)/100) - Pmin)/(Pmax - Pmin);
+u2 = ((data.(1)/100) - Pmin)/(Pmax - Pmin);
+u3 = ((data.(1)/100) - Pmin)/(Pmax - Pmin);
+u4 = ((data.(1)/100) - Pmin)/(Pmax - Pmin);
+u = [u1; u2; u3; u4];
+Vb = [data.(6); data.(6); data.(6); data.(6)];
+omega = [rotor_speed_0; rotor_speed_1; rotor_speed_2; rotor_speed_3];
+A = omega.^2;
+b = ((u.*Vb)/(Rm*Kq))-omega./(Rm*Kq*Kv)-If/Kq;
+
+Kd = ((A'*A)^-1)*A'*b
+end
+

controls/MATLAB/thrustAndDragConstant/calcThrustConstant.m

+function [ Kt ] = calcThrustConstant( data )
+%CALC_THRUST_CONSTANT
+%   Calculates the thrust constant (Kt) given experimental data. The thrust
+%   constant is described in detail in sections 4.2.1.1 and 5.5.1 of 
+%   "Model development, system identification, and control of a quadrotor
+%   helicopter" by Matt Rich.
+
+%   Input Arguments:
+%   data: experimental data
+
+% Convert RPM to angular speed of each rotor.
+rotor_speed_0 = data.(2) * (pi/30);
+rotor_speed_1 = data.(3) * (pi/30);
+rotor_speed_2 = data.(4) * (pi/30);
+rotor_speed_3 = data.(5) * (pi/30);
+
+% Define the A matrix as the sum of each rotors speed squared.
+A = (rotor_speed_0.^2 + rotor_speed_1.^2 + rotor_speed_2.^2 + rotor_speed_3.^2);
+
+% Convert weight (g) to thrust force (N) by converting grams to kilograms and
+% multiplying by the acceleration of gravity.
+T = (data.Scale_g_/1000)*9.8;
+
+% Calculate the thrust constant (Kt) through the following
+% equation: Kt = (A'A)^-1.*A'.*T as defined on page 65 of "Model
+% development, system identification, and control of a quadrotor
+% helicopter" (Matt Rich's Thesis).
+Kt = ((A'*A)^(-1))*(A'*T) 
+
+end
+

controls/MATLAB/thrustAndDragConstant/thrustDragConstantCalculations.m

+% Add path to zeroLoadCurrent() function
+addpath('C:\Users\Andy\Documents\School\MicroCART\GitRepo\MicroCART_17-18\controls\MATLAB\zeroLoadCurrent')
+
+% Import data as a table.
+data = readtable('C:\Users\Andy\Documents\MATLAB\MicroCART\Thrust and Drag Constant\Thrust and Drag Constant Data.xlsx');
+
+filePath_noLoadCurrentData = 'C:\Users\Andy\Documents\MATLAB\MicroCART\Zero Load Current\No Load Friction Current.csv';
+
+[I_0, I_1, I_2] = zeroLoadCurrent(filePath_noLoadCurrentData);
+
+rotor_speed_0 = data.(2) * (pi/30);
+rotor_speed_1 = data.(3) * (pi/30);
+rotor_speed_2 = data.(4) * (pi/30);
+rotor_speed_3 = data.(5) * (pi/30);
+
+% Define Pmin and Pmax, as well as the motor parameters Rm, Kv, Kq, and If
+Pmin = 0.40;
+Pmax = 0.8;
+Rm = 0.2308;
+Kv = 96.3422;
+Kq = 96.3422;
+If0 = I_0 * sign(rotor_speed_0) + I_1 * rotor_speed_0 + I_2 * rotor_speed_0.^2;
+If1 = I_0 * sign(rotor_speed_1) + I_1 * rotor_speed_1 + I_2 * rotor_speed_1.^2;
+If2 = I_0 * sign(rotor_speed_2) + I_1 * rotor_speed_2 + I_2 * rotor_speed_2.^2;
+If3 = I_0 * sign(rotor_speed_3) + I_1 * rotor_speed_3 + I_2 * rotor_speed_3.^2;
+If = [If0; If1; If2; If3];
+
+% Call the calc_thrust_constant() function.
+Kt = calcThrustConstant(data);
+
+% Call the calc_drift_constant() function.
+Kd = calcDragConstant(data, Pmin, Pmax, Rm, Kv, Kq, If );
+

controls/MATLAB/zeroLoadCurrent/Kd_error.m

+% Read in data 
+data = readtable('C:\Users\Andy\Documents\MATLAB\MicroCART\Thrust and Drag Constant\Thrust and Drag Constant Data.xlsx');
+
+filePath_noLoadCurrentData = 'C:\Users\Andy\Documents\MATLAB\MicroCART\Zero Load Current\No Load Friction Current.csv';
+
+[I_0, I_1, I_2, dutyCycle_error, error, residual_error, residual_error_ConstantIf] = zeroLoadCurrent(filePath_noLoadCurrentData);
+
+Pmin = 0.40;
+Pmax = 0.8;
+Rm = 0.2308;
+Kv = 96.3422;
+Kq = 96.3422;
+
+% Convert RPM to angular speed of each rotor.
+rotor_speed_0 = data.(2) * (pi/30);
+rotor_speed_1 = data.(3) * (pi/30);
+rotor_speed_2 = data.(4) * (pi/30);
+rotor_speed_3 = data.(5) * (pi/30);
+
+duty_cycle_percentage = data.(1);
+Vb = data.(6);
+
+u = ((data.(1)/100) - Pmin)/(Pmax - Pmin);
+
+If = I_0 * sign(rotor_speed_0) + I_1 * rotor_speed_0 + I_2 * rotor_speed_0.^2;
+w_num = -1 + sqrt( 1 - 4*Rm*Kv*Kq*Kd*(Kv*Rm*If - Kv*u.*Vb));
+w_den = 2*Rm*Kv*Kq*Kd;
+w = w_num / w_den;
+
+figure()
+hold on
+plot(duty_cycle_percentage, w); grid()
+plot(duty_cycle_percentage, rotor_speed_0);
+%plot(duty_cycle_percentage, rotor_speed_1);
+%plot(duty_cycle_percentage, rotor_speed_2);
+%plot(duty_cycle_percentage, rotor_speed_3);
+
+xlabel('Duty Cycle Percentage');
+ylabel('Rotor Speed (rad/s)')
+
+residual_error
+
+residual_error_ConstantIf
\ No newline at end of file

controls/MATLAB/zeroLoadCurrent/zeroLoadCurrent.m

+function [I_0, I_1, I_2, dutyCycle, error, residual_error, residual_error_constantIf] = zeroLoadCurrent(absoluteFilePath)
+%This function takes in an absolute file path to a .csv or Excel file with 
+%the following format:
+%   Second column should be the motor speed in revolutions per minute
+%   Third column should be the current (A)
+%absoluteFilePath = 'C:\Users\Andy\Documents\MATLAB\MicroCART\Zero Load Current\No Load Friction Current.csv';
+
+    %Get the data from the file
+    data = readtable(absoluteFilePath);
+
+    %Extract the motor speed column and convert it to radians per second
+    dutyCycle = data.(1);
+    motorSpeed_rpm = data.(2);
+    motorSpeed_radPerSec = motorSpeed_rpm * pi/30;
+
+    %Extract the current column
+    current = data.(3);
+
+    %Define the three columns of the A matrix for the least squares 
+    %approximation, using this equation:
+    %   
+    %   I_f = sgn(motorSpeed)*I_0 + motorSpeed*I_1 + (motorSpeed^2)*I_2
+    %
+    %   for this least squares equation: A*x = b (where A and b are matrices)
+    %
+    % where: sgn(motorSpeed) -- first column of the A matrix,
+    %        motorSpeed      -- second column of the A matrix
+    %        motorSpeed^2    -- third column of the A matrix
+    %        I_f             -- b matrix
+    %
+    %  and: [I_0; I_1; I_2]  -- x vector that we are solving for
+    A_col1 = sign(motorSpeed_radPerSec);
+    A_col2 = motorSpeed_radPerSec;
+    A_col3 = motorSpeed_radPerSec .^ 2;
+
+    %Create the A matrix from its three columns
+    A = [A_col1, A_col2, A_col3];
+
+    %Define the b matrix
+    b = current;
+
+    %Least squares approximation -- solving for x, which is I_vector
+    I_vector = (A' * A)^(-1) * A' * b;
+
+    %Get the error vector
+    error = A*I_vector - b;
+    
+    If_constant = 0.511;
+    error_withIfConstant = A*[If_constant; 0; 0] - b;
+    
+    %Get the residual error
+    residual_error = sum(error.^2);
+    
+    residual_error_constantIf = sum(error_withIfConstant.^2);
+
+    %Extract the components of the vector
+    I_0 = I_vector(1);
+    I_1 = I_vector(2);
+    I_2 = I_vector(3);
+
+end %function