输出结果
fo =
50.0990
ao =
0.5000
pho =
127.0243
再来个9点的
clear all;clc
Fs=5000;
N=1024;
freq = 50.1;
amp1 = 0.5;
ph1= 127.0;
t = (0:N);
y = amp1*sin(2*pi*freq*t/Fs + ph1*pi/180.0);
xfft=fft(y,N);
y=abs(xfft);
=max(y(2:N/2));
z1=xfft(k)/576-(xfft(k-1)+xfft(k+1))/720+(xfft(k-2)+xfft(k+2))/1440-(xfft(k-3)+xfft(k+3))/5040+(xfft(k-4)+xfft(k+4))/40320;
k=k+1;
z2=xfft(k)/576-(xfft(k-1)+xfft(k+1))/720+(xfft(k-2)+xfft(k+2))/1440-(xfft(k-3)+xfft(k+3))/5040+(xfft(k-4)+xfft(k+4))/40320;
k=k-1;
ka=abs(z1)/abs(z2);
r=(5-4*ka)/(1+ka);
fo=(k-1+r)*Fs/N
ao=2*pi*r*(1-r*r)*(4-r*r)*(9-r*r)*(16-r*r)*(abs(z1))/(N*sin(r*pi))
pho=mod((angle(z1)-pi*r)*180/pi+90, 360)
结果：
fo =
50.1000
ao =
0.5000
pho =
127.0000
50.099999991413036
0.499999999774201
1.270000003234529e+02
看到精度没？
最近开始研究ZOOM-FFT，如何在不增加计算点数提高频谱分辨率，实现如1HZ或者更高分辨率的选频测量
clear all;clc
Fs=1000;
N=10000;
freq = 50.1;
amp1 = 0.5;
ph1= 127.0;
t = (0:N);
y = amp1*sin(2*pi*freq*t/Fs + ph1*pi/180.0);
xfft=fft(y,N);
y=abs(xfft);
=max(y(2:N/2));
k=k+1;
z1=xfft(k)-(xfft(k-1)+xfft(k+1))/2;
z2=xfft(k+1)-(xfft(k)+xfft(k+2))/2;
ka=abs(z1)/abs(z2);
r=(2-ka)/(1+ka);
fo=(k-1+r)*Fs/N
ao=2*pi*r*(1-r*r)*(abs(z1))/(N*sin(r*pi))
pho=(angle(z1)-pi*r)*180/pi+90
结果：
fo = 50.1000
ao = 0.5000
pho =127.0000
--------------------------
知道以上说明什么问题吗？
Fs=1000;
N=10000;
哎，我无语了！看来你真的没有实战经验！
如果freq = 123.45678呢？哎
回复 36 # dsp2008 的帖子
这个实现得必须配合锁相倍频电路来实现了
这样还校正个P，现实应用是不可能的
况且你这样10S才能出一个数，这表你能卖的出去吗？
回复 41 # a.gain 的帖子
你们一般要求多长时间，0.2s够不够？
回复 36 # dsp2008 的帖子
为什么要k=k+1？用max求的不就是最大值了嘛？你们有没密集型频谱校正的好方法？
仔细看求最大浦线的那行，把直流踢了，密集的只要频普分辩虑高点就行，频间隔有3根以上即可校正的很准确，这也是我研究zoom fft的原因
你可以研究现在数字示波器的模式，我们基本类似它，前端adc不变，改变数据抽取率
a.gain 发表于 2012-7-1 22:03 static/image/common/back.gif
Fs=1000;
N=10000;
哎，我无语了！看来你真的没有实战经验！
我只讨论“全相位”的理论谬误问题，不关心你的工程。