%83764
%Tomáš Heitel

function [a]=tonality(S,sAll1,cfs)

cfs2=cfs.^2;

%Prahování tvrdé
for j=1:S
    for i=1:sum(sAll1)
        if cfs2(i,j)>=0.001
            cfs2(i,j)=cfs2(i,j);
        else
            cfs2(i,j)=0.001;
        end
    end
end
%Výpočet PSD
for j=1:S
    for i=2:length(sAll1)
        while i <= length(sAll1)
        d=1:(i-1);
        Pyy(sum(sAll1(i-d))+1:sAll1(i)+sum(sAll1(i-d)),j) = (cfs2(sum(sAll1(i-d))+1:sAll1(i)+sum(sAll1(i-d)),j))/(max(abs(cfs2(sum(sAll1(i-d))+1:sAll1(i)+sum(sAll1(i-d)),j))));
        i=i+1;
        end
    end
end

%Výpočet spektrální plochosti sfm a tonálního parametru "a"
ug=zeros(length(sAll1),S);
ua=zeros(length(sAll1),S);
for j=1:S
    for i=2:length(sAll1)
        d=1:(i-1);
        ug(i,j)=prod(Pyy((sum(sAll1(i-d))+1):sAll1(i)+sum(sAll1(i-d)),j))^(1/sAll1(i));
        ua(i,j)=sum(Pyy((sum(sAll1(i-d))+1):sAll1(i)+sum(sAll1(i-d)),j))/sAll1(i);
        sfm(i,j)=ug(i,j)/ua(i,j);
        sfmdB(i,j)=10*log10(sfm(i,j));
        sfmm=-25;
            if sfmdB(i,j)/sfmm >= 1
                lambda(i,j) = 1;
            else
                lambda(i,j)=sfmdB(i,j)/sfmm;
            end
        a(i,j)=lambda(i,j)*(-0.275*(i-1)-15.025)+(1-lambda(i,j))*(-9); 
    end
end

a=a(2:26,:);