clear all;
close all;

%intro
fprintf('\n*******************************************************\n');
fprintf('*                                                     *\n');
fprintf('*      Algoritmus na odhalenie period funkcie         *\n');
fprintf('*                                                     *\n');
fprintf('*******************************************************\n');

% *** SETTNIGS *** 
cas=10*pi;        %kolko casu budeme vidiet funckiu

% *** zistime presnot a potrebny pocet vzoriek ***
pocet_vzoriek=100;
max_rozdiel=100;
rozkmit=1;
while (max_rozdiel>(0.0005*rozkmit))
    pocet_vzoriek=floor(pocet_vzoriek*1.05);
    vzorky=linspace(0,cas,pocet_vzoriek);
    % -- par skusobnych funkcii --%
    %f=abs(cos(vzorky));
    %f=sin(vzorky)+cos(4*vzorky);
    f=sin(2*vzorky)+cos(4*vzorky);
    %f=5*vzorky;
    max_rozdiel=0; %v podstate max rozdiel
    rozkmit=abs(max(f)-min(f)); %rodiel minimum a maximum
    for i=1:pocet_vzoriek-1,
        tmp=abs(f(i)-f(i+1));
        if (tmp>max_rozdiel)
               max_rozdiel=tmp;
        end;
    end;
end;

%nasledujucich 6 riadkov odkomentovat ak je potrebne zadavat funkciu rucne
% max_rozdiel=0;
% f=[ones(1,100),zeros(1,100),ones(1,100),zeros(1,100),ones(1,100)];
% % f=[1,2,3,4,5,6,1,2,3,4,5,6,1,2,3,4,5,6];
% pocet_vzoriek=length(f);
% cas=pocet_vzoriek;
% rozkmit=1;

% *** Doplnujuce vypisy *** %
fprintf('Tolerancia %.7f\n',max_rozdiel);
fprintf('Rozkmit funkcie %.2f\n',rozkmit);
fprintf('Pocet vzoriek %d\n',pocet_vzoriek);

subplot(2,1,1); %prvy graf
g1=stem(f);
title('Cela funkcia'); %titulok grafu
set(g1, 'Marker', 'none'); %koncove body nebudu zobrazene

max_rozdiel=(max_rozdiel*2); %tolerancia
sedi_perioda=1; %stav zistovania periody
najdena_perioda=0;

% v zadanom signale by nemali byt viac ako 400 period signalu
% perioda sa musi opekovat aspon 2x aby bola datekovatelna
% perioda musi byt vacsia ako nase 2 vzorky
% algoritmus nepocita s tym, zeby kvazi-periodicka funkcia bola na zaciatku neperiodicka
% *** Hlavny algoritmus *** %
for perioda=((pocet_vzoriek)/400):0.01:(floor((pocet_vzoriek-1)/2)),
    if (sedi_perioda==0)
        sedi_perioda=1;
    end;
    nasobna_perioda=floor(pocet_vzoriek/perioda)-1;
    while (sedi_perioda==1 && nasobna_perioda>=1)
        vzorka=floor(perioda*nasobna_perioda);
        i=0;
        while (sedi_perioda==1 && i<floor(perioda-1))
            debug_vzroka_1=f(1+i);
            debug_vzroka_2=f(i+vzorka);
            if (  abs( f(1+i)-f(i+vzorka)  )>max_rozdiel)
                sedi_perioda=0;
            end;
            i=i+1;
        end;
        nasobna_perioda=nasobna_perioda-1;
    end;
    if (sedi_perioda==1 && perioda>2)
        fprintf('Perioda %.7f\n',cas/pocet_vzoriek*(perioda-1));
        najdena_perioda=perioda;
        sedi_perioda=2;
    end;      
end;

% *** koncove vypisy *** %
if (najdena_perioda>0)
    tmp=(pocet_vzoriek/(najdena_perioda-1))-floor((pocet_vzoriek/(najdena_perioda-1)));
    if ( tmp <0.01 || tmp>0.99) %presnost 1 percento
        fprintf('-> Periodicka funkcia\n');
    else
        fprintf('-> Kvázy-periodicka funkcia\n');
    end;
    fprintf('Funkcia precnieva o %.3f percenta\n',tmp*100);
    
    prva_perioda=f(1:floor(najdena_perioda)-1);
    subplot(2,1,2); %druhy graf
    g2=stem (prva_perioda); 
    set(g2, 'Marker', 'none');
    title ('Prva perioda'); %nadpis grafu
else
    fprintf('-> Neperiodicka funkcia\n');
end;