I don't think your formula is quite right.
The two complex samples are the locations (at two successive sampling instants)
of the tip of the rotating phasor that represents the analog signal. The information that you need is the angle between phasor positions at these
two successive time instants. If $S_{n+1} = r_{n+1}e^{j\theta_{n+1}}$ and
$S_n = r_ne^{j\theta_n}$, then you want the value of $\theta_{n+1}-\theta_n$.
Now, $S_{n+1}S_n^* = r_{n+1}r_ne^{j(\theta_{n+1}-\theta_n)}$ and if you
express these complex numbers in rectangular coordinates, then you can
get the desired angle as
$$\theta_{n+1}-\theta_n
= \arctan\left(\frac{\text{Im}(S_{n+1}S_n^*)}{\text{Re}(S_{n+1}S_n^*)}\right)
= \arctan\left(\frac{\text{Re}(S_{n})\text{Im}(S_{n+1})
- \text{Im}(S_n)\text{Re}(S_{n+1})}{\text{Re}(S_{n})\text{Re}(S_{n+1})
+ \text{Im}(S_{n})\text{Im}(S_{n+1})}\right)$$
though usually the four-quadrant arctangent function atan2(y,x) is used
instead of atan which returns an answer between $-\pi/2$ and $\pi/2$.
In digital communications applications such as DQPSK demodulation,
the real interest is not in the actual value of the angle, but in which
of the four intervals $(-\pi/4,\pi/4)$, $(\pi/4, 3\pi/4)$, $(3\pi/4, 5\pi/4)$,
and $(5\pi/4, 7\pi/4)$ the angle belongs -- in other words,
did the phase not change at all, or change by $\pi/2$, or by $\pi$, or
by $3\pi/2$ -- and this can be readily determined
by comparing the values of the numerator and denominator in
the argument of the arctangent above, thus saving computational
effort and speeding up the demodulation process. But the
more general formula gives the actual angle, and thus can be
used to demodulate a general FM signal as well. The output
will be a sequence of samples of the original continuous-time
signal that modulated the carrier and thus formed the signal
that was transmitted. D/A conversion will give the reconstructed
signal.
