How to find distorted signal frequency with CMSIS-DSP?

Sampled signal frequency calculation can be very simple but when signal is noisy or distorted it is an entirely different story.
Since there is more than one frequency and/or distortions, counting time between or number of zero-crossings is no longer a feasible solution.
The problem probably is generic but let's assume I want to measure main voltage 50/60Hz frequency and calculate RMS value - RMS is out of scope here.

I gather 4096 samples at about 3kHz every 1s and apply Hanning window. Using FFT transformation I correctly find a peek at 50Hz (33th bin), that is, at 1.6% of the full 1.5kHz scale. The problem is that the accuracy of such obtained measurement is limited since frequency bins are located every ~1.5Hz while I need at least 0.1Hz accuracy. Using weighted average of the adjacent bins does not seem to improve accuracy enough. To obtain the required precision this way, I would have to lower sample frequency to about 200Hz but then single measurement would take (1/200)*4096=20s instead of 0.6s - value I would rather not to exceed.

There must be a better way. Maybe once I know the dominant frequency, I should filter everything else out and then find zero-crossings? How do I go about it in practice?
Or maybe I should use phase-shifts returned by arm_rfft_q15() function as imaginary numbers at the beginning and at the end of the sampled interval? This would require running arm_rfft_q15() again using (e.g.) last 256 samples.

Is there any known, standard, robust solution? Please advise.