How to get smoothen 50 Hz frequency calculated from angle of the complex signal?

load('c.mat') % see attached file
fs=1200; 
ts=1/fs; 
t=0:ts:1-ts;
N=length(t);
f(1)=50;
for k=2:N
  a1(k)=angle(c(k));
  theta=unwrap(a1);
  f(k)=abs((theta(k)-theta(k-1))/(2*pi*ts));
end
plot(1:N,f)

The obtained frequency has very minute fluctuations due to taking discrete derivative or difference of angle at between two consecutive time steps.

My question is: How to get rid of these fluctuations to have a smooth frequency

I have already tried using smoothing filters; moving average, Savitzky Golay filters, IIR low pass filters, and rounding off the obtained frequency to one decimal point.

It would be great if any one can suggest me any useful technique to smooth this signal (50 Hz)