求实测信号FFT之后得到的正确幅值和相位,很困惑啊

期待着一天

file:///C:/Documents%20and%20Settings/Administrator/桌面我用传感器实测了一组数据，转速是6000转/分，FFT之后应该得到信号的频率是100HZ，因为我的采样频率是20000HZ，采样数一般为2的整数幂，所以就使信号不是整周期采样，要想得到正确的频率，幅值和相位，必须通过频率校正，我在下面的程序中用的是比值校正法，但是程序运行的结果，频率只有97和90多一点，并且校正之后还不如不校正的好，我不知道程序那里不对还是这个频谱校正的方法不适合？下面是程序，希望大侠们仔细阅读提出修改，不胜感激。数据在附件里） 1024x.txt (66.02 KB, 下载次数: 8)

2011-8-29 14:45 上传

点击文件名下载附件

function bizhi1()
close all;clear all;clc  
Fs=20000;N=2048;
t =(0:N-1)/Fs;
tt=0:N-1;
hanning=0.5-0.5*cos(2*pi*tt/N);%汉宁窗的表达式
x=load('1024x.txt');
x=x(1:N)
x_hanning=x'.*hanning;
y = fft(x);
yh=fft(x_hanning);
Y=y(1:N/2)/N*2;
Yh=yh(1:N/2)/N*2;
f = Fs/2*linspace(0,1,N/2);
A=abs(Y);
Ah=abs(Yh);
subplot(221);stem(f,A(1:N/2));title('加矩形窗未校正');
subplot(222);stem(f,2*Ah(1:N/2));title('加汉宁窗未校正');%作图函数中的2倍是幅值恢复系数
%加矩形窗
[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)*Fs/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)*Fs/N
    am=Amax/sinc(deltf)
    ph=(phmax-pi*deltf)*180/pi
end
A(k)=am;f(k)=f0;
subplot(223);stem(f,A(1:N/2));title('加矩形窗校正后');
%加汉宁窗
[Amaxh,kh]=max(Ah);
phmaxh=angle(Yh(kh));
if Ah(kh-1)>Ah(kh+1);
    deltfh=(Ah(kh)/Ah(kh-1)-2)/(1+Ah(kh)/Ah(kh-1));
    f0h=(kh-1+deltfh)*Fs/N
    amh=2/sinc(deltfh)*Amaxh*(1-deltfh^2)
    phh=(phmaxh-pi*deltfh)*180/pi
else Ah(kh-1)<Ah(kh+1);
    deltfh=(Ah(kh)/Ah(kh+1)-2)/(1+Ah(kh)/Ah(kh+1));
    f0h=(kh-1-deltfh)*Fs/N
    amh=2/sinc(deltfh)*Amaxh*(1-deltfh^2)
    phh=(phmaxh+pi*deltfh)*180/pi
end
Ah(kh)=amh;f(kh)=f0h;
subplot(224);stem(f,Ah(1:N/2));title('加汉宁窗校正后');

f0 =
   96.1969
am =
    0.5659
ph =
  116.4117
f0h =
   96.0570
amh =
    0.6109
phh =
  121.8247