function [ T , P ] = psycho( frameT , frameType )

%PSYCHO Calculates the psychoacoustic thresholds and the energy spectrum

%   Inputs: frameT  : the frame array with the samples in time

%                   : the type of the frame

%   Outputs: T: a <numberofbands X numberofchannels> matrix with the

%               threshold for every band.

%            P: a <numberofbands X numberofchannels> with the energy

%               spectrum of every band





%We check what type of frame this is

if strcmp(frameType,'ESH')==1

    %defining the bjs, keep in mind that I added another value in the end

    %of the array b(43)=127

    b=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 19 21 23 25 27 29 31 34 37 40 43 46 50 54 58 63 68 74 80 87 95 104 114 126 128];

    %defining the number of bands

    numberofbands=42;

    %defining B, keep in mind that I added another value in the end

    %of the array B(43)=24

    B=[0 1.88 3.70 5.39 6.93 8.29 9.49 10.53 11.45 12.26 12.96 13.59 14.15 14.65 15.11 15.52 15.90 16.56 17.15 17.66 18.13 18.54 18.93 19.28 19.69 20.14 20.54 20.92 21.27 21.64 22.03 22.39 22.76 23.13 23.49 23.5 24.00 24.00 24.00 24.00 24.00 24.00 24.00];

    %defining Tquiet

    Tquiet=[27.28 22.28 14.28 12.28 12.28 12.28 12.28 12.28 12.28 12.28 12.28 12.28 12.28 12.28 12.28 12.28 12.28 15.29 15.29 15.29 15.29 15.29 15.29 15.29 17.05 20.05 20.05 20.05 20.05 23.30 28.30 28.30 29.27 39.27 40.06 40.06 50.73 51.31 51.82 52.28 53.07 53.07];

else

    %defining the bjs, keep in mind that I added another value in the end

    %of the array b(70)=1023

    b=[0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 41 44 47 50 53 56 59 62 66 70 74 78 82 87 92 97 103 109 116 123 131 139 148 158 168 179 191 204 218 233 249 266 284 304 325 348 372 398 426 457 491 528 568 613 663 719 782 854 938 1024];

    %defining the number of bands

    numberofbands=69;

    %defining B, keep in mind that I added another value in the end

    %of the array B(43)=24

    B=[0.24 0.71 1.18 1.65 2.12 2.58 3.03 3.48 3.92 4.35 4.77 5.19 5.59 5.99 6.37 6.74 7.10 7.45 7.80 8.20 8.68 9.13 9.55 9.96 10.35 10.71 11.06 11.45 11.86 12.25 12.62 12.96 13.32 13.70 14.05 14.41 14.77 15.13 15.49 15.85 16.20 16.55 16.91 17.25 17.59 17.93 18.28 18.62 18.96 19.30 19.64 19.97 20.31 20.95 20.99 21.33 21.66 21.99 22.32 22.66 23.00 23.33 23.67 24.00 24.00 24.00 24.00 24.00 24.00 24.00];

    %defining Tquiet

    Tquiet=[40.29 40.29 35.29 35.29 32.29 32.29 27.29 27.29 27.29 25.29 25.29 25.29 25.29 25.29 25.29 25.29 25.29 25.29 25.29 27.05 27.05 27.05 27.05 27.05 27.05 27.05 27.05 28.30 28.30 28.30 28.30 28.30 29.27 29.27 29.27 30.06 30.06 30.73 30.73 31.31 31.31 31.82 32.28 32.28 32.69 33.07 33.42 33.74 34.04 34.32 34.58 34.83 38.29 38.50 38.89 41.08 41.43 41.75 47.19 47.59 47.96 58.30 28.81 69.27 69.76 70.27 70.85 71.52 70.20];

end





%Calling filterbank

%define the wintype

winType='KAIZER';

frameF=filterbank(frameT,frameType,winType);





%repeat for each channel

channels=2;

for channel=1:1:channels

    



    %Calculating the energy spectrum Pi(j)

    if strcmp(frameType,'ESH')==1

        for i=0:1:7

            for j=1:1:numberofbands

                P(j+i*42,channel)=0;

                for k=(b(j)+1):1:b(j+1)

                    P(j+i*42,channel)=P(j+i*42,channel)+ frameF(channel,k+i*128)*frameF(channel,k+i*128); 

                end

            end

        end

    else

        for j=1:1:numberofbands

        P(j,channel)=0;

            for k=(b(j)+1):1:b(j+1)

                P(j,channel)=P(j,channel)+ frameF(channel,k)*frameF(channel,k);

            end

        end

    end

end



%Call the TNS level

[frameF,TNScoeffs]=TNS(frameF,frameType,P);





%repeat for each channel

for channel=1:1:channels

    



    %ReCalculating the energy spectrum Pi(j)

    if strcmp(frameType,'ESH')==1

        for i=0:1:7

            for j=1:1:numberofbands

                P(j+i*42,channel)=0;

                for k=(b(j)+1):1:b(j+1)

                    P(j+i*42,channel)=P(j+i*42,channel)+ frameF(channel,k+i*128)*frameF(channel,k+i*128); 

                end

            end

        end

    else

        for j=1:1:numberofbands

        P(j,channel)=0;

            for k=(b(j)+1):1:b(j+1)

                P(j,channel)=P(j,channel)+ frameF(channel,k)*frameF(channel,k);

            end

        end

    end

end



%repeat for each channel

for channel=1:1:channels

    %Calculate the acoustic threshold Ti(j)

    if strcmp(frameType,'ESH')==1

        for i=0:1:7

            for j=2:1:numberofbands

                DB=(B(j+1)+B(j-1))/2;

                sl(j)=20*DB;

                sh(j)=15*DB;

                Tsc(j+i*42)=10*log10(P(j+i*42,channel))-29;

                T(j+i*42,channel)=max(Tsc(j+i*42),Tsc(j+i*42-1)-sh(j));

            end

        end

    else

        for j=2:1:numberofbands

            DB=(B(j+1)+B(j-1))/2;

            sl(j)=30*DB;

            sh(j)=20*DB;

            Tsc(j)=10*log10(P(j,channel))-29;

            T(j,channel)=max(Tsc(j),Tsc(j-1)-sh(j));

        end

    end

    

    %Calculate the Ti,spread

    if strcmp(frameType,'ESH')==1

        for i=0:1:7

            for j=1:1:(numberofbands-1)

                Tspread(j+i*42)=max(T(j+i*42,channel),T(j+i*42+1,channel)-sl(j));

            end

            Tspread(numberofbands+i*42)=T(numberofbands+i*42,channel);

        end

    else

        for j=1:1:(numberofbands-1)

            Tspread(j)=max(T(j,channel),T(j+1,channel)-sl(j));

        end

        Tspread(numberofbands)=T(numberofbands,channel);

    end

    

    

    if strcmp(frameType,'ESH')==1

        for i=0:1:7

            for j=1:1:numberofbands

                T(j+i*42,channel)=max(Tspread(j+i*42),Tquiet(j));

            end

        end

    else

        for j=1:1:numberofbands

            T(j,channel)=max(Tspread(j),Tquiet(j));

        end

    end



end



end