% FT pouze pro CpG ostrůvky a umělé rozptýlené repetice
clear all; close all; clc

%% načtení fasta souboru  
[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('umele_rozptylene_repetice_dlouhe.fasta'); % umělé rozptýlené repetice dlouhé 3x71 bp
% [header seq] = fastaread('umele_rozptylene_repetice.fasta'); % umělé rozptýlené repetice krátké 31x5 bp

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

num_repr = zeros(4,L); 
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í spektrogramu
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

% 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: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'
%         x = (rectwin(length(x1)))'.*x1; % vynásobení pravoúhlým oknem
        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
    
    % mapování do RGB prostoru:
    R1 = ar.*(spektrum_celkem(1,:))+tr.*(spektrum_celkem(2,:))+cr.*(spektrum_celkem(3,:))+gr.*(spektrum_celkem(4,:));
    G1 = ag.*(spektrum_celkem(1,:))+tg.*(spektrum_celkem(2,:))+cg.*(spektrum_celkem(3,:))+gg.*(spektrum_celkem(4,:));
    B1 = ab.*(spektrum_celkem(1,:))+tb.*(spektrum_celkem(2,:))+cb.*(spektrum_celkem(3,:))+gb.*(spektrum_celkem(4,:));
    
    R(:,pom)=R1'; % transponování vektorů
    G(:,pom)=G1';
    B(:,pom)=B1';
    
    pom = pom+1;
end

%% normalizace hodnot pixelů
% 1. normalizace
str_hodn = mean([mean2(R) mean2(G) mean2(B)]); % výpočet průměru středních hodnot z matic
R = R./(3*str_hodn); % str_hodn...normalizační konstanta pro snížení rozptylu hodnot v maticích
G = G./(3*str_hodn);
B = B./(3*str_hodn);
    
[m,n] = 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),n(i)); % hodnoty R>1 uložíme do prom. 'cervena'
    R(m(i),n(i))= 1; % nastavení R na max. hodnotu 1
    G(m(i),n(i))= G(m(i),n(i))/cervena; % snížení G o hodnoty R>1
    B(m(i),n(i))= B(m(i),n(i))/cervena; % snížení B o hodnoty R>1
end
    
[m,n] = find(G>1);
pocet = length(m);
for i = 1:pocet
    zelena = G(m(i),n(i)); % hodnoty G>1 uložíme do prom. 'zelena'
    G(m(i),n(i))= 1; % nastavení G na max. hodnotu 1
    R(m(i),n(i))= R(m(i),n(i))/zelena; % snížení R o hodnoty G>1
    B(m(i),n(i))= B(m(i),n(i))/zelena; % snížení B o hodnoty G>1
end
    
[m,n] = find(B>1);
pocet = length(m);
for i = 1:pocet
    modra = B(m(i),n(i)); % hodnoty B>1 uložíme do prom. 'modra'
    B(m(i),n(i))= 1; % nastavení B na max. hodnotu 1
    R(m(i),n(i))= R(m(i),n(i))/modra; % snížení R o hodnoty B>1
    G(m(i),n(i))= G(m(i),n(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(:,:,1) = nnR;
RGB(:,:,2) = nnG;
RGB(:,:,3) = nnB;
                          
%% vykreslení spektrogramu
image(RGB)
periody = round((w./(1:w/2))*100)/100; % převedení svislé osy z frekvence na periodu
krok_osy_y = 10;
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,L,krok_osy_x));
set(gca,'xtick',[linspace(1,colms,krok_osy_x)],'xTicklabel',osa_x)

title('human chr22') % název spektrogramu
xlabel('pozice v sekvenci')
ylabel('perioda')
