%% *************************************************************************************************
%  File: phase_shifting.m
% 
%  Brief: subor obsahuje algoritmus posuvania faze pre detekciu stredu 
%         interferencnych pruzkov
% 
%  Version: 1.0  Revision: 0.0   Name: peter novotny
% 
%  Date: 2011-01-23
%
%  Matlab version: R2009a
% **************************************************************************************************

% !! TENTO SUBOR NIEJE POUZITY KVOLI NEVHODNOSTI PRE PREVEDENE MERANIA !!

% clear all; close all; clc; 

I_psi = I_1D;

%% -------------------------------------------------------------------------------------------------
%  Automaticke zistenie hodnoty S(pocet vzoriek/pruzok)
I_temp = smoothts(I_psi,'b',11); % filter bude asi treba ladit(najlepsie zatial sedi WSIZE = 10)
[max_I_value, max_I_position] = max(I_temp);
[min_I_left, min_I_left_position] = min(I_temp(1:max_I_position));
[min_I_right, min_I_right_position] = min(I_temp(max_I_position:N));
S = (min_I_right_position) + (max_I_position - min_I_left_position);
% figure; plot(I_temp);
%---------------------------------------------------------------------------------------------------

%% *************************************************************************************************
%  5 step phase shift algorithm - podla T. Ailing
% pre simulaciu posunu zrkadla o 1/8 strednej vlnovej dlzky, ktora sposoby zmenu faze interferogramu
% o pi/2 vychadzam z toho, ze pocet vzorkov na jednu vlnovu dlzku interferogramu je S a preto 
% vyberiem kazdy S/4 vzorok co odpoveda OPD 1/4 lambda interferogramu (== fazy pi/2) a to je 
% rovne posunu zrkadla po kroku 1/8 lambda zdroja
% pi_pol - je mnozstvo vzoriek sposobujuce zmenu faze interferogramu o pi/2
% **************************************************************************************************
 
pi_pol = round(S/4);    

% predalokacia na radu matlabu kvoli rychlosti 
I_pi_pol = zeros(1,floor(N/pi_pol)); 

% vyber kazdeho n-teho vzorku odpovedajucemu posunutiu o pi/2 
i=1;k=1;
while i <= numel(I)
   if ~mod(i,pi_pol)
       I_pi_pol(k) = I(i);
       k = k + 1;
   end
   i = i + 1;
end
N_pi_pol = numel(I_pi_pol);

% figure, plot(I_pi_pol), grid

%%  Krok 1: Vypocet modulacie(viditelnosti) M a najdenie jej maximalnej hodnoty
%---------------------------------------------------------------------------------------------------

% zacina a konci sa s posunom o 3 vzorky kvoli pouzitiu 5 krokoveho algoritmu. M sa hlada pre I(k) 
% a je treba mat k dispozicii aj Ik-1 a Ik-2 ako aj Ik+1 a Ik+2 co sposobuje dany posuv
M_na_2 = zeros(1,N_pi_pol);
for k = 3:1:N_pi_pol-2            
    M_na_2(k) = ((I_pi_pol(k-1) - I_pi_pol(k+1))^2 - (I_pi_pol(k-2)...
                - I_pi_pol(k))*(I_pi_pol(k) - I_pi_pol(k+2)));
end

M_na_2_filt = smoothts(M_na_2,'b',5);

% figure, plot(M_na_2,'r'), grid, hold
% plot(M_na_2_filt,'g'),
% plot(I_pi_pol,'b'),

%%  Krok 2: Vypocet mikro-fazoveho posuvu delta_fi pre zaznam s maximalnou M
%---------------------------------------------------------------------------------------------------

[max_M_value, max_M_position] = max(M_na_2);
[max_M_filt_value, max_M_filt_position] = max(M_na_2_filt);

% nasledujuce priradenei sluzi len ako zmena znacenia kvoli zkrateniu zapisu nasledujuceho vzorca 
% a je treba zvolit bud filtrovany alebo nefiltrovany  stred
% k = max_M_position;         % stred z nefiltrovanej modulacie
k = max_M_filt_position;    % stred z filtrovanej modulacie

delta_fi_ailing = atan((2*I_pi_pol(k-1) - I_pi_pol(k+1))/...
                       (2*I_pi_pol(k) - I_pi_pol(k+2) - I_pi_pol(k-2)));
% delta_fi_larkin = atan(sqrt(((4*(I_pi_pol(k-1) - I_pi_pol(k+1))^2) - (I_pi_pol(k-2) - ...
%                   I_pi_pol(k+2))^2))/(2*I_pi_pol(k) - I_pi_pol(k+2) - I_pi_pol(k-2)));

% disp('delta fi ratana podla Ailing:'), disp(delta_fi_ailing)
% disp('delta fi ratana podla Larkina:'), disp(delta_fi_larkin)

%%  Krok 3: Najdenie stredu pruzkov 
%---------------------------------------------------------------------------------------------------

fringe_center_ailing = round(max_M_position*pi_pol + delta_fi_ailing*(S/(2*pi)));
% fringe_center_larkin = round(max_M_position*pi_pol + delta_fi_larkin*(S/(2*pi)));

fringe_center = z(round(fringe_center_ailing));

display('stred pruzkov je :');
disp(fringe_center);


% **************************************************************************************************
%  END OF FILE: phase_shifting.m
% **************************************************************************************************