FSK carrier frequency estimation
To perform carrier recovery for FSK, if the data is random and can assumed to be equiprobable, then you can take the mean of the derivative of the phase versus time to determine the carrier offset.

Optionally the phase change from one sample to the next can be estimated using using a complex conjugate product of successive samples (for small angles, the imaginary result of the complex conjugate product of two normalized vectors is approximately the phase). Since this is the change in phase over a time interval, this result will be directly proportional to frequency: 

$$f \propto Im(s(n)s(n-1)^*)$$

Where s(n) is the normalized complex FSK time domain signal.

Note in terms of I and Q samples this is simply:

$$I_{n-1}Q_n-I_nQ_{n-1}$$

The average of this result is proportional to the frequency of
the carrier and the deviation from average is the demodulated FSK signal. (The carrier component assuming a fixed offset can therefore be removed by using a simple DC Nulling Filter such as the one shown in this post: What does correcting IQ do?)

Since a small angle approximation was used, this provides a linear demodulation when the sampling rate is sufficiently high such that the phase change between two samples is a small angle such that $sin(\theta) \approx \theta$. For FSK that linearity is generally not as critical so would also work with much larger phase angles.