One possible problem with converting to a square wave first is that the effect of the noise may actually be worse, especially if your conversion to a square wave involves comparators with hysteresis. If you have any hysteresis in your comparators then a very short noise pulse that would have had a small effect on the FFT of the sine wave, will instead cause a longer-lasting level change and therefore show up more prominently in the frequency spectrum. But frankly I think you are doing about as well as can be expected, getting within 0.213 Hz with only 50 msec. samples. Are you sure there is no way you can cheat? If the requirement is only to produce an updated answer every 50 msec., you might not need to limit yourself to FFTs covering only 50 msec. You might form a FIFO of samples, perhaps 500 msec. long. With each new 50 msec. sample, add that sample on to the end of the FIFO and drop the first 50 msec. Then FFT the whole thing every 50 msec. That will give you an answer every 50 msec. (assuming the signal is continuous). The only trouble is when the frequency does change it will take your FIFO a full 500 msec. to fully switch over to the new frequency. In the interim there will be calculations that gradually transition from the old frequency to the new frequency. But if this signal just appears out of nowhere, lasts only 50 msec., the disappears, then you are stuck. No cheating is possible.
I noticed you said you used parabolic interpolation. I used that on my musical instrument tuner app initially. But I found that I could do better by constructing a look-up table post-processor to the parabolic interpolation. I developed the lookup table by generating signals of various frequencies in software and then applying parabolic interpolation to them. Since I generated the signals in software, I know exactly what their frequencies are. I divided the range of +/- 1 frequency bin up into 64 table entries and then just searched the table to find the range that best fit any particular parabolic interpolation. Also, things like padding with zeroes and windowing the time series before the FFT do have an effect on the accuracy of parabolic interpolation, so if you are going to develop a lookup table, make sure you develop it using the same padding and windowing as you plan to use on the real signals. FWIW, I tried windowing, but I don't do it anymore, as it seems to affect only the low frequency end of the FFT.