How to perform carrier phase recovery in software?
To demodulate a phase-shift keyed signal, of which BPSK is the simplest, you have to recover the carrier frequency, phase, and symbol timing.

Bursty Signals
Some signals are bursty and provide a known data sequence called a preamble or mid-amble (depending on whether it shows up at the beginning or middle of the burst).  Demodulators can use a matched filter that "looks" for the known data sequence and use it to figure out the burst's frequency, phase, and symbol timing.  They do this for every burst and don't need to bother with "drift" (the gradual distance that builds up between the demodulator and the received signal as small errors accumulate) because the bursts are generally short enough that the drift isn't a problem.

Continuous Signals
Then there are continuous signals.  They are in many ways tougher to recover than bursty signals because you normally do not have a known data sequence to aid in locking up with the signal, and you do have to worry about drift even after you lock onto the signal.  I'll try to describe at a high level the main steps that are typically used to recover continuous signals.

Carrier Recovery

Usually you know what frequency the signal you are looking for will be at, or at least what set of frequencies it could be at.  Even with this knowledge, though, you usually have to be able to correct frequency offset because no two transmitters transmit at the same frequency.  There is always some error.  The usual method, then, is to mix the frequency that you think that the signal will be at, and then correct the residual frequency error.  This can be done with a Costas Loop, or by taking the fourth power of the baseband signal data and looking for a frequency spike.  There should be a frequency spike at the carrier offset * 4 (e.g. if you take an FFT of the data to the 4th power and see a frequency spike at 8300 Hz, that means the carrier offset is 8300/4 = 2075 Hz).  This is a very effective means of getting an initial lock on the frequency offset.  You can also use it to compensate for drift if you re-do it occasionally.  There is another way to compensate for drift which I will touch upon later.

Carrier Phase

At this point if you plotted your complex data in the complex plane (x-axis is real, y-axis is imaginary) it should look something like the following-



If you look closely you can see two dense areas towards the ends of the blurry line.  Those are the BPSK constellation points.  The points in between are the transitions between the constellation points.  Those will clear up once we get the symbol timing.  The reason that the line is at an angle is because of the carrier phase.  This can be measured by mirroring the signal by multiplying all points that have negative real values by $e^{j*\pi}$



and then taking the mean of the angles of the points.  Once you have calculated that, subtract that angle from all points that have had the carrier offset removed by multiplying the points by $e^{j*-\omega}$.  You can also compensate for carrier offset drift with this technique by dynamically updating the phase offset.  Once the data is phase corrected it should look something like this-



Once the data is phase corrected you can drop the imaginary part of the data because it does not add any information.

Symbol Timing

You should normally know, a priori, the symbol period of the signal you are trying to demodulate.  If you do need to determine the symbol period/frequency though, you can do it in a similar way to how the carrier offset was detected.  You can square the data which will cause a frequency spike at twice the symbol frequency.

Like with the carrier you will have to get the phase (timing) right and then compensate for drift.  The usual method for both of these problems is to look for the zero crossings.  Unless the noise is quite bad it should only cross the zero-point in the middle of a symbol transition from -1 to 1 or from 1 to -1.  Even if noise does cause this to happen in the middle of a symbol it won't happen very often.



The picture above is usually called an "eye diagram" or "eye pattern".  It is two symbol periods wide, and has many symbols "stacked" on top of one another.  I don't know if you are familiar with oscilloscopes or not, but you can get an oscilloscope to show a picture like this.  Anyway, the two "X's" are symbol transitions.  The high to low lines in the X's are when the symbol transitions from a 1 to a -1, and the low to high lines in the X's are when the symbol transitions from a -1 to a 1.  The point in between, where the pink line is, is the optimum place to sample the data to see if the symbol is a 1 or -1.

This same technique can be used to handle symbol timing drift.  Take a running average of the distance from the prior zero crossing and following zero crossing.  If the two averages are about the same, all is well.  If one is bigger than the other then you need to shift where you are grabbing your sample.

Once you have sampled the symbols at the correct points your data points should get something that looks like the classic BPSK constellation points.



I hope this helps.