% normalizace ve FOR cyklu
clear all; close all; clc

%% načtení fasta souboru 
% [header seq] = fastaread('X64775.fasta'); % oryza sativa s mikrosatelity 
[header seq] = fastaread('AL133493_human_chr21.fasta'); % human chr21 s minisatelity
% [header seq] = fastaread('randomseq_repetice_z_oryzy.fasta'); % náhodná sekvence s mikrosatelity z o. sativa
% [header seq] = fastaread('neisseria gonorrhoeae clone pNG273.fasta'); % rozptýlené repetice bakterie n. gonorrhoeae 
% [header seq] = fastaread('Rabbit short interspersed C repeat.fasta'); % rozptýlené repetice rabbit
% [header seq] = fastaread('Caenorhabditis elegans Cosmid F56F11.fasta'); % gen F56F11 C. elegans s pěti exony 
% [header seq] = fastaread('NT_011519.9 human chr22.fasta'); % CpG ostrůvky v chr22
% [header seq] = fastaread('Homo sapiens chromosome 21.fasta'); % CpG ostrůvky v chr21
% [header seq] = fastaread('Mycoplasma pneumoniae.fasta'); % část genomu 1-100000bp
% [header seq] = fastaread('C. carsonella ruddii.fasta'); % celý genom 160000bp
% [header seq] = fastaread('umele_rozptylene_repetice.fasta'); % umělé rozptýlené repetice

% 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

% seq(10000:3:10150) = 'C'; % trinukleotidová repetice CGG opakující se v délce 150 bp - mikrosatelit
% seq(10001:3:10151) = 'G';
% seq(10002:3:10152) = 'G';

%% 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 = 180; % window size - délka okna
o = 10;  % overlap - posun oken
pom = 1; % pomocná pro vkládání segmentů do sloupců vedle sebe ve spektrogramu
p = 50;  % řád AR modelu

% 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 = 1/3; gg = 1/3; gb = 1/3; % Grey  [1/3 1/3 1/3]

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)-mean(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 = p_burg(x,p); % volání funkce 'p_burg.m'
        spektrum_celkem(v,:)= spektrum.^2; % čtyři spektra pro A,T,C,G uložíme do matice o čtyřech řádcích
    end
    
   % mapování do RGB prostoru:
    R = ar.*(spektrum_celkem(1,:))+tr.*(spektrum_celkem(2,:))+cr.*(spektrum_celkem(3,:))+gr.*(spektrum_celkem(4,:));
    G = ag.*(spektrum_celkem(1,:))+tg.*(spektrum_celkem(2,:))+cg.*(spektrum_celkem(3,:))+gg.*(spektrum_celkem(4,:));
    B = ab.*(spektrum_celkem(1,:))+tb.*(spektrum_celkem(2,:))+cb.*(spektrum_celkem(3,:))+gb.*(spektrum_celkem(4,:));
    
    R=R'; % transponování vektorů
    G=G';
    B=B';
    
    % normalizace hodnot pixelů:       
    % 1. normalizace: 
    str_hodn = mean([mean(R) mean(G) mean(B)]); % výpočet průměru středních hodnot z vektorů
    R = R./(3*str_hodn); % str_hodn...normalizační konstanta pro snížení rozptylu hodnot ve vektorech
    G = G./(3*str_hodn);
    B = B./(3*str_hodn);
    
    [m] = find(R>1); % vyhledání hodnot přesahujících max. povolenou hodnotu 1
    pocet = length(m); % počet vyšších hodnot než 1
    for i = 1:pocet
        cervena = R(m(i)); % hodnoty R>1 uložíme do prom. 'cervena'
        R(m(i))= 1; % nastavení R na max. hodnotu 1
        G(m(i))= G(m(i))/cervena; % snížení G o hodnoty R>1
        B(m(i))= B(m(i))/cervena; % snížení B o hodnoty R>1
    end
    
    [m] = find(G>1);
    pocet = length(m);
    for i = 1:pocet
        zelena = G(m(i)); % hodnoty G>1 uložíme do prom. 'zelena'
        G(m(i))= 1; % nastavení G na max. hodnotu 1
        R(m(i))= R(m(i))/zelena; % snížení R o hodnoty G>1
        B(m(i))= B(m(i))/zelena; % snížení B o hodnoty G>1
    end
    
    [m] = find(B>1);
    pocet = length(m);
    for i = 1:pocet
        modra = B(m(i)); % hodnoty B>1 uložíme do prom. 'modra'
        B(m(i))= 1; % nastavení B na max. hodnotu 1
        R(m(i))= R(m(i))/modra; % snížení R o hodnoty B>1
        G(m(i))= G(m(i))/modra; % snížení G o hodnoty B>1
    end
    
    % 2. normalizace:
    MaxR = max(max(R)); MinR = min(min(R)); % normalizace hodnot RGB na rozsah 0-1
    nR = R - MinR; nnR = nR/(MaxR-MinR);
    
    MaxG = max(max(G)); MinG = min(min(G));
    nG = G - MinG; nnG = nG/(MaxG-MinG);
    
    MaxB = max(max(B)); MinB = min(min(B));
    nB = B - MinB; nnB = nB/(MaxB-MinB);
      
    % vkládání segmentů (sloupců) do spektrogramu:
    RGB(:,pom,1) = nnR; 
    RGB(:,pom,2) = nnG;
    RGB(:,pom,3) = nnB;
         
    pom = pom+1; % přechod na další segment
    
end

%% vykreslení spektrogramu
image(RGB)
periody=round((w./linspace(1,w/2,129))*100)/100; % převedení svislé osy z frekvence na periodu
krok_osy_y = 10;
y_tick = [linspace(1,129,129)];
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('AL133493 human chr21') % název spektrogramu
xlabel('pozice v sekvenci')
ylabel('perioda')                     
