% funkce pro zjištění průběhu základního kmitočtu analyzované nahrávky
%
% význam vstupních proměnných:
% ============================
% y vstupní signál analyzované nahrávky
% Fvz je vzorkovací kmitčet nahrávky
% reference určuje kmitočet referenční noty a1, které odpovídá index n = 0 a od něj se odvíjí ostatní noty
%
% význam výstupních proměnných:
% =============================
% kmitocet je průběh detekovaného kmitočtu
% noty je odpovídající průběh detekovaných not
% t je odpovídající časová osa
% pulsy je vektor pořadí vzorků hlasivkových pulsů v analyzovaném signálu

function [kmitocet, noty, t, pulsy] = Analyza(y, Fvz, reference)

    bar = waitbar(0, 'Probíhá analýza ...', 'WindowStyle', 'modal'); % ukazatel fáze dokončení analýzy

    %% Dolní propust
    fm = 1100;  % mezní kmitočet
    rad = 6;    % řád filtru

    [b,a] = butter(rad, fm/(Fvz/2));

    yfilt = filter(b,a,y);
    
    %% Detekce částí signálu se zpěvem a bez
    N = length(yfilt); % celkový počet vzorků v nahrávce
    ts = 0.001; % [s] délka segmentu, na kterém se bude určovat energie signálu
    
    Nseg = round(Fvz * ts); % počet vzorků v segmentu
    seg = floor(N/Nseg);  % počet celých segmentů, který lze získat
    
    for i = 1:seg
        energie = sum(yfilt((i-1)*Nseg+1:(i-1)*Nseg+Nseg).^2); % výpočet krátkodobé energie signálu
        if energie < 0.02
            yfilt((i-1)*Nseg+1:(i-1)*Nseg+Nseg) = 0; % vynulování signálu s energií pod definovanou úrovní
        end
    end
   
    %% Segmentace
                
    N = length(y); % celkový počet vzorků v nahrávce 
    ts = 0.025; % [s] délka segmentu
    krok = ts/2; % [s] krok posuvu okénka
    
    Nseg = round(Fvz * ts); % počet vzorků v segmentu
    Nkrok = round(Fvz * krok);  %počet vzorků na jeden krok mezi segmenty
    seg = floor((N-Nseg)/Nkrok);  %počet segmentů s překryvem, který lze získat
    
    y = y(1:(seg*Nkrok+Nseg)); %ořezání signálu, aby to vyšlo na celé segmenty
    yfilt = yfilt(1:(seg*Nkrok+Nseg)); %ořezání signálu, aby to vyšlo na celé segmenty

    segmenty = zeros(Nseg,floor((length(y)-Nseg)/Nkrok)); % vytvoření prázdné matice pro naplnění nasegmentovaným signálem 
    segmentyfilt = zeros(Nseg,floor((length(yfilt)-Nseg)/Nkrok)); % vytvoření prázdné matice pro naplnění nasegmentovaným signálem dolní propustí
    kmitocet = zeros(1, seg); % vytvoření prázdného vektoru pro naplnění výsledným průběhem kmitočtu
    noty = zeros(1, seg); % vytvoření prázdného vektoru pro naplnění výsledným průběhem not

    % Segmentace s pravúhlým okénkem
    for i = 1:seg
        segmenty(:,i) = y((i-1)*Nkrok+1:(i-1)*Nkrok+Nseg); % signál neproběhlý dolní propustí
    end                                                    % pro hledání polohy pulsů, protože
                                                           % poloha pulsů se musí hledat v nijak nepoškozeném signálu
    
    for i = 1:seg
        segmentyfilt(:,i) = yfilt((i-1)*Nkrok+1:(i-1)*Nkrok+Nseg); % signál filtorvaný dolní propustí
    end                                                            % pro analýzu
    
    %% Váhovaná autokorelační funkce
    
    waitbarkrok = 0.8/seg; % obsluha waitbaru
    aktualnikrok = 0;

    VACFseg = zeros(Nseg, seg);

    for i = 1:seg
        VACFseg(:,i) = ACF(segmentyfilt(:,i))./(AMDF(segmentyfilt(:,i))+1);
        
        aktualnikrok = aktualnikrok + waitbarkrok; % obsluha waitbaru
        waitbar(aktualnikrok);
    end
    

    %% Vyčtení frekvence z průběhu váhované autokorelační funkce každého segmentu
    waitbarkrok = 0.05/seg; % obsluha waitbaru    
    
    for i = 1:seg  
        [pks,locs] = findpeaks(VACFseg(:,i)); % najde v průběhu autokorelační funkce každého segmentu maxima
        minpeakweight = 0.3 * max(pks);
        [pks,locs] = findpeaks(VACFseg(:,i), 'minpeakheight', minpeakweight);
        
        if length(locs) < 1         % pokud je počet maxim v segmentu menší než 2, nelze vypočítat kmitočet (pozice prvníhoo maxima není v locs obsaženo)
            f0 = 0;
        else                        % výpočet časové úseku mezi pozicí prvního maxima a prvním vzorkem
            Tpeak = (locs(1)-1)/Fvz;

            f0 = 1/Tpeak;       %výsledná frekvence základního tonu daného segmentu
        end
        
        kmitocet(1,i) = f0; % naplnění hodnotou základního kmitočtu pro každý segment
        noty(1,i) = log2((f0/reference)^12); % výpočet indexu noty odpovídajícímu kmitočtu
    
        if f0 == 0 % pokud nebylo možné určit kmitočet, f0 = 0 a hodnota hledaného kmitočtu bude mimo graf
            kmitocet(1,i) = -50;
            noty(1,i) = -50;
        end 
        
        aktualnikrok = aktualnikrok + waitbarkrok; % obsluha waitbaru
        waitbar(aktualnikrok);
    end
    
    %% Filtrace mediánovým filtrem
    kmitocet = medfilt1(kmitocet,5); % pro vyfiltrování prudkých změn detekovaného kmitočtu
    noty = medfilt1(noty,5);         % zajišťující větší přesnost detekce

    %% Detekce hlasivkových pulsů v analyzovaném signálu
    waitbarkrok = 0.15/seg; % obsluha waitbaru
    
    pulsy = zeros(4 * floor(Nseg/(Fvz/max(kmitocet))) * seg, 1); % vytvoření prázdného vektoru pro konečné naplnění pozicemi jednotlivých pulsů
    
    konstanta = 0.8; % určuje šířku oblasti pro vyhledávání dalšího pulsu ( 0.8 => šířka oblasti 0.2 periody na každou stranu)
                                                                          % 0.5 => šířka oblasti 0.5 periody na každou stranu)    
    for i = 1:seg
        if kmitocet(i) > 0
            [pks,locs] = findpeaks(segmenty(:,i)); % nalezení všech maxim v segmentu
            
            if ~isempty(pks)
                pulsyseg = zeros((4 * floor(Nseg/(Fvz/kmitocet(i)))), 1);   % vytoření prázdného vektoru o velikosti 4*počet možných pulsů při nejvyšším kmitočtu
                                                                 % úmyslně větší, než je očekávaný počet možných pulsů
                
                % --- nalezení jednoho maxima ---
                index = locs(pks == max(pks)); % zjištění pozice maximálního pulsu v daném segmentu
                index = index(1); % omezení se pouze na jednu pozici, pokud v signálu budou dvě maxima se stejnou hodnotou
                pocatecniindex = index;
                
                pulsyseg(2*floor(Nseg/(Fvz/kmitocet(i)))) = pocatecniindex;   % vložení pozice prvního nalezeného pulsu do středu prázdného vektoru
                
                pocatecnipozice = 2*floor(Nseg/(Fvz/kmitocet(i))); % nastavení počáteční hodnotu pro indexování ve vektoru
                
                % --- označení všech pulsů nad počátečním pulsem ---
                j = 1;  % nastavení počáteční pozice indexování
                                                                      
                while locs(end) > index + (2 - konstanta)*(Fvz/kmitocet(i));     % dokud se šířka vybírací oblasti vleze do vektoru locs
                    spodnimez = index + konstanta*(Fvz/kmitocet(i)) - 1;         % nastavení mezí vybírací oblasti
                    hornimez = index + (2 - konstanta)*(Fvz/kmitocet(i)) - 1;

                    moznepozice = locs((locs >= spodnimez) & (locs <= hornimez)); % pozice možných vrcholů [poradi vzorku] v očekávatelném úseku signálu
                    ocekavanapozice = index + (Fvz/kmitocet(i)) - 1;      % ocekavaná pozice vzorku ( + 1 perioda za aktuální pozicí)
                    
                    if isempty(moznepozice) % pokud se v hledaném rozsahu nenachází puls
                        skutecnapozice = ocekavanapozice; % určí se očekávaná pozice
                    else
                        rozdil = moznepozice - ocekavanapozice;    % nalezení nejbližšího nalezeného pulsu k očekávanému
                        skutecnapozice = ocekavanapozice + rozdil(abs(rozdil) == min(abs(rozdil)));
                        skutecnapozice = skutecnapozice(1); % omezení se pouze na jednu pozici, když bude rozdíl na obě strany stejný
                    end

                    index = skutecnapozice; % zapsání nové aktuální pozice
    
                    pulsyseg(pocatecnipozice + j) = skutecnapozice; % zapsání nové pozice do výsledného vektoru
    
                    j = j + 1;  % posun indexu
                end
                
                if locs(end) >= index + (Fvz/kmitocet(i)); % označení posledního pulsu (pro poslední puls se musí zhora omezit vybírací meze)
                    spodnimez = index + konstanta*(Fvz/kmitocet(i)) - 1;     % nastavení mezí vybírací oblasti
                    hornimez = index + (Fvz/kmitocet(i));                    % horni mez není ovlivněna konstantou ( shora je vybírací oblast nulová)

                    moznepozice = locs((locs >= spodnimez) & (locs <= hornimez)); % pozice možných vrcholů [poradi vzorku] v očekávatelném úseku signálu
                    ocekavanapozice = index + (Fvz/kmitocet(i)) - 1;
                    
                    if isempty(moznepozice) % pokud se v hledaném rozsahu nenachází puls
                        skutecnapozice = ocekavanapozice; % určí se očekávaná pozice
                    else
                        rozdil = moznepozice - ocekavanapozice;    % nalezení nejbližšího nalezeného pulsu k očekávanému
                        skutecnapozice = ocekavanapozice + rozdil(abs(rozdil) == min(abs(rozdil)));
                        skutecnapozice = skutecnapozice(1); % omezení se pouze na jednu pozici, když bude rozdíl na obě strany stejný
                    end

                    pulsyseg(pocatecnipozice + j) = skutecnapozice;    % zapsání nové pozice do výsledného vektoru
                end

                % --- označení všech pulsů pod počátečním pulsem ---
                j = 1;  % nastavení počáteční pozice indexování
                index = pocatecniindex;

                while locs(1) < index - (2 - konstanta)*(Fvz/kmitocet(i));       % dokud se šířka vybírací oblasti vleze do vektoru locs
                    spodnimez = index - (2 - konstanta)*(Fvz/kmitocet(i));       % nastavení mezí vybírací oblasti
                    hornimez = index - konstanta*(Fvz/kmitocet(i));

                    moznepozice = locs((locs >= spodnimez) & (locs <= hornimez)); % pozice možných vrcholů [poradi vzorku] v očekávatelném úseku signálu
                    ocekavanapozice = index - (Fvz/kmitocet(i));

                    if isempty(moznepozice) % pokud se v hledaném rozsahu nenachází puls
                        skutecnapozice = ocekavanapozice; % určí se očekávaná pozice
                    else
                        rozdil = moznepozice - ocekavanapozice;    % nalezení nejbližšího nalezeného pulsu k očekávanému
                        skutecnapozice = ocekavanapozice + rozdil(abs(rozdil) == min(abs(rozdil)));
                        skutecnapozice = skutecnapozice(1); % omezení se pouze na jednu pozici, když bude rozdíl na obě strany stejný
                    end
    
                    index = skutecnapozice;     % zapsání nové aktuální pozice
    
                    pulsyseg(pocatecnipozice - j) = skutecnapozice;    % zapsání nové pozice do výsledného vektoru
    
                    j = j + 1;  % posun indexu
                end

                if locs(1) <= index - (Fvz/kmitocet(i)); % označení posledního pulsu (pro poslední puls se musí zdola omezit vybírací meze)
                    hornimez = index - konstanta*(Fvz/kmitocet(i));      % nastavení mezí vybírací oblasti
                    spodnimez = index - (Fvz/kmitocet(i));               % spodni mez není ovlivněna konstantou ( zdola je vybírací oblast nulová)

                    moznepozice = locs((locs >= spodnimez) & (locs <= hornimez)); % pozice možných vrcholů [poradi vzorku] v očekávatelném úseku signálu
                    ocekavanapozice = index - (Fvz/kmitocet(i));

                    if isempty(moznepozice) % pokud se v hledaném rozsahu nenachází puls
                        skutecnapozice = ocekavanapozice; % určí se očekávaná pozice
                    else
                        rozdil = moznepozice - ocekavanapozice;    % nalezení nejbližšího nalezeného pulsu k očekávanému
                        skutecnapozice = ocekavanapozice + rozdil(abs(rozdil) == min(abs(rozdil)));
                        skutecnapozice = skutecnapozice(1); % omezení se pouze na jednu pozici, když bude rozdíl na obě strany stejný
                    end

                    pulsyseg(pocatecnipozice - j) = skutecnapozice;    % zapsání nové pozice do výsledného vektoru
                end
                
                pulsy(length(pulsy(pulsy~=0))+1:length(pulsy(pulsy~=0))+length(pulsyseg(pulsyseg~=0))) = pulsyseg(pulsyseg~=0)+Nkrok*(i-1);
            end
        end
        
        aktualnikrok = aktualnikrok + waitbarkrok; % obsluha waitbaru
        waitbar(aktualnikrok);
    end
    
    pulsy = unique(pulsy(pulsy~=0)); % vymazání vícenásobných pulsů vzniklých při skládání
                                     % výsledků ze segmentů
      
    % korekce nepřesné detekce pulsů
    for i = 1:length(pulsy)-3
        if (pulsy(i) ~= 0) & (ceil(pulsy(i+1)/Nkrok) <= length(kmitocet))
            zmerenaperioda =  pulsy(i+1) - pulsy(i);
       
            if zmerenaperioda < 0.4*Fvz/(kmitocet(ceil(pulsy(i+1)/Nkrok))) % příliš blízké pulsy zprůměruje
                pulsy(i) = round((pulsy(i+1) + pulsy(i) )/2);
                pulsy(i+1) = 0;
            elseif (zmerenaperioda < 0.9*Fvz/(kmitocet(ceil(pulsy(i+1)/Nkrok)))) & ((pulsy(i+2)-pulsy(i)) < 1.1*Fvz/(kmitocet(ceil(pulsy(i+1)/Nkrok))))
                pulsy(i+1) = 0; % pokud je druhý pulz blízko a třetí lépe odpovídá detekované periodě
            end                 % druhý pulz se vymaže
        end
    end
    
    pulsy = pulsy(pulsy~=0); % odstraní vymazané pulzy
    
    %% Zpřesnění kmitočtu pomocí detekovaných pulzů - spíš zhoršuje vlivem ménší přesnosti určení polohy pulzů
%     kmitocet(:) = -50; % vynulování předchozí detekce
%     noty(:) = -50; % vynulování předchozí detekce
%     
%     for i = 1:length(pulsy)-1
%         Tpeak = (pulsy(i+1) - pulsy(i))/Fvz; % změřená perioda základního kmitočtu
%         f0 = 1/Tpeak;       % výsledná frekvence základního tonu daného segmentu
%         
%         if floor(pulsy(i)/Nkrok) >= 1
%             kmitocet(floor(pulsy(i)/Nkrok)) = f0; % naplnění hodnotou základního kmitočtu pro každý segment
%             noty(floor(pulsy(i)/Nkrok)) = log2((f0/reference)^12); % výpočet indexu noty odpovídajícímu kmitočtu
%         end
%     end
%     
%     %% Druhá filtrace mediánovým filtrem
%     kmitocet = medfilt1(kmitocet,5); % filtrace mediánovým filtrem
%     noty = medfilt1(noty,5); % filtrace mediánovým filtrem
    
    %% časová osa pro zobrazení v grafu
    t = 0:krok:krok*(seg-1);
    
    %% ukončení ukazatele fáze dokončení analýzy   
    waitbar(1);
    close(bar);

end