function SMR = psycho(frameT, frameType, frameTprev1, frameTprev2)
%Implementation of Psychoacoustic Model
%   Usage SMR = psycho(frameT, frameType, frameTprev1, frameTprev2), where:
%       Inputs
%       - frameT is a frame in the time domain, containing only one of the
%       audio channels stored in a vector of length 2048
%       - frameType is the type of the current frame in string
%       representation, can be one of "OLS" (ONLY_LONG_SEQUENCE), "LSS"
%       (LONG_START_SEQUENCE), "ESH" (EIGHT_SHORT_SEQUENCE), "LPS"
%       (LONG_STOP_SEQUENCE)
%       - frameTprev1 is the previous frame in the time domain, containing
%       only one of the audio channels stored in a vector of length 2048
%       - frameTprev2 is the frame before frameTprev1 in the time domain,
%       containing only one of the audio channels stored in a vector of
%       length 2048
%
%       Output
%        - SMR is the signal to mask ratio array of dimensions 42X8 for
%       EIGHT_SHORT_SEQUENCE frames and 69X1 otherwise

    % Declares constant numbers of bands for long and short windows
    LONG_WINDOW_NUMBER_OF_BANDS = 69;
    SHORT_WINDOW_NUMBER_OF_BANDS = 42;

    % Declares persistent variable holding the TNS tables and initializes if empty
    persistent TNSTables spreadingLong spreadingShort hannLong hannShort;
    if isempty(TNSTables) || isempty(spreadingLong) || ...
            isempty(spreadingShort) || isempty(hannLong) || isempty(hannShort)
        TNSTables = load('TableB219.mat');

        spreadingLong = zeros(LONG_WINDOW_NUMBER_OF_BANDS);
        spreadingShort = zeros(SHORT_WINDOW_NUMBER_OF_BANDS);

        % Calculates all possible spreading function results for long
        % windows
        tmpx = repmat(TNSTables.B219a(:, 5)', LONG_WINDOW_NUMBER_OF_BANDS, 1);
        tmpx = tmpx - repmat(TNSTables.B219a(:, 5), 1, LONG_WINDOW_NUMBER_OF_BANDS);
        indeces = logical(tril(ones(LONG_WINDOW_NUMBER_OF_BANDS)));
        tmpx(indeces) = tmpx(indeces) .* 3;
        indeces = ~indeces;
        tmpx(indeces) = tmpx(indeces) .* 1.5;

        tmpz = 8 * min((tmpx - 0.5) .^ 2 - 2 * (tmpx - 0.5), 0);
        tmpy = 15.811389 + 7.5 * (tmpx + 0.474) - 17.5 * sqrt(1 + (tmpx + 0.474) .^ 2);

        tmpSum = tmpz + tmpy;
        spreadingLong(tmpy >= -100) = 10 .^ (tmpSum(tmpy >= -100) ./ 10);

        % Calculates all possible spreading function results for short
        % windows
        tmpx = repmat(TNSTables.B219b(:, 5)', SHORT_WINDOW_NUMBER_OF_BANDS, 1);
        tmpx = tmpx - repmat(TNSTables.B219b(:, 5), 1, SHORT_WINDOW_NUMBER_OF_BANDS);
        indeces = logical(tril(ones(SHORT_WINDOW_NUMBER_OF_BANDS)));
        tmpx(indeces) = tmpx(indeces) .* 3;
        indeces = ~indeces;
        tmpx(indeces) = tmpx(indeces) .* 1.5;

        tmpz = 8 * min((tmpx - 0.5) .^ 2 - 2 * (tmpx - 0.5), 0);
        tmpy = 15.811389 + 7.5 * (tmpx + 0.474) - 17.5 * sqrt(1 + (tmpx + 0.474) .^ 2);

        tmpSum = tmpz + tmpy;
        spreadingShort(tmpy >= -100) = 10 .^ (tmpSum(tmpy >= -100) ./ 10);

        hannLong = 0.5 - 0.5 * cos((pi * (0:2047 + 0.5)) / 1024);
        hannShort = 0.5 - 0.5 * cos((pi * (0:255 + 0.5)) / 128);
    end

    if ~strcmp(frameType, 'ESH')
        % Applies window to the frames
        windowedFrameT = frameT .* hannLong';
        windowedFrameTprev1 = frameTprev1 .* hannLong';
        windowedFrameTprev2 = frameTprev2 .* hannLong';

        % Calculates the FFT of each frame
        frameF = fft(windowedFrameT);
        frameFMag = 2 * abs(frameF(1:1024));
        frameFPhase = angle(frameF(1:1024));

        frameFrameFprev1 = fft(windowedFrameTprev1);
        frameFrameFprev1Mag = 2 * abs(frameFrameFprev1(1:1024));
        frameFrameFprev1Phase = angle(frameFrameFprev1(1:1024));

        frameFrameFprev2 = fft(windowedFrameTprev2);
        frameFrameFprev2Mag = 2 * abs(frameFrameFprev2(1:1024));
        frameFrameFprev2Phase = angle(frameFrameFprev2(1:1024));

        % Calculates the predicted magnitude and phase compontents of each
        % frequency
        magPred = 2 .* frameFrameFprev1Mag - frameFrameFprev2Mag;
        phasePred = 2 .* frameFrameFprev1Phase - frameFrameFprev2Phase;

        % Calculates this frame's predictability
        framePredictability = sqrt((frameFMag .* cos(frameFPhase) - ...
            magPred .* cos(phasePred)) .^ 2 + ...
            (frameFMag .* sin(frameFPhase) - ...
            magPred .* sin(phasePred)) .^ 2) ./ ...
            (frameFMag + abs(magPred));

        % Calculates the energy and weighted predictability in the
        % threshold calculation partitions
        bandEnergy(LONG_WINDOW_NUMBER_OF_BANDS) = 0;
        bandPredictability(LONG_WINDOW_NUMBER_OF_BANDS) = 0;
        for band = 1:LONG_WINDOW_NUMBER_OF_BANDS
            bandEnergy(band) = sumsqr(frameFMag(TNSTables.B219a(band, 2) + 1: ...
                TNSTables.B219a(band, 3) + 1));
            bandPredictability(band) = sumsqr(frameFMag( ...
                TNSTables.B219a(band, 2) + 1:TNSTables.B219a(band, 3) + 1) .* ...
                framePredictability(TNSTables.B219a(band, 2) + 1: ...
                TNSTables.B219a(band, 3) + 1));
        end

        % Convolves the partitioned energy and predictability with the
        % spreading function
        ecb = sum(bandEnergy .* spreadingLong', 2);
        ct = sum(bandPredictability .* spreadingLong', 2);

        % Renormalizes values
        cb = ct ./ ecb;
        en = ecb ./ sum(spreadingLong, 1);
        
        % Calculates the tonality index
        tb = -0.299 - 0.43 .* log(cb);
    end
end