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