% I have used some functions which are at the bottom of this file. I placed
% comments in the function code as well, so in the 'main' code I will not
% explain how the function works, just what it returns.

% I will give comments AFTER every line of code.

% I used the matlab-given function 'arrayfun' a few times. What this
% function does is that it passes every point in a given array to a
% function, and places the result in it's respective place in a new array.
% It's syntax is "NewArray = arrayfun(@FunctionName, DataArray".

%% ****No for loops were used or harmed in the creation of this code.****
%% Code Below

Fs = 44100;
% Fs = The number of samples per second in which the input file is divided.
% See lecture notes for how/why this method is used for digitalizing sound.
% The industry standard is 44100.

[arr_A,Fs] = audioread('klaatu.wav');
% Read the audio file and store it's data in array form.

arr_B = arrayfun(@absolute,arr_A);
% We make the input data absolute (the root-square step of the RMS) so that
% when we take the 'integral' of the data it doesn't result in 0.

K = 100001;
arr_C = smoothdata(arr_B, 'movmean', K);
% We use the 'smoothdata' function supplied by matlab to find the average
% ('movmean' flag) 'loudness' over a certain timeperiod (in this case the
% timeperiod is K = 100 001 samples long). For every point in arr_B it
% looks X values back and X values forward, adds these values together and
% divides this by X+X+1 to get the average over that period, and places it
% in arr_C. This is the 'Mean' step of RMS.

% If K were to be too small, let's say K = 1, then later on when
% normalizing the sound, the result for every data point would be exactly
% the same, removing the variations needed for carrying the data (The
% result would probably be a continuous single tone, if it results in
% anything at all).

global MaxValRMS;
MaxValRMS = max(arr_C);
% We find the highest point of the RMS. We will use this to find the
% average amplification that every datapoint will strive towards. See the
% explanation at the amplitudeMod function for more info.

% The variable MaxValRms is declared global because matlab works in such a
% way that any variable in any scope is NOT accesable in any other scope.
% In this case, we're in the 'Base' scope, but any functions have their own
% 'Function' scope, and cannot see any variables in the 'Base' scope.
% Declaring the variable as global in both 'Base' and 'Function' allows for
% sharing a single variable accross multiple scopes.

arr_AmpMod = arrayfun(@amplitudeMod,arr_C);
% We make an array with values that are the INVERSE of the RMS. The loud
% parts will have a low value, and the quit parts will have a high value.

arr_Normalized = arr_A .* arr_AmpMod;
% We now multiply the original input signal with the inverse-RMS. Since the
% quit parts have a high inverse-RMS value, multiplying them with the input
% signal will increase the 'quit part' amplitude. Since the loud parts have
% a low inverse-RMS value, multiplying the input with this inverse RMS
% value results in a lower amplitude.
% Basically, this is the part where we make the quit parts loud and the
% loud parts quit, as to make them all approximately the same level.

global MaxValInput;
MaxValInput = max(arr_B);
arr_Output = arrayfun(@removePeaks,arr_Normalized);
% This part above is completely optional. After running the normalization
% part above, I saw that there were some data spikes caused by the fact
% that the RMS is the value over a certain period, but in this period could
% be very high peaks that now also get amplified the same as its quiter
% neighbour values. This results in an output with some spikes at some
% parts, which could be unenjoyable. This part just removes the peaks.

audiowrite('output.flac',arr_Output,Fs);
% Write the final output in a sound file to be able to listen to the
% result.

%% The part below plots all the graph of every intermediate step.

figure(1);
plot(arr_A);
% "Representation of Input Data"
figure(2);
plot(arr_B);
% "Absolute-ed input data"
hold on
plot(arr_C);
% RMS of absolute'd input array, overlayed on the absolute'd array to show
% the level of RMS at any point in time, and easily visually compare to the
% input data.
figure(3);
plot(arr_AmpMod);
% "Inverse of RMS"
figure(4);
plot(arr_Normalized);
% "Input data multiplied by the inverse amplitudes, normalized"
figure(5);
plot(arr_Output);
% "Output array, without spikes"

%% Functions
function output = absolute(x)
    output = sqrt(x^2);
end
% Takes the input value, squares it and takes the root of this new value.
% This is now the absolute value of the original value.

function output = amplitudeMod(x)
    global MaxValRMS;
    output = (MaxValRMS/2) / x;
end
% As explained above, MaxValRMS is global to be accesable in this scope.
% MaxValRms/2 is half the loudness of the loud parts, which should also be
% about double the loudness of quit parts. By dividing this by low (quit)
% value, it will return a high value (1 / 0.5 = 2), and by dividing the
% MaxValRMS/2 by a high value, it will return a lower value (1 / 2 = 0.5).

% Basically this returns the amplification needed to get a datapoint to
% the average loudness.

function output = removePeaks(x)
    global MaxValInput;
    if x > ((MaxValInput/3)*2)
        output = ((MaxValInput/3)*2);
    elseif x < ((-MaxValInput / 3) * 2)
        output = ((-MaxValInput/3)*2);
    else
        output = x;
    end
end
% This function checks if the data point is above 2/3rd of the original
% 'loud' part (meaning it is a spike), and if that is the case it changes
% its value to be 2/3rd of the original input loud part. This checks for
% both positive and negative datapoints since sound in this case is an
% osscilation around zero.