加窗FFT估计频率值的问题

第一曌

我的输入信号是一个频率随时间变化的正弦波，我需要实时测出它的频率，尝试过小波和HHT，效果都不理想。目前想到的方法是加窗FFT，我的采样频率是20k，一秒钟要分析出100个频率值，也就是说200个点需要出一个频率值。我的信号频率在1000赫兹上下波动，大致可以认为一个周期20个点，也就是说我的一个窗内包含10个周期的信号。因此，我进行仿真的时候，按照最理想的情况，一种频率的信号维持了十个周期，紧接着再是下一种频率的信号，按照时间顺接下去。为了提高分辨率，又进行了基于幅值比的插值算法，得出的结果还算符合预期，如图1。只是我再对频率值进行插值，也就是缩小频率变化的时候，效果就差的很多，如图2所示，想问一下有何改进方法，能再次提高频率分辨率？下面附了程序和数值，希望得到大家的帮助，谢谢~
clear all
close all
clc

sf=20000;
dx0=[980;984;988.1;992.3;996.4;1000.7;1002.6;1006.4;1008.4;1009.9;1011.5;1012.5;1013.2;1013.2;1013;1012.9;1010.5;1008.4;1006.2;1004.3;996.4;995.2;991.2;987.2;983.6;979.6;975.6;972.1;967.8;964;960.4;957;953.4;951.6;948.9;948.2;947;945.8;945.1;945;945.4;946.7;948.9;951.4;954.1;957.5;961.4;965.3;969.7;974.3;978.8;983.8];%dx0为频率值
dx3=0:51;
dx4=0:0.25:51;%对dx0四倍插值
dx5=interp1(dx3,dx0,dx4);
n=length(dx5);
dx1=10;%dx1为一个频率值对应的周期数
t1=cell(1,n);
x1=cell(1,n);
for dx2=1:n%dx2为频率值的个数
     t1{dx2}=0:1/sf:dx1/dx5(dx2);
     x1{dx2}=cos(2*pi*dx5(dx2)*t1{dx2});
end
C=cell2mat(x1);
N=length(C);
t=0:1/sf:(N-1)/sf;
data1=C';
figure
plot (t,data1)
title('输入信号');
xlabel('t/s');
ylabel('幅值');

length=200;%窗内点的个数
for(m=1:floor(size(data1,1)/length))
count = 1;
for(j=(m-1)*length+1:m*length)
    data2(count,1)=data1(j);
    count=count+1;
end
    x=data2';
    n=count-1;
   
y = fft(x);
Y=y(1:n/2)/n*2;
f = sf/2*linspace(0,1,n/2);
A=abs(Y);
phi=angle(Y);

%加矩形窗
[Amax,k]=max(A);
phmax=angle(Y(k));
if A(k-1)>A(k+1);
   deltf=1/(1+A(k)/A(k-1));
   f0=(k-1-deltf)*sf/n ;   %校正后频率
   am=Amax/sinc(deltf);     %校正后幅值
   ph=(phmax+pi*deltf)*180/pi ; %校正后相位
else A(k-1)<A(k+1);
     deltf=1/(1+A(k)/A(k+1));
     f0=(k-1+deltf)*sf/n;
     am=Amax/sinc(deltf);
     ph=(phmax-pi*deltf)*180/pi;
end
A(k)=am;f(k)=f0;
amp_rec=A(1:n/2);
max_where=find(amp_rec==max(amp_rec));
frequency0=f(max_where);
frequency1(m)=frequency0;
frequency2=frequency1';
end   
t=0:length/sf:(m-1)*length/sf;
figure
plot(t,frequency2,'b*');
title('时频曲线');
xlabel('t/s');
ylabel('f/Hz');

C:\Users\dx\Desktop\1
C:\Users\dx\Desktop\3

第一曌

呃。。。。第一次发，图怎么发啊？？囧。。。

第一曌

C:\Users\dx\Desktop\1