% 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.
K = 100001;
%Smoothness factor to be used later. This is the amount of data points
%around a data point over which the average is calculated.

[s,Fs] = audioread('klaatu.wav');
%Read audio file and place it in array s

Output_Array = thatwayRMSConverter(s,K);
%Normalize it and place the output in Output_Array

audiowrite('output.flac',Output_Array, Fs);
%Write Output_Array to a file


function arr_Output = thatwayRMSConverter(arr_A,K)

squaredData = arr_A .* arr_A;
% Square the input data
SummedData = smoothdata(squaredData, '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 and will result in the squaredRMS
% array.

% 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).

rootSummedData = arrayfun(@RootTaker,SummedData);
% This takes the root of the squared RMS array.

global MaxValRMS;
MaxValRMS = max(rootSummedData);

% 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,rootSummedData);
% 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_A);
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.

%We return the arr_output which will then be written to a file.

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

figure(1);
plot(arr_A);
% "Representation of Input Data"
figure(2);
plot(squaredData);
% "Squared input data"
hold on
plot(SummedData);
% This is the RMS-Squared
hold on
plot(rootSummedData);
% This is the root of the RMS-squared. This is the RMS.

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"


end


%% Functions

function output = RootTaker(datapoint)
    output = sqrt(datapoint);
end
% This returns the root of the data point.

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 quiet parts. By dividing this by low (quiet)
% 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.