Digital modulation feature extraction using unwrapped phase

I created a BPSK signal in Matlab and I'm wanting to extract the unwrapped phase from the instantaneous phase obtained. The detail is that the result will apply an equation to get the centered non-linear component of the instantaneous phase which is given for "phi_NL" in the code below. This nonlinear phase should be between a value close to zero and pi, which is in accordance with the signal, as occurs in the BPSK phase shift of pi in each transition from 0 to 1 or from 1 to 0, in the encoding to be transmitted.

I've used the function unwrap Matlab, but the nonlinear phase ends infinitely increasing or decreasing. Researched more about the unwrapped phase in the internet and got the function below. However, the same result was obtained earlier. But changing the limits of "difference" to pi/2 and using a sine to create the modulated signal, the result is a phase that varies between zero and Pi, as expected.

I wonder if I can modify this threshold unwrap that way, since I can not get the expected result with the original code. Another thing I would like to know what is the relationship between this threshold with the cosine and sine used in signal modulation of the BPSK signal.

The full paper is in the address below:
http://www.radioeng.cz/fulltexts/2011/11_01_025_030.pdf

I found this code for unwrapping phase in this link below:
http://www.ljmu.ac.uk/GERI/CEORG_Doc...ping_Final.pdf

Thanks!

__________________________________________________ __________________________________

ns = 100; % Number of symbols
fs = 100*10^6; % Sample frequency
fc = 10*10^6; % Carrier frequency
rs = 500*10^3; % Symbol rate
nb = fs/rs; % Number of samples per symbol duration
Ts = 1/fs; % Sample interval


% Creating the modulated signal
% -------------------------------------------------
r = rand(ns,1); a = zeros(ns,1); n = zeros(ns*nb,1);

for i=1:ns
if (r(i)<0.5)
a(i) = 0;
else
a(i) = 1;
end
end

ind=0;
for i=1:1:ns
for j=nb*ind+1:nb*ind+nb
n(j) = a(i);
end
ind=ind+1;
end

z = zeros(ns*nb,1); % BPSK signal

for i=1:ns*nb
fase = n(i)*pi;
z(i) = sin(2*pi*fc*i*Ts + fase); % BPSK signal
end
% -------------------------------------------------

x = hilbert(z);
phi = angle(x); % Instantaneous phase
c = zeros(ns*nb,1);

% The linear component of the instantaneous phase
% -------------------------------------------------
for k = 1:ns*nb
c(k) = 2*pi*fc*k/fs;
end
% -------------------------------------------------


% Unwrapped phase function
% -------------------------------------------------
difference = zeros(ns*nb,1);
phi_uw = phi;

for i = 2:ns*nb
difference = phi(i) - phi(i-1);

if difference > pi/2
phi_uw(i:end) = phi_uw(i:end) - 2*pi;
end

if difference < -pi/2
phi_uw(i:end) = phi_uw(i:end) + 2*pi;
end
end
% --------------------------------------------------

% The centered non-linear component of the instantaneous phase
% --------------------------------------------------
phi_NL = phi_uw(1:1000) - c(1:1000);
% --------------------------------------------------

% PLOTS
figure; plot(n(1:1000));
figure; plot(phi_NL);
figure; plot(phi_uw(1:1000));