function [frameFout, TNScoeffs] = TNS(frameFin, frameType)
%Implementation of the TNS step
%   Usage [frameFout, TNScoeffs] = TNS(frameFin, frameType), where:
%       Inputs
%       - frameFin is the frame in the frequency domain, in MDCT coefficients
%       representation containing both channels of the audio stored in an
%       array of dimensions 1024X2
%       - 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)
%
%       Output
%       - frameFout is the frame in the frequency domain after Temporal Noise
%       Shaping, in MDCT coefficients representation containing both channels
%       of the audio stored in an array of dimensions 1024X2
%       - TNScoeffs is the quantized TNS coefficients array of dimensions
%       4X8 for EIGHT_SHORT_SEQUENCE frames and 4X1 otherwise

    % Declares constant numbers of bands for long and short windows
    LONG_WINDOW_NUMBER_OF_BANDS = 69;
    SHORT_WINDOW_NUMBER_OF_BANDS = 42;
    % Declares constant order of the linear prediction filter
    LPF_ORDER = 4;
    % Declares constant coefficients' resolution
    COEF_RES = 4;

    % Declares persistent variable holding the TNS tables and initializes if empty
    persistent TNSTables;
    if isempty(TNSTables)
        TNSTables = load('TableB219.mat');
    end

    if ~strcmp(frameType, 'ESH')
        % Calculates the energy per band
        bandEnergy(LONG_WINDOW_NUMBER_OF_BANDS) = 0;
        for band = 1:LONG_WINDOW_NUMBER_OF_BANDS - 1
            bandEnergy(band) = sumsqr(frameFin(TNSTables.B219a(band, 2) + 1:TNSTables.B219a(band + 1, 2)));
        end
        bandEnergy(LONG_WINDOW_NUMBER_OF_BANDS) = sumsqr( ...
            frameFin(TNSTables.B219a(LONG_WINDOW_NUMBER_OF_BANDS, 2) + 1:...
            TNSTables.B219a(LONG_WINDOW_NUMBER_OF_BANDS, 3) + 1));

        % Calculates the normalization factors
        
        % PROBABLY WRONG, THE ONE BELLOW SHOULD BE RIGHT AND BETTER
