There are several methods that give approximate results, but the only way I know to get exact results is the arctangent function, also available in complex math libraries as the carg (complex argument) function. Applying this to a complex sample gives you its phase, and since frequency is the derivative of phase you look at how much the phase changes from one sample to the next.
You can do this by taking the carg of each sample and subtracting the carg of the previous sample from it, but that requires some extra code to keep you in the range of -pi to +pi. Another way is to multiply the current sample by the complex conjugate of the previous sample and then take the carg of the result; that eliminates the need for a special test. This works because multiplying two complex numbers gives you a third complex number with phase equal to the sum of the two phase angles, so multiplying a complex number by the complex conjugate of another number gives you a complex number whose phase is equal to the difference between their phase angles, and that tells us the instantaneous frequency. Then you determine this actual angle with the carg() function.
The arctangent function is fairly complex, and at high sample rates it may be too much of a load on your CPU. But if you know that the peak FM deviation is a small fraction of the sample rate, i.e., you know that the phase won't change much from one sample to the next, then you can just use the small angle approximations sin(x) ~ x, cos(x) ~ 1 and tan(x) ~ x and get nearly the same result. So you can just do the multiply by the conjugate of the previous sample and look at the imaginary component of the result. As long as the samples are normalized (have magnitude 1) this will be the sine of the phase angle, which for small angles is approximately equal to the angle itself in radians.
Several methods apply this approximation, the most common being the cross product approach. You take the cross product of two consecutive samples, divide by their magnitudes (if they're not already normalized), and the magnitude of the result will be the change in phase. Again it's actually equal to the sine of the phase change, and for small angles this is approximately equal to the angle itself, in radians.
Hope this helps.
Reason: No reason