The signal spectrum only extends up to 4.04kHz, so 40kHz will be more than enough as explained by the Nyquist Theorem. Note that the sampling interval resolution does not limit the frequency resolution you can detect in the signal unless you are quantising it at 1-bit. As long as the quantisation errors do not swamp your signal, twice the highest frequency component plus some anti-aliasing filter margin is all you need. Anything higher is not actually adding any useful information.

Regarding your algorithm, it seems quite excessive to do an FFT simply to extract a narrowband frequency offset. If the signal frequency is stable enough I would look at having a couple of notch filters at 4.04kHz and 4.00kHz, calculating the output magnitude from them, and using the largest one to determine where the signal is. You can use a multi-rate filter to produce an extremely efficient algorithm.

Note that this is actually mathematically equivalent to doing a DFT with a window size equivalent to the filter width and then discarding all the coefficients except the two that correspond to the notch filter locations. As you have so much blank spectrum there is little point in calculating out all the other coefficients which is why a time domain solution should be more efficient.

The signal spectrum only extends up to 4.04kHz, so 40kHz will be more than enough as explained by the Nyquist Theorem. Note that the sampling interval resolution does not limit the frequency resolution you can detect in the signal unless you are quantising it at 1-bit. As long as the quantisation errors do not swamp your signal, twice the highest frequency component plus some anti-aliasing filter margin is all you need. Anything higher is not actually adding any useful information.

Regarding your algorithm, it seems quite excessive to do an FFT simply to extract a narrowband frequency offset. If the signal frequency is stable enough I would look at having a couple of notch filters at 4.04kHz and 4.00kHz, calculating the output magnitude from them, and using the largest one to determine where the signal is. You can use a multi-rate filter to produce an extremely efficient algorithm.

Note that this is actually mathematically equivalent to doing a DFT with up-sampling and windowing to give the required frequency discrimination, then discarding all the coefficients except the two that correspond to the notch filter locations. This is the algorithm you are essentially talking about doing, but as you have so much blank spectrum there is little point in calculating out all the other coefficients which is why a time domain solution should be more efficient.