I'm not sure what you mean by the three steps you describe in your post, especially what "shape into a NRZ signal" means or why you'd do it to something you just modulated as BPSK, or what your "sub-carrier" is a sub-carrier of or why it's at 64kHz, etc.   So I can't comment on whatever it is that you're doing because I don't understand it.

However, if you happen to have a BPSK signal that is timing-locked and phase-locked on the vertical axis of an I/Q plane, then "imag()" will recover the phase-locked signal and taking the sign of that (signum) will recover the raw bit stream.   angle() recovers the phase angle, which is not quite the same.

Hi, thanks for your reply.

Sorry if it was unclear. I tried to update my original post by removing redundant/useless information. I hope my issue would be understandable now :). 

sin(x) == x for x <<< 1  so its only used for very small values.  for large values it's incorrect.  you are multiplying e^j by h a scalar by a signal with large 180 degree phase changes (bpsk).  do you expect that you will get a large change in angle sometimes?  of course.  angle takes into account both real and imaginary - so you are removing the variability with taking only imag and not tan-1 (I/R).  So the first graph might be the most accurate?  And you are measuring the recovery of your phase modulated signal then the recovery of your bpsk signal.  why would it look like standard bpsk?