% Demonstracia vyuzitia Matlabu na detekciu signalu v zasumenom
% signale. "Vstupny signal" je vytvoreny zmesou harmonickeho signalu s
% nahodnym posunutim a aditivneho sumu s gausovskym rozlozenim hustoty
% pravdepodobnosti.
% Detekcia je realizovana pomocou:
% a) "korelacneho" principu,
% b) pomocou konvolucie s vyuzitim tzv. prisposobeneho filtra.
%
% Hlavnym cielom je naznacenie moznosti vyuzitia specializovanych
% m-funkcii (konkretne funkcie filter) a tiez demonstracia suvislosti
% medzi korelaciou a konvoluciou.
% Zaroven je naznaceny sposob vyuzitia prace s externymi subormi.
%
% PPET, M.D. 20-03-2007, v.1.0

fvz = 10000;       % vzorkovacia frekvencia
fsig = 200;        % fekvencia "uzitocneho" harmonickeho signalu
Asig = 1.0;        % amplituda "uzitocneho" harmonickeho signalu
Asum = 2.0;        % parameter urcujuci amplitudu sumu (!!! nie je to
                   % priamo amplituda)

N = 1024;          % celkova dlzka (pocet vzoriek) analyzovaneho signalu
T = 100;           % dlzka (pocet vzoriek) uzitocneho signalu

% generovanie a zobrazenie "uzitocneho" signalu
sig = Asig*sin(2*pi*fsig*(0:T-1)/fvz);	% generovanie

figure(1);         % zobrazenie bude v samostatnom okne c.1
plot(sig);         % prikaz na zobrazenie
title('Uzitocny signal');               % popis obrazku
xlabel('Diskretny cas');
ylabel('Amplituda');
grid;              % zapnutie mriezky v zobrazenom obrazku

pause;             % docasne zastavenie vypoctu, caka sa Enter

% definovanie "nahodnej" polohy uzitocneho signalu v celkovom (N-vzorkovom)
% zazname
pol = round(1+ rand(1)*(N-T));	% poloha zaciatku
in = [ zeros(1,pol) sig ];      % "vlozenie" uzitocneho signalu (od pol+1)
in(N) = 0;                      % dodefinovanie celkovej dlzky zaznamu

figure(2);         % zobrazenie bude v samostatnom okne c.2
plot(in);          % prikaz na zobrazenie
title('Nahodne posunuty uzitocny signal');  % popis obrazku
xlabel('Diskretny cas');
ylabel('Amplituda');
grid;              % zapnutie mriezky v zobrazenom obrazku

pause;

% generovanie zasumeneho signalu (aditivna zmes in + Gaussov sum)
sum = Asum*randn(1,N); % vygenerovanie gausovskeho sumu
in_sum = in + sum; % aditivna zmes posunuteho signalu a sumu

figure(3);         % zobrazenie bude v samostatnom okne c.3
plot(in_sum);          % prikaz na zobrazenie
title('Zasumeny posunuty uzitocny signal');  % popis obrazku
xlabel('Diskretny cas');
ylabel('Amplituda');
grid;              % zapnutie mriezky v zobrazenom obrazku

pause;

% Vypocet KORELACIE medzi (vektormi) sig a in_sum. Vypocet je realizovany
% velmi neefektivne len s vyuzitim zakladnych prikazov Matlabu.
% Pocas prednasok a cviceni budu analyzovane efektivnejsie metody,
% ako je mozne tento vypocet v Matlabe realizovat.

% klasicky sposob vypoctu s vyuzitim KOREALCIE dvoch vektorov
out = zeros(1,N-T+1);  % nie je to nutne, ale alokacia moze zrychlit vypocet
for k=1:(N-T+1)      % cyklus v Matlabe
   out(k) = sig*in_sum(k:k+T-1)';   % vypocet korelacie pre posun "k"
end

% vypocet s vyuzitim konvolucie (pomocou prisposobeneho filtra)
N_sig = length(sig);
for k=1:N_sig      % impulzova odpoved prisposobeneho filtra je zrkadlovo
   h_imp(k) = sig(N_sig-k+1); % otocena verzia detekovaneho signalu
end
% filtracia (KONVOLUCIA!!!) zasumeneho signalu prisposobenym filtrom
% fumkcie filter standardnou funckiou Matlabu (je znacne optimalizovana)
in_sum2 = [in_sum zeros(1,T)];     % predlzenie signalu
out2 = filter(h_imp, 1, in_sum);
out2 = out2(N_sig:length(out2)); % odstranenie prechodoveho deja na zaciatku


figure(4);         % zobrazenie bude v samostatnom okne c.4
plot(out);         % prikaz na zobrazenie
title('Vystup korelacneho prijimaca');  % popis obrazku
xlabel('Posun');
ylabel('Amplituda');
grid;              % zapnutie mriezky v zobrazenom obrazku

pause;             % docasne zastavenie vypoctu, caka sa Enter

% zobrazenie s preddefinovanym skalovanim na osiach x a y
figure(5);         % zobrazenie bude v samostatnom okne c.5
plot(out-out2);    % prikaz na zobrazenie
title('Rozdiel medzi "konvolucnym" a "korelacnym" vypoctom');  % popis obrazku
xlabel('Posun');
ylabel('Chyba');
grid;              % zapnutie mriezky v zobrazenom obrazku

% zobrazenie s uzivatelom definovanym rozsahom v x-ovej suradnici
% (demonstruje vyuzitie prikazu axis)
roz = axis;        % rozlisenie v x-ovej a y-ovej suradnici obr. Figure(5)

figure(6);         % zobrazenie bude v samostatnom okne c.6
plot(out-out2);    % prikaz na zobrazenie
title('Rozdiel s definovanym rozsahom x-vej suradnice');  % popis obrazku
xlabel('Posun');
ylabel('Chyba');
axis([1 N-T+1 roz(3) roz(4)]);  % roz(3), roz(4) je rozlisenie v y-ovej suradnici
grid;              % zapnutie mriezky v zobrazenom obrazku

pause

[m,odhad] = max(out);  % najdenie maxima
pol = pol+1;       % realny zaciatok je posunuty o 1
pol                % presna poloha zaciatku (definovana pri generovani)
odhad              % odhad polohy pomocou korelacneho prijimaca

% demonstracie prace s externymi subormi. Filozofia je podobna ako v
% jazyku C (stdio kniznica), t.j. existuje podstatne vacsia flexibilita
% pri definovani vystupneho formatu ako to je pri vyuziti napr. prikazu
% save -ascii.

% v nasledujucej casti zapiseme vystup (po zaokruhleni) out ako cele cisla
% jednak v klasickom ASCII formate a tiez v binarnom formate (2 bajty/vzorku).
% Binarny format moze byt vyuzity napr. pri priamej komunikacii s technickymi
% prostreidkami (napr. DA prevodniky, mikrokontrolery a pod.).
x = round( out );   % zaokruhlenie cisel s plnou presnostou na cele cisla
save x.txt x -ascii % ulozenie vzoriek do suboru v ASCII formate (pre referencny vypocet)
outFile = fopen('x.dat','wb'); % vytvorenie identifikatora (handle) suboru (x.dat)
                               % pre zapis binarnych dat (wb)

fwrite(outFile,x,'short'); % ulozenie vzoriek v binarnom formate (2 bajty/vzorku)
fclose(outFile);

input('a teraz porovnajte subory x.dat a x.txt...');


