close all
clear all

fs=20000;       %vzorkovaci kmitocet 20kHz (dulezite)
diftype=[0.5 -0.5]*fs/pi;
gama=2.6752e8;  %gyromagneticka konstanta pro vodu
sirka=0.046;      %vzdalenost vrstev

%konstanta pro prepocet procent je 10/0.34
%tedy pokud v matlabu naměřím 0.05% z gradientu pak nastavuju 10/0.34*0.05=1.47% 

%načtení minus fidu
load bezpreGx0minusGy10kTR400msTAU15msf0Hz.mat
Gm10km=DATA;

%načtení plus fidu
load bezpreGx0plusGy10kTR400msTAU15msf0Hz.mat
Gm10kp=DATA;
[Nop,N]=size(Gm10kp);

%prah pro protišumovou filtraci, zpravidla volit pro plus i minus fid
%zvlášť
prahp=10;
prahm=10;

%Mazání vzorků Nmz na začátku a Nmk na koci fidu 
Nmz=20;
Nmk=50;

%nastavení parametrů sekvence pre_adj.pll
tg=0.5; %v ms, doba mezi koncem gradientu a prvni exitaci "graddeley" (nastavena +400us)
Top=15; %v ms, doba opakovani exitace a aqvizice, (nastavena +1.6ms)

Gin=70e-3; %T/m Velikost gradientu 10000 DAC
proc1=-0.02; %amplituda 1. exponenciály v procentech z Gin
proc2=0.0; %amplituda 2. exponenciály v procentech z Gin
T1=150;      %časová konstanta 1. exponenciály
T2=2000;     %časová konstanta 2. exponenciály
G0=-0.013e-3;


frekvp=[];
frekvm=[];

for k=1:Nop
    frekvp1(k,:)=filtIF(Gm10kp(k,:),prahp,fs,'db5',5);
    frekvm1(k,:)=filtIF(Gm10km(k,:),prahm,fs,'db5',5);
end

for k=1:Nop
    frekvp(k,:)=frekvp1(k,1:N);
    frekvm(k,:)=frekvm1(k,1:N);
end
    

for k=1:Nop
    grad(k,:)=((frekvp(k,:)-frekvm(k,:))/gama/sirka);
    B0(k,:)=(frekvp(k,:)+frekvm(k,:))/gama;
end



taq=linspace(0,(N-1)*1/fs,N)*1e3; %v ms

t=[];
for k=1:Nop
    t(k,:)=taq+Top*(k-1)+tg;
end

for k=1:Nop
    tm(k,:)=t(k,Nmz+1:end-Nmk);
    gradm(k,:)=grad(k,Nmz+1:end-Nmk)-G0;
    B0m(k,:)=B0(k,Nmz+1:end-Nmk);
end

figure(1)
for k=1:Nop
    plot(tm(k,:),gradm(k,:)*1000,'k')
    hold on
    title('grad')
    xlabel('t (ms)')
end

figure(2)
for k=1:Nop
    plot(tm(k,:),B0m(k,:),'k')
    hold on
    title('B0')
    xlabel('t (ms)')
end

tlin=reshape(tm',(N-Nmz-Nmk)*Nop,1)';
Glin=reshape(gradm',(N-Nmz-Nmk)*Nop,1)';
Blin=reshape(B0m',(N-Nmz-Nmk)*Nop,1)';


%ručne definovaný exponenciální pokles
expi=(Gin*proc1/100*exp(tlin./-T1)+Gin*proc2/100*exp(tlin./-T2));

