I was going to throw this out there to see if there are any comments or ideas.
The EU1KY analyzer uses the STM32F7 discovery board with an SAI audio codec. The audio codec provides 16 bits ADC sampling. To get the benefit of this sampling, however, one needs to know the phase of the intermediate frequency (approximately 10 kHz) to 2^-16 radians of accuracy. If one has a IF bandwidth say, of 50 Hz, that corresponds to 20 ms sampling time or (20 ms) X 2^-16 radians X 10 kHz = 0.03 radians or 1.7 degrees of tolerable phase drift of the intermediate frequency before the phase drift limits the accuracy of the demodulated signal rather than the bit sampling depth.
Now because one samples two signals, and one is only really interested in the difference in phase between those two signals, one should be able to use the phase drift of one to cancel the other. If one uses the FFT to demodulate of the signal, then one may be depending on the intermediate frequency to always land in a particular frequency bin. So if you have two signals you are trying to measure the phase difference between, you take the FFT of both of the signals, take the complex-valued amplitude and phase of each signal in the frequency bin you know where the demodulated signal is, and then multiple one signal's sample by the complex conjugate of the other, and you get the relative phase of the two. However, if there is phase noise, the signal's FFT is split "between" adjacent frequency bins, and does not land perfectly on the expected frequency bin, and so one does not perfectly cancel this phase noise. The nonlinear demodulation step here is the multiplication of the two frequency samples from the two FFT's.
Another possible approach is to do this demodulation step in the time domain rather than in the frequency domain. One can form the "complex analytic signal" versions of each of the two signals captured by the audio codec. This is done by using the Hilbert transform to infer the imaginary part of each signal from each real-valued signal. In practice, this may also be done by using the FFT of a real-valued signal. The real-valued signal has a FFT that is Hermitian conjugate, that is, the negative frequencies are the complex conjugate of the positive frequencies. If you take the FFT of a signal, zero out the negative frequencies, and take the inverse FFT, you get a complex signal that is the complex analytic signal version of the real signal. So the idea is to form the complex analytic signal version of each of the two signals, and in the time domain multiply the complex sample in one signal by the complex conjugate of the sample in the other signal, and then sum over all of the samples. By doing the multiplication in the time domain rather that in the frequency domain, any phase drift between the two signals is cancelled. One can get better noise immunity by, rather than just zeroing out negative frequencies, use a passband filter in the frequency domain only around those frequencies one knows the intermediate frequency should be in, but it is a larger bandwidth than just using one frequency bin.
I am not sure if this is already done in the EU1KY analyzer but I think that phase noise is a potential problem to obtain the highest accuracy, and there are separate crystals on the codec and the SI5351A and so some drift will occur between these. In the analyzer I designed, I actually generate a trigger signal using a comparator from the intermediate frequency which I use to sample the ADC which eliminates that kind of phase noise, but it may be possible to eliminate it through software as well. My design uses a comparatively low power processor without floating point and therefore is too slow to do real-time floating-point demodulation of the intermediate frequency signal.