clear all; close all; clc

%% načtení fasta souboru
[header seq] = fastaread('X64775.fasta'); % oryza sativa s mikrosatelity

% seq = randseq(60000); % vygenerování náhodné sekvence o délce 60000 nukleotidů
% seq(15:15:end) = 'A'; % každý 15. nukleotid je A
% seq(13:13:end) = 'T'; % každý 13. nukleotid je T
% seq(9:9:end) = 'G';   % každý 9. nukleotid je G
% seq(11:11:end) = 'C'; % každý 11. nukleotid je C

%% numerická reprezentace
L = length(seq); % délka sekvence

num_repr = zeros(4,length(seq)); 
num_repr(1,:)= seq == 'A'; % převod symbolického zápisu bází na numerickou formu
num_repr(2,:)= seq == 'T';
num_repr(3,:)= seq == 'C';
num_repr(4,:)= seq == 'G';

%% vytvoření spekrogramu
w = 120; % window size - délka okna
o = 1;  % overlap - posun oken
pom = 1;  % pomocná pro vkládání segmentů do sloupců vedle sebe ve spektrogramu

for u = 1:o:L-w % procházení numerickou reprezentací 'num_repr' po sloupcích
    for v = 1:4 % procházení numerickou reprezentací 'num_repr' po řádcích = každý pro jiný nukleotid
        x =(num_repr(v,u:u+w)); % procházíme vždy jeden řádek num. repr. ze všech čtyř pro všechny hodnoty sloupců 'u'
        spektrum = DFT(x); % volání funkce 'DFT.m'
        spektrum_celkem(v,:)= spektrum.^2; % čtyři spektra pro A,T,C,G uložíme do matice o čtyřech řádcích
    end
      
    s1 = spektrum_celkem(1,:);
    s2 = spektrum_celkem(2,:);
    s3 = spektrum_celkem(3,:);
    s4 = spektrum_celkem(4,:);
    
    % normalizace hodnot pixelů:
    Maxs1 = max(max(s1)); Mins1 = min(min(s1)); % normalizace hodnot na rozsah 0-1
    ns1 = s1 - Mins1; nns1 = ns1/(Maxs1-Mins1);
    
    Maxs2 = max(max(s2)); Mins2 = min(min(s2));
    ns2 = s2 - Mins2; nns2 = ns2/(Maxs2-Mins2);
    
    Maxs3 = max(max(s3)); Mins3 = min(min(s3));
    ns3 = s3 - Mins3; nns3 = ns3/(Maxs3-Mins3);
    
    Maxs4 = max(max(s4)); Mins4 = min(min(s4));
    ns4 = s4 - Mins4; nns4 = ns4/(Maxs4-Mins4);
    
    % mapování do RGB prostoru:
    tr = 1; tg = 0; tb = 0; % Red   [1 0 0] % barevné mapovací vektory
    cr = 0; cg = 1; cb = 0; % Green [0 1 0]
    ar = 0; ag = 0; ab = 1; % Blue  [0 0 1]
    gr = 0; gg = 0; gb = 0; % Black [0 0 0]
    
    R = ar.*nns1 + tr.*nns2 + cr.*nns3 + gr.*nns4;
    G = ag.*nns1 + tg.*nns2 + cg.*nns3 + gg.*nns4;
    B = ab.*nns1 + tb.*nns2 + cb.*nns3 + gb.*nns4;
    
    R=R'; % transponování vektorů
    G=G';
    B=B';
    
    % vkládání segmentů (sloupců) do spektrogramu:
    RGB(:,pom,1) = R;
    RGB(:,pom,2) = G;
    RGB(:,pom,3) = B;
                          
    pom = pom+1; % přechod na další segment
end

%% vykreslení spektrogramu
image(RGB)
periody = round((w./(1:w/2))*100)/100; % převedení svislé osy z frekvence na periodu
krok_osy_y = 3;
y_tick = [linspace(1,w/2,w/2)];
set(gca,'ytick',y_tick(1:krok_osy_y:end),'yTicklabel',periody(1:krok_osy_y:end)) 

[rows colms layers] = size(RGB); % upravení vodorovné osy zobrazující pozice bází
krok_osy_x = 10;
osa_x = fix(linspace(1,length(seq),krok_osy_x));
set(gca,'xtick',[linspace(1,colms,krok_osy_x)],'xTicklabel',osa_x)

title('Spektrogram sekvence X64775 kroky 1243') % název spektrogramu
xlabel('pozice v sekvenci')
ylabel('perioda')
