Recover frequencies from IQ samples

I have a set of IQ data. From those data I'm trying to get the amplitudes and the frequencies of my signal as I then want to plot them vs. time.

I am able to obtain the amplitudes by squaring both my I and Q values, summing them up and taking the square root of the sum; however, I am struggling to obtain the frequency.

I understand that I need to take a FFT to get from the time domain to the frequency domain but I'm not sure on what values I should apply the FFT. Should I do it individually on Is and j*Qs and then sum them up. Should I do it on the sum Is + j*Qs (which just gave me a new array of complex numbers of the form x + j*y). What role my center frequency / sampling frequency are going to play into that?

For context: I'm doing this using Python. (And I'm obviously pretty new to all of this.)

Since you took 1000 point FFT, each of those complex output values correspond to frequencies at multiples of $$\frac{F_s}{1000}\text{Hz}$$ where $$F_s$$ is the Sampling Frequency. For example, the first value will be corresponding to frequency at $$0\text{Hz}$$, next value is the frequency content at $$\frac{F_s}{1000}\text{Hz}$$, 3rd value corresponds to frequency $$\frac{2F_s}{1000}\text{Hz}$$ and so on..
1. calculate the phase of your samples by simply performing $$\phi(k) =\tan^{-1}(s_I(k)/s_Q(k))$$.
2. Do not forget to unwrap $$\phi(k)$$ since it is limited in $$[-\pi/2,\pi/2)$$ interval. (Both Python and MATLAB have unwrap function.)
3. Take the derivative of $$\phi(k)$$ and divide to $$2\pi$$. You can use forward difference technique for this $$\frac{1}{2\pi}\frac{d\phi(k)}{dk} \approx \frac{ \phi(k+1)-\phi(k)}{2\pi}$$. The result will be your inst. frequency for your I/Q data.