Instantaneous frequency is the time derivative of phase. Since the OP already has the analytic representation of the chirp, every sample can easily be used for the frequency estimate (unlike zero crossing frequency estimators for real signals, where the better solution in that case is to use the Hilbert Transform to arrive at this point that the OP already has).    

The concern with using the derivative that was raised in the comments under the question is the discontinuity in phase. This is resolved by unwrapping the phase first and then use the derivative to compute the instantaneous frequency. Rick Lyons has a [nice blog post][1] on alternate discrete time derivative computations that can have better noise performance (considering the high frequency noise enhancement of a true derivative).

If using MATLAB, Octave or Python phase unwrapping is a built in function (np.unwrap(phase) in Python or unwrap(phase) in MATLAB/Octave.  Otherwise phase unwrapping can be done manually by detecting the discontinuity (easily done since the step change far exceeds any other sample to sample variation). 


Also note that for small phase steps and normalized amplitude, the frequency can be well approximated using cross product frequency discriminator due to the small angle approximation $\sin(\theta) \approx \theta$.  Note that $\sin(\theta)$ is the imaginary component for the conjugate product of two successive complex samples:

$$(I_1-jQ_1)(I_2+jQ_2)$$

$$= (I_1I_2+Q_1Q2) + j(I_1Q_2-Q_1I_2)$$

Given the frequency as a change in phase over a change in time, knowing the time step between the complex samples, if sufficiently oversampled (to maintain the small angle approximation) the frequency is given as:

$$f = \frac{\Delta \theta}{\Delta T} \approx f_s(I_1Q_2-Q_1I_2) $$


  [1]: https://www.dsprelated.com/showarticle/814.php