%         normalizationFactor(1024) = 0;
%         for normIndex = 1:1024
%             bandIndices = find(TNSTables.B219a(:, 2) >= (normIndex - 1), 1);
%             if ~isempty(bandIndices)
%                 normalizationFactor(normIndex) = sqrt(bandEnergy(bandIndices));
%             else
%                 normalizationFactor(normIndex) = sqrt(bandEnergy(find(TNSTables.B219a(:, 3) >= normIndex - 1, 1)));
%             end
%         end

        bandIndices = quantiz(0:1023, TNSTables.B219a(:, 2) - 1);
        normalizationFactor = sqrt(bandEnergy(bandIndices));
        
        % Smooths normalization factors
        for normIndex = 1023:-1:1
            normalizationFactor(normIndex) = (normalizationFactor(normIndex) + normalizationFactor(normIndex + 1)) / 2;
        end
        for normIndex = 2:1024
            normalizationFactor(normIndex) = (normalizationFactor(normIndex) + normalizationFactor(normIndex - 1)) / 2;
        end

        % Normalizes MDCT coefficients according to each band energy
        normalizedFrameFin = frameFin ./ normalizationFactor';

        % Calculates the linear prediction coefficients
        [linPredCoeff, ~] = lpc(normalizedFrameFin, LPF_ORDER);

        % Quantizes these coefficients
        quantizedLinPredCoeff(LPF_ORDER) = 0;
        for coeffIndex = 2:length(linPredCoeff)
            quantizedLinPredCoeff(coeffIndex - 1) = asin(linPredCoeff(coeffIndex));
            if linPredCoeff(coeffIndex) >= 0
                quantizedLinPredCoeff(coeffIndex - 1) = round(quantizedLinPredCoeff(coeffIndex - 1) * ...
                    (bitshift(1, COEF_RES - 1) - 0.5) / (pi / 2));
                quantizedLinPredCoeff(coeffIndex - 1) = quantizedLinPredCoeff(coeffIndex - 1) / ...
                    (bitshift(1, COEF_RES - 1) - 0.5) / (pi / 2);
                quantizedLinPredCoeff(coeffIndex - 1) = sin (quantizedLinPredCoeff(coeffIndex - 1));
            else
                quantizedLinPredCoeff(coeffIndex - 1) = round(quantizedLinPredCoeff(coeffIndex - 1) * ...
                    (bitshift(1, COEF_RES - 1) + 0.5) / (pi / 2));
                quantizedLinPredCoeff(coeffIndex - 1) = quantizedLinPredCoeff(coeffIndex - 1) / ...
                    (bitshift(1, COEF_RES - 1) + 0.5) / (pi / 2);
                quantizedLinPredCoeff(coeffIndex - 1) = sin (quantizedLinPredCoeff(coeffIndex - 1));
            end
        end
        
        % Filters MDCT coefficients
        if ~isstable(1, [1 (quantizedLinPredCoeff * (-1))])
            error('TNS, l[79]: Inverse filter not stable!');
        else
            TNScoeffs = quantizedLinPredCoeff;
            frameFout = filter([1 (quantizedLinPredCoeff * (-1))], 1, frameFin);
        end
    else
        % Initializes output vectors
        TNScoeffs(LPF_ORDER * 8) = 0;
        frameFout(1024) = 0;
        
        bandEnergy(SHORT_WINDOW_NUMBER_OF_BANDS) = 0;
        for subFrameIndex = 1:8
            subFrame = frameFin((subFrameIndex - 1) * 128 + 1:subFrameIndex * 128);
            
            % Calculates the energy per band
            for band = 1:SHORT_WINDOW_NUMBER_OF_BANDS - 1
                bandEnergy(band) = sumsqr(subFrame(TNSTables.B219b(band, 2) + 1:TNSTables.B219b(band + 1, 2)));
            end
            bandEnergy(SHORT_WINDOW_NUMBER_OF_BANDS) = sumsqr( ...
                subFrame(TNSTables.B219b(SHORT_WINDOW_NUMBER_OF_BANDS, 2) + 1:...
                TNSTables.B219b(SHORT_WINDOW_NUMBER_OF_BANDS, 3) + 1));

            % Calculates the normalization factors
            normalizationFactor(128) = 0;
            for normIndex = 1:128
                bandIndices = find(TNSTables.B219b(:, 2) >= normIndex - 1, 1);
                if ~isempty(bandIndices)
                    normalizationFactor(normIndex) = sqrt(bandEnergy(bandIndices));
                else
                    normalizationFactor(normIndex) = sqrt(bandEnergy(find(TNSTables.B219b(:, 3) <= normIndex - 1, 1)));
                end
            end

            % Smooths normalization factors
            for normIndex = 127:-1:1
                normalizationFactor(normIndex) = (normalizationFactor(normIndex) + normalizationFactor(normIndex + 1)) / 2;
            end
            for normIndex = 2:128
                normalizationFactor(normIndex) = (normalizationFactor(normIndex) + normalizationFactor(normIndex - 1)) / 2;
            end

            % Normalizes MDCT coefficients according to each band energy
            normalizedFrameFin = subFrame ./ normalizationFactor';

            % Calculates the linear prediction coefficients
            [linPredCoeff, ~] = lpc(normalizedFrameFin, LPF_ORDER);

            % Quantizes these coefficients
            quantizedLinPredCoeff(LPF_ORDER) = 0;
            for coeffIndex = 1:length(linPredCoeff)
                quantizedLinPredCoeff(coeffIndex - 1) = asin(linPredCoeff(coeffIndex));
                if linPredCoeff(coeffIndex) >= 0
                    quantizedLinPredCoeff(coeffIndex - 1) = round(quantizedLinPredCoeff(coeffIndex - 1) * ...
                        (bitshift(1, COEF_RES - 1) - 0.5) / (pi / 2));
                    quantizedLinPredCoeff(coeffIndex - 1) = quantizedLinPredCoeff(coeffIndex - 1) / ...
                    (bitshift(1, COEF_RES - 1) - 0.5) / (pi / 2);
                    quantizedLinPredCoeff(coeffIndex - 1) = sin (quantizedLinPredCoeff(coeffIndex - 1));
                else
                    quantizedLinPredCoeff(coeffIndex - 1) = round(quantizedLinPredCoeff(coeffIndex - 1) * ...
                        (bitshift(1, COEF_RES - 1) + 0.5) / (pi / 2));
                    quantizedLinPredCoeff(coeffIndex - 1) = quantizedLinPredCoeff(coeffIndex - 1) / ...
                    (bitshift(1, COEF_RES - 1) + 0.5) / (pi / 2);
                    quantizedLinPredCoeff(coeffIndex - 1) = sin (quantizedLinPredCoeff(coeffIndex - 1));
                end
            end

            % Filters MDCT coefficients
            if ~isstable(1, [1 (quantizedLinPredCoeff * (-1))])
                error('TNS, l[79]: Inverse filter not stable!');
            else
                TNScoeffs((subFrameIndex - 1) * 4 + 1:subFrameIndex * 4) = quantizedLinPredCoeff;
                frameFout((subFrameIndex - 1) * 128 + 1:subFrameIndex * 128) = ...
                    filter([1 (quantizedLinPredCoeff * (-1))], 1, subFrame);
            end
        end
    end
end