diff --git a/BPMCalculate.m b/BPMCalculate.m
new file mode 100644
index 0000000..52d44f6
--- /dev/null
+++ b/BPMCalculate.m
@@ -0,0 +1,95 @@
+function [BPM1, BPM2, BPM3] = BPMCalculate(data, fs, sSize, threshold, setting)
+%BPMCalculate reads in a song in a .wav file and calculates its BPM 
+%   Detailed explanation goes here
+
+%% Process Data
+if setting == 1
+    ampValues = ampData(data,fs,sSize);
+
+else
+    for i = 1:(length(data) / sSize)
+        % Much of the antics here are due to 1 indexing...
+        sample = data([(((i - 1) * sSize) + 1):(i * sSize)], 1);
+        sample = abs(sample);
+        ampValues(:,i) = sample;
+    end
+end
+
+sample2Time = sSize / fs;
+
+%% Find Peaks
+sumVec = sum(ampValues, 1);
+
+sumChangeVec = zeros(length(sumVec),1);
+for i = 2:length(sumVec)
+    sumChangeVec(i) = sumVec(i) - sumVec(i - 1);
+end
+
+peaks = zeros(length(sumChangeVec));
+for i = 1:length(sumChangeVec)
+    if sumChangeVec(i) > threshold
+        peaks(i) = 1;
+    end
+end
+
+peakIndices = find(peaks);
+if isempty(peakIndices) || length(peakIndices) == 1 || length(peakIndices) == 2
+    BPM1 = nan;
+    BPM2 = nan;
+    BPM3 = nan;
+    return;
+end
+
+extraDistance = 0;
+nextIndex = 1;
+for i = 2:length(peakIndices)
+    % Special case where distance between is 1
+    if peakIndices(i) - peakIndices(i-1) == 1
+        extraDistance = extraDistance + 1;
+        continue
+    end
+    
+    % Special case where extra distance needs to be added on
+    if extraDistance ~= 0
+        % Add on extra distance from previous 1s
+        peakIndices(i) = peakIndices(i) + extraDistance;
+        peakSampleDist(nextIndex) = peakIndices(i) - peakIndices(i-1);
+
+        % Remove that extra distance so it doesn't affect future points
+        peakIndices(i) = peakIndices(i) - extraDistance;
+        extraDistance = 0;
+        nextIndex = nextIndex + 1;
+        continue
+    end
+
+    peakSampleDist(nextIndex) = peakIndices(i) - peakIndices(i-1);
+    nextIndex = nextIndex + 1;
+end
+
+%% Convert peakSampleDist to peakBPMs
+peakTimeDist = peakSampleDist .* sample2Time;
+
+peakBPMs = zeros(length(peakTimeDist),1);
+for i = 1:length(peakTimeDist)
+    peakBPMs(i) = (peakTimeDist(i) ^ (-1)) * 60;
+end
+
+% Convert BPMs to BPMs between 75 and 150
+i = 1;
+while i < length(peakBPMs) + 1
+    if peakBPMs(i) > 150
+        peakBPMs(i) = peakBPMs(i) / 2;
+
+    elseif peakBPMs(i) < 75
+        peakBPMs(i) = peakBPMs(i) * 2;
+
+    else
+        i = i + 1;
+    end
+end
+
+BPM1 = mean(peakBPMs);
+BPM2 = median(peakBPMs);
+BPM3 = mode(peakBPMs);
+
+end
\ No newline at end of file
diff --git a/ampData.m b/ampData.m
new file mode 100644
index 0000000..6bea50c
--- /dev/null
+++ b/ampData.m
@@ -0,0 +1,14 @@
+function [ampData] = ampData(data, fs, sSize)
+%ampData takes in .wav or .m4a file data and outputs amplitude data for samples
+%with size sSize
+%   Detailed explanation goes here
+
+%% Make a fft for each time sample and store in an array
+for i = 1:(length(data) / sSize)
+   % Much of the antics here are due to 1 indexing...
+   sample = data([(((i - 1) * sSize) + 1):(i * sSize)], 1);
+   [~, singleSided] = fftProcess(sample, length(sample), fs);
+   ampData(:,i) = singleSided;
+end
+
+end
\ No newline at end of file
diff --git a/bpmMain.m b/bpmMain.m
new file mode 100644
index 0000000..4c4cd7c
--- /dev/null
+++ b/bpmMain.m
@@ -0,0 +1,93 @@
+clear; close all;
+% Optimal Setting: Threshold: 0.0920 Resolution: 960
+
+% %% Resolution Test 1
+% [data, fs] = audioread("102ByAndBy.m4a");
+% actBPM = 102;
+% 
+% resolutions = [1080, 1200, 1320, 1440, 1560, 1680, 1800, 2040, 2160, 2280, 2400];
+% for j = 1:length(resolutions)
+%
+% Setting 1
+%    ampValues = ampData(data,fs,sSize);
+% 
+% Setting 2
+%     ampValues = zeros(resolutions(j), floor(length(data) / resolutions(j)));
+%     for i = 1:(length(data) / resolutions(j))
+%         % Much of the antics here are due to 1 indexing...
+%         sample = data([(((i - 1) * resolutions(j)) + 1):(i * resolutions(j))], 1);
+%         sample = abs(sample);
+%         ampValues(:,i) = sample;
+%     end
+% 
+%     % Data Processing
+%     sumVec = sum(ampValues, 1);
+% 
+%     sumChangeVec = zeros(length(sumVec),1);
+%     for i = 2:length(sumVec)
+%         sumChangeVec(i) = sumVec(i) - sumVec(i - 1);
+%     end
+% 
+%     sample2Time = resolutions(j) / fs;
+%     sampleTimeVec = [1:size(ampValues, 2)] .* sample2Time;
+% 
+%     % This is solely for testing and not part of the actual algorithm
+%     sampleBPMVec = sampleTimeVec .* actBPM ./ 60;
+% 
+%     figure();
+%     plot(sampleBPMVec, sumChangeVec);
+%     axis([2.5, 12, 0, 20]);
+% end
+
+%% Resolution Test 2
+resolutions = [1080, 1200, 1320, 1440, 1560, 1680, 1800, 2040, 2160, 2280, 2400];
+songList = ["75LeanOnMe", "77Power", "78NoWoman", "82Hopeless", "86Blinding", ...
+    "98TheMan", "102ByAndBy", "104InMyLife", "107Snow", "114Rainbow", ...
+    "116IceIceBaby", "120ProveIt", "121Talk", "128September", "137LetItBe", ...
+    "140RadioNW", "142AllLights"];
+numSongs = length(songList);
+for i = 1:numSongs
+    songList(i) = songList(i) + ".m4a";
+end
+
+% Threshold Vector
+starts = [3, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4];
+ends = [8, 8, 8, 12, 10, 10, 12, 15, 12, 12, 12];
+numThreshValues = 11;
+threshold = zeros(numThreshValues, length(resolutions));
+for j = 1:numThreshValues
+    for k = 1:length(resolutions)
+        threshold(j, k) = starts(k) + ((ends(k) - starts(k)) * (j - 1) / (numThreshValues - 1));
+    end
+end
+
+meanExpBPM = zeros(numSongs,numThreshValues,length(resolutions));
+medExpBPM = meanExpBPM;
+modeExpBPM = medExpBPM;
+
+% Actual BPM Vector
+actBPM = zeros(numSongs, numThreshValues, length(resolutions));
+actBPMplane = zeros(numSongs,numThreshValues);
+for j = 1:numThreshValues
+    actBPMplane(:,j) = [75, 77, 78, 82, 86, 98, 102, 104, 107, 114, 116, 120, 121, 128, 137, 140, 142];
+end
+
+for j = 1:length(resolutions)
+    actBPM(:,:,j) = actBPMplane;
+end
+
+for i = 1:numSongs
+    [data, fs] = audioread(songList(i));
+    %data = data(1:2000000);
+    for j = 1:numThreshValues
+        for k = 1:length(resolutions)
+            [meanExpBPM(i,j,k), medExpBPM(i,j,k), modeExpBPM(i,j,k)] = ...
+                BPMCalculate(data, fs, resolutions(k), threshold(j, k), 1);
+        end
+    end
+end
+
+%% Data Analysis
+[meanError, medError, modeError, meanBigError, medBigError, modeBigError, ...
+    settingTest, averageError, difficulty] = ...
+    dataAnalysis(meanExpBPM, medExpBPM, modeExpBPM, actBPM);
diff --git a/dataAnalysis.m b/dataAnalysis.m
new file mode 100644
index 0000000..4a214ee
--- /dev/null
+++ b/dataAnalysis.m
@@ -0,0 +1,171 @@
+function [meanError, medError, modeError, meanBigError, medBigError, modeBigError, settingTest, averageError, difficulty] = dataAnalysis(meanExpBPM, medExpBPM, modeExpBPM, actBPM)
+%dataAnalysis Takes experimental and theoretical BPM data for various
+%settings and returns analysis on data
+%   Detailed explanation goes here
+%% Calculate Error
+% Row: song
+% Column: Threshold
+% Plane: Resolution
+
+meanDifferences = meanExpBPM - actBPM;
+medDifferences = medExpBPM - actBPM;
+modeDifferences = modeExpBPM - actBPM;
+
+meanError = 100 .* abs(meanDifferences) ./ actBPM;
+medError = 100 .* abs(medDifferences) ./ actBPM;
+modeError = 100 .* abs(modeDifferences) ./ actBPM;
+
+% Adjust Error that is for values twice too big
+for i = 1:size(meanError, 1)
+    for j = 1:size(meanError, 2)
+        for k = 1:size(meanError, 3)
+            if meanError(i,j,k) > 80
+                meanExpBPM(i,j,k) = 0.5 * meanExpBPM(i,j,k);
+                meanDifferences(i,j,k) = meanExpBPM(i,j,k) - actBPM(i,j,k);
+                meanError(i,j,k) = 100 .* abs(meanDifferences(i,j,k)) ./ actBPM(i,j,k);
+            end
+
+            if medError(i,j,k) > 80
+                medExpBPM(i,j,k) = 0.5 * medExpBPM(i,j,k);
+                medDifferences(i,j,k) = medExpBPM(i,j,k) - actBPM(i,j,k);
+                medError(i,j,k) = 100 .* abs(medDifferences(i,j,k)) ./ actBPM(i,j,k);
+            end
+
+            if modeError(i,j,k) > 80
+                modeExpBPM(i,j,k) = 0.5 * modeExpBPM(i,j,k);
+                modeDifferences(i,j,k) = modeExpBPM(i,j,k) - actBPM(i,j,k);
+                modeError(i,j,k) = 100 .* abs(modeDifferences(i,j,k)) ./ actBPM(i,j,k);
+            end
+        end
+    end
+end
+
+%% Find Where the Function was innaccurate or Failed
+% Row: song
+% Column: Threshold
+% Plane: Resolution
+
+meanBigError = zeros(size(meanError, 1), size(meanError, 2), size(meanError, 3));
+medBigError = meanBigError;
+modeBigError = meanBigError;
+bigErrorPercent = 5;
+
+for i = 1:size(meanError, 1)
+    for j = 1:size(meanError, 2)
+        for k = 1:size(meanError, 3)
+            if meanError(i, j, k) > bigErrorPercent || isnan(meanError(i, j, k))
+                meanBigError(i, j, k) = 1;
+            end
+            if medError(i, j, k) > bigErrorPercent || isnan(medError(i, j, k))
+                medBigError(i, j, k) = 1;
+            end
+            if modeError(i, j, k) > bigErrorPercent || isnan(modeError(i, j, k))
+                modeBigError(i, j, k) = 1;
+            end
+        end
+    end
+end
+
+%% Count the NaN Values for Each Resolution & Threshold
+% Row: Threshold
+% Column: Resolution
+% Plane 1: Mean
+% Plane 2: Median
+% Plane 3: Mode
+
+nanCounts = zeros(size(meanError, 2), size(meanError, 3), 3);
+for j = 1:size(meanError, 2)
+    for k = 1:size(meanError, 3)
+        nanCount1 = 0;
+        nanCount2 = 0;
+        nanCount3 = 0;
+
+        for i = 1:size(meanError, 1)
+            if isnan(meanError(i, j, k))
+                nanCount1 = nanCount1 + 1;
+            end
+            if isnan(medError(i, j, k))
+                nanCount2 = nanCount2 + 1;
+            end
+            if isnan(modeError(i, j, k))
+                nanCount3 = nanCount3 + 1;
+            end
+        end
+        nanCounts(j, k, 1) = nanCount1;
+        nanCounts(j, k, 2) = nanCount2;
+        nanCounts(j, k, 3) = nanCount3;
+    end
+end
+
+%% Find Percent of Time that Song Failed
+% Row: Song
+% Column 1: Mean Errors
+% Column 2: Median Errors
+% Column 3: Mode Errors
+
+difficulty = zeros(size(meanError, 1), 3);
+for i = 1:size(meanError, 1)
+    difficulty(i, 1) = 100 * sum(sum(meanBigError(i, :, :))) / (size(meanError, 2) * size(meanError, 3));
+    difficulty(i, 2) = 100 * sum(sum(medBigError(i, :, :))) / (size(meanError, 2) * size(meanError, 3));
+    difficulty(i, 3) = 100 * sum(sum(modeBigError(i, :, :))) / (size(meanError, 2) * size(meanError, 3));
+end
+
+%% Find Average Error and Test Each Setting
+% Row: Threshold
+% Column: Resolution
+% Plane 1: Mean
+% Plane 2: Median
+% Plane 3: Mode
+
+settingTest = zeros(size(meanError, 2), size(meanError, 3), 3);
+averageError = settingTest;
+for i = 1:size(meanError, 2)
+    for j = 1:size(meanError, 3)
+        settingTest(i, j, 1) = sum(meanBigError(:, i, j));
+        settingTest(i, j, 2) = sum(medBigError(:, i, j));
+        settingTest(i, j, 3) = sum(modeBigError(:, i, j));
+        averageError(i, j, 1) = sum(meanError(:, i, j), 'omitnan') / (size(meanError, 1) - nanCounts(i, j, 1));
+        averageError(i, j, 2) = sum(medError(:, i, j), 'omitnan') / (size(medError, 1) - nanCounts(i, j, 2));            
+        averageError(i, j, 3) = sum(modeError(:, i, j), 'omitnan') / (size(modeError, 1) - nanCounts(i, j, 3));
+    end
+end
+
+% %% Find the Optimal Setting for each Song
+% % Row: Song
+% % Column 1: Threshold
+% % Column 2: Resolution
+% % Column 3: Percent Error
+% % Plane 1: Mean Error
+% % Plane 2: Median Error
+% % Plane 3: Mode Error
+% 
+% bestSetting = zeros(size(meanError, 1), 3, 3);
+% for i = 1:size(meanError, 1)
+%     bestError = 100;
+%     for j = 1:size(meanError, 2)
+%         for k = 1:size(meanError, 3)
+%             if meanError(i, j, k) < bestError
+%                 bestSetting(i, 1, 1) = j;
+%                 bestSetting(i, 2, 1) = k;
+%                 bestSetting(i, 3, 1) = meanError(i, j, k);
+%                 bestError = meanError(i, j, k);
+%             end
+% 
+%             if medError(i, j, k) < bestError
+%                 bestSetting(i, 1, 2) = j;
+%                 bestSetting(i, 2, 2) = k;
+%                 bestSetting(i, 3, 2) = medError(i, j, k);
+%                 bestError = medError(i, j, k);
+%             end
+% 
+%             if modeError(i, j, k) < bestError
+%                 bestSetting(i, 1, 3) = j;
+%                 bestSetting(i, 2, 3) = k;
+%                 bestSetting(i, 3, 3) = modeError(i, j, k);
+%                 bestError = modeError(i, j, k);
+%             end
+%         end
+%     end
+% end
+% 
+end
\ No newline at end of file
diff --git a/fftProcess.m b/fftProcess.m
new file mode 100644
index 0000000..dcc58bb
--- /dev/null
+++ b/fftProcess.m
@@ -0,0 +1,9 @@
+function [freqVec, singleSided] = fftProcess(data, recordLength, fs)
+%fftProcess Takes in sound data and returns fast fourier transformed data
+dataFFT = fft(data);
+twoSided = abs(dataFFT ./ recordLength); % Two sided spectrum
+singleSided = twoSided(1:floor(recordLength/2 + 1)); % Remove negative frequency peaks
+singleSided(2:(end - 1)) = 2 .* singleSided(2:(end - 1));
+
+freqVec = fs .* (0:recordLength/2) / recordLength;
+end
\ No newline at end of file