%DEMODULATESOUND  Takes raw recorded audio and attempts to recreate a message created using the modulateData function
%
% PARAMETERS:
%  sampleRate: 
%       The sample rate of the recorded sound. ALl the symbols need to match this sample rate
%  audio:
%     The raw audio data to demodulate
%  dataSymbols:
%     Symbols used to represent different data values within the message. If
%     a data symbol with the index 1 is found, a 1 is added to the decoded
%     message, a symbol with index 2 adds a 2 and so on.
%  pulseDistance:
%     Amount of samples between two data symbols.
%  dataLength:
%     The amount of symbols to look for in the message
%  startSymbol (Optional):
%     The symbol used to indicate the start of the message
%  startPadding (Optional):
%     The amount of samples between the start symbol and the first data symbol
%   
% RETURNS:
%  The final decoded message
%   
%  See also: MODULATEDATA
function decodedMsg = demodulateSound(sampleRate, audio, dataSymbols, pulseDistance, dataLength, startSymbol, startPadding)
    DEFAULT_PADDING = 4000;


    %Default parameters
    switch nargin
        case 5 %No start symbol
            startSymbol = zeros(1,1);
            startPadding = 0;
        case 7
        otherwise
            error('Failed to demodulate sound, wrong amount of parameters');
    end

    %Symbols for displaying points
    locationSymbol = 'gx'; %Change the g if you want other colors
    startLocationSymbol = 'dk'; %change d to set color

    %Analyzing the start symbol
    startFrequency = getFrequencyOfSample(sampleRate, startSymbol);
    startConv = conv(audio, fliplr(startSymbol), 'same');
    startAmplitude = getAmplitudeOverTime(sampleRate, startConv, 100);

    %Doing the same for all the data symbols
    convData = {};
    symbolAmplitudes = {};
    for n = 1:length(dataSymbols);
        %Calculating the frequency of a sample
        symbolFreq = getFrequencyOfSample(sampleRate, dataSymbols{n});

        %Calculating the correlation between the different symbols and the
        %recorded audio
        convData(n) = {conv(audio, dataSymbols{n}, 'same')};
        
        %Finding the 'amplitude' of each symbol over time
        symbolAmplitudes(n) = {getAmplitudeOverTime(sampleRate, convData{n}, symbolFreq)};
    end

    %Figguring out where the signal starts
    [maxStartAmp, signalStart] = max(startAmplitude);
    
    %The starting position of the actual data
    symbolStart = signalStart + startPadding;

    %Expected locations of the symbols
    symbolLocations = symbolStart:pulseDistance:symbolStart + pulseDistance * dataLength - 1;

    %Decoding the message
    %This is done by checking which symbol has the highest amplitude in each
    %locaction where a bit should be found
    decodedMsg = [];
    for n = symbolLocations
        %Making sure all the bit locations were recorded
        if(length(symbolAmplitudes{1}) < n)
            fprintf('Some expected bit locations are outside recorded audio. This can be caused by recording time being to short or the start symbol not being found propperly\n');
            fprintf('Look at the raw recorded audio figure to see what caused the issue.\n');
            fprintf('Setting remaining symbols to 1\n');
            break;
        end

        highestAmp = getHighestAmplitudeInPoint(symbolAmplitudes, round(n));

        decodedMsg(end + 1) = highestAmp;
    end

    %If there are missing bits
    remainingBits = dataLength - length(decodedMsg);
    if(remainingBits > 0)
        appended = ones(1, remainingBits);
        decodedMsg = [decodedMsg, appended];
    end

    %Subtract one from the result to make it start at 0
    decodedMsg = decodedMsg - 1;
    
    %Plot the audio data along with the expected location of each symbol
    figure(10);
    plot(audio);
    hold on;
        plot(signalStart, 0, startLocationSymbol);
        plot(symbolLocations, 0, locationSymbol);
    hold off;
    title('Raw recorded audio');
    

    %Plot the data after convolution
    plotLegends = {};

    figure(11);
    plot(signalStart, 0, startLocationSymbol);
    plotLegends(end + 1) = {'Start location'};
    hold on;
        for n = 1:length(convData)
            plot(convData{n});
            plotLegends(end + 1) = {['Symbol ', num2str(n - 1)]};
        end
        plot(startConv);
        plotLegends(end + 1) = {'Start symbol'};
        plot(symbolLocations, 0, locationSymbol);
        plotLegends(end + 1) = {'Bit location'};
    hold off;
    title('Symbol correlation');
    legend(plotLegends);
    
    plotLegends = {};
    %Plotting the amplitude of the convoluted data
    figure(12);
    plot(signalStart, 0, startLocationSymbol);
    plotLegends(end + 1) = {'Start location'};
    hold on;
        for n = 1:length(symbolAmplitudes)
            plot(symbolAmplitudes{n});
            plotLegends(end + 1) = {['Symbol ', num2str(n - 1)]};
        end
        plot(startAmplitude);
        plotLegends(end + 1) = {'Start symbol'};
        plot(symbolLocations, 0, locationSymbol);
        plotLegends(end + 1) = {'Bit location'};
    hold off;
    title('Amplitude of symbol correlation');
    legend(plotLegends);
end

%Returns the main frequency of a sample
%sampleRate: the sample rate of the sample
%sample: the sample to check
function frequency = getFrequencyOfSample(sampleRate, sample)
    [maxValue, maxIndex] = max(abs(fft(sample)));

    frequency = maxIndex * sampleRate / length(sample);
end

%Returns the amplitude of a signal over time. Normalized
%sampleRate: The sample rate of the data
%data: The signal to calculate the amplitude of
%frequency: The main frequency of the signal, used to calculate the cutoff
%       frequency for the filter
function result = getAmplitudeOverTime(sampleRate, data, frequency)
    flatData = abs(data);
    
    Fnorm = frequency / 50 / (sampleRate / 2);
    lpFilt = designfilt('lowpassfir', 'FilterOrder', 40, 'CutoffFrequency', Fnorm);
    ampFiltered  = filtfilt(lpFilt, flatData);
    
    result = ampFiltered;
end

function result = getRawAmplitudeOverTime(sampleRate, data, frequency)
    flatData = abs(data);

    Fnorm = frequency / 25 / (sampleRate / 2);
    lpFilt = designfilt('lowpassfir', 'FilterOrder', 30, 'CutoffFrequency', Fnorm);
    ampFiltered  = filtfilt(lpFilt, flatData);
    
    result = ampFiltered;
end

%Returns the index of the curve with the highest amplitude in a specific ponit.
%curves: Array of curves to compare
%point: X coordinate to check
function [highestIndex, highestValue] = getHighestAmplitudeInPoint(curves, point)
    highestIndex = 1;
    highestValue = -1;

    for n = 1:length(curves)

        if curves{n}(point) > highestValue 
            highestIndex = n;
            highestValue = curves{n}(point);
        end
    end
end

