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';

%% mapování do RGB prostoru
% barevné mapovací vektory:
tr = 1; tg = 0; tb = 0; % Red   [1 0 0]
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.*(num_repr(1,:))+tr.*(num_repr(2,:))+cr.*(num_repr(3,:))+gr.*(num_repr(4,:)); % redukce 4D reprezentace do 3D
G = ag.*(num_repr(1,:))+tg.*(num_repr(2,:))+cg.*(num_repr(3,:))+gg.*(num_repr(4,:));
B = ab.*(num_repr(1,:))+tb.*(num_repr(2,:))+cb.*(num_repr(3,:))+gb.*(num_repr(4,:));

RGB(1,:)= R; % uložení tří vektorů do jedné matice
RGB(2,:)= G;
RGB(3,:)= B;

%% 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í maticí 'RGB' po sloupcích
    for v = 1:3 % procházení maticí 'RGB' po řádcích = každý pro jinou barvu
        x =(RGB(v,u:u+w)); % procházíme vždy jeden řádek 'RGB' ze všech tří pro všechny hodnoty sloupců 'u'
        spektrum = DFT(x); % volání funkce 'DFT.m'
        spektrum_celkem(v,:)= spektrum.^2; % tři spektra pro R,G,B uložíme do matice o třech řádcích
    end
    
    s1 = spektrum_celkem(1,:)'; % transponování vektorů
    s2 = spektrum_celkem(2,:)';
    s3 = spektrum_celkem(3,:)';
    
    % normalizace hodnot pixelů:
    Maxs1 = max(max(s1)); Mins1 = min(min(s1)); % normalizace hodnot RGB 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);
    
    % vkládání segmentů (sloupců) do spektrogramu:
    spektrogram(:,pom,1) = nns1;
    spektrogram(:,pom,2) = nns2;
    spektrogram(:,pom,3) = nns3;
                          
    pom = pom+1; % přechod na další segment
end

%% vykreslení spektrogramu
image(spektrogram)
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(spektrogram); % 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 1324') % název spektrogramu
xlabel('pozice v sekvenci')
ylabel('perioda')
