回复 楼上 czk108 的帖子
p22(k)=p2(ff+1) 是 弧度 p22(k)=p2(ff+1)*180/pi是 度. 程序中ph1=60 度, 所以用 度.频率取整数时误差为0.5还正不好介决, 到不是判断语句的范围取的不对, 无噪声p1-p2应是零,有噪时p1-p2可能是正一点点,也可能是负一点点, 当负一点点时输出就多了1
对fft/apfft法, 可不用判断语句, 因为fft的相位谱不是水平相位特性,频偏>0.5或<0.5相位不同, 自行判断了. 如下程序
close all;clc;clear all;
T=1;k=1;
while T>0;
T=T-1;
N=1024;
f1=155.72;
A1=1;
ph1=60;
w=2*pi;
t=-N+1:N-1;
y=A1*cos(1*(2*pi*t*f1/N+ph1*pi/180));
y=awgn(y,15);
y1 = y(N:2*N-1);%后N个输入数据
win =hanning(N)';;
win1 = win/sum(win);%窗归1
y11= y1.*win1;
y11_fft = fft(y11,N);
a1 = abs(y11_fft);%FFT振幅谱
p1 = mod(phase(y11_fft)*180/pi,360);;%FFT相位谱
y2 = y(1:2*N-1);%2N-1个输入数据
win =hanning(N)';;
winn =conv(win,win);%apFFT须要卷积窗
win2 = winn/sum(winn);%窗归1
y22= y2.*win2;
y222=y22(N:end)+;%构成长N的apFFT输入数据
y2_fft = fft(y222,N);;
a2 = abs(y2_fft);%apFFT振幅谱
p2=mod( phase(y2_fft)*180/pi,360);%apFFT相位谱
ee=(p1-p2)/180/(1-1/N);%频率偏离校正值
aa=(a1.^2)./a2*2;%振幅校正值
rr=round(f1);
p22(k)=p2(rr+1);
aaa(k)=aa(rr+1);
f22(k)=rr+ee(rr+1);
k=k+1;
end
disp('phase correction')
p22=sqrt(mean(p22.^2))
disp('frequency correction')
f22=sqrt(mean(f22.^2))
disp('amplitude correction')
aaa=sqrt(mean(aaa.^2))
你试用上程序,有否问题
ff=round(f1) 中round语句是4舍5入,ff+1点是振幅峰值点.+1是因为Matlab数值从1始,
当频率分辨率不为1时, 如采样频率为fs时, ff=round(f1*N/fs);
[ 本帖最后由 zhwang554 于 2010-1-15 03:22 编辑 ]
回复 16楼 zhwang554 的帖子
您好感谢您的指点您上面的程序我试了下 发现当小数部分落在0.5到0.7之间时误差会达到2Hz，其他区间的误差都很小达到10的-4次幂的级别，这是什么原因呢
[ 本帖最后由 czk108 于 2010-1-14 22:42 编辑 ]
回复 楼上 czk108 的帖子
本帖最后由 zhwang554 于 2010-8-23 14:15 编辑很抱歉, 16楼程序中
p1 = mod(phase(y11_fft)*180/pi,360); 改为 p1 = phase(y11_fft)*180/pi;
p2 = mod( phase(y2_fft)*180/pi,360); 改为 p2 = phase(y2_fft)*180/pi;
这样不会出现你所述的问题, 一取mod, 会将360或其倍数减去了, 影响频率1-2.Hz
你再试试, 有什么问题
[ 本帖最后由 zhwang554 于 2010-1-15 06:45 编辑 ]
回复 18楼 zhwang554 的帖子
您好 程序中p1，p2改完之后所有点校正结果都不准，误差大小为零点几到2点多，回复 楼上czk108 的帖子
喔!找到你在13楼上fft/apfft加判断语句程序在155.0时频率估计有时多1的原因了.ee=mod(dph/pi/(1-1/N),1); 改为 ee=dph/pi/(1-1/N);判断后不能再取mod(1),不然 -0.01 变成 0.99,频率估计就多1了
f22(k)=ee+floor(f1); 改为 f22(k)=ee+ff;
close all;clc;clear all;
T=100;k=1;
while T>0;
T=T-1;
N=1024;
w=2*pi;
t=-N+1:N-1;
f1=155.0;
ph1=60;
y=cos(2*pi*t*f1/N+ph1*pi/180);
y=awgn(y,15);
y1=y(N:2*N-1);
win=hanning(N)';
win1=win/sum(win);
y11=y1.*win1;
y11_fft=fft(y11,N);
a1=abs(y11_fft);
p1=angle(y11_fft);
y2=y(1:2*N-1);
win=hanning(N)';
winn=conv(win,win);%apFFT须要卷积窗
win2=winn/sum(winn);
y22=y2.*win2;
y222=y22(N:end)+;%构成长N的apFFT输入数据
y2_fft=fft(y222,N);
a2=abs(y2_fft);
p2=angle(y2_fft);
ff=round(f1);
dph=mod(p1-p2,w);
dph=dph(ff+1);
if dph>pi
dph=dph-w;
elseif dph<-pi
dph=dph+w;
end
ee=dph/pi/(1-1/N);
aa=(a1.^2)./a2*2;
p22(k)=p2(ff+1)*180/pi;
f22(k)=ee+ff;
a22(k)=aa(ff+1);
k=k+1;
end
disp('phase correction')
p22=sqrt(mean(p22.^2))
disp('frequency correction')
f22=sqrt(mean(f22.^2))
disp('amplitude correction')
aaa=sqrt(mean(a22.^2))
运行结果
phase correction
p22 =59.99897470440486
frequency correction
f22 =1.550001836767962e+002
amplitude correction
aaa =0.99785079144301
上程序中由於dph=mod(p1-p2,w)中,dph 已对w取mod, 前面p1和p2不必对w取mod;
至於16楼的不加判断的程序, 前面p1,p2不取mod后,后面ee=(p1-p2)/180/(1-1/N)必取mod. 不管它了.
你试一下上程序有何问题, 这程序在峰值取值,效果好
[ 本帖最后由 zhwang554 于 2010-1-15 16:21 编辑 ] 这些程序看了好几遍，感觉还是有点晕呼呼的！
回复 楼上dsfire 的帖子
本帖最后由 wdhd 于 2016-9-20 13:51 编辑最好先弄清<书>中附录2的fft/apfft谱分析程序, 弄清fft/apfft校正原理,
附录3的apfft/apfft谱分析程序,弄清apfft/apfft校正原理
它不加判断语句,循环语句,噪声,不搜索峰值,不会弄得晕呼呼.
若无<书>,可见下网址
http://forum.vibunion.com/UChome/space-62061-do-blog-id-17741.html
然后
10楼程序是apfft/apfft 加判断语句程序, 它在峰值读值,振幅校正效果好,无频偏时频率校正正确
12楼程序是apfft/apfft 加判断语句,加噪音及T=100次仿真结果
20楼程序是fft/apfft 加判断语句,加噪音及T=100次仿真结果
[ 本帖最后由 zhwang554 于 2010-1-21 22:43 编辑 ]
谢谢王老师的回复，您的书我已经拿到了，正在看呢！
回复 20楼 zhwang554 的帖子
您好，我在用时移相位差法对实测数据处理时，只是按照书中原理介绍的求出两序列apFFT谱后主谱线的相位值，然后做差再除以时延N得到信号的数字角频率w,然后根据模拟角频率和数字角频率的关系计算出信号频率。不知道这样做正确与否，但是这样做的结果不正确。另外，由apFFT谱得到的频率和直接由FFT谱得到的频率是否一样呢？
回复 楼上 czk108的帖子
本帖最后由 wdhd 于 2016-9-20 13:52 编辑由於相位360度循环, 测整个频率k+dk不易. 一般先测出频偏dk,
p1=phase(Out1)
p2=phase(Out2)
dk=mod(p2-p1)/(2*pi),1)
必对1取模
再估算频率整数值k
各种校正方法中最关键是测频偏, 由频偏可校正振幅, 在FFT中可校正相位
由apFFT谱得到的校正频率和直接由FFT谱得到的校正频率是一样的, 但精度有差别
[ 本帖最后由 zhwang554 于 2010-1-27 17:12 编辑 ]
回复 25楼 zhwang554 的帖子
由于受我的实测数据的影响，实际采样1536点，采样率为21900Hz，所以作谱点数最大只能取N=512了，所要估计的频率在800HZ左右，这样用您上面求dk的方法应该会超出它的范围，可不可以求频偏时就用两个谱最大值处的相位求，即p1=phase(A1)p2=phase(A2)
dk=mod(p2-p1)/(2*pi),1)
其中A1，A2为两个谱最大值的谱线号。
回复 楼上 czk108的帖子
本帖最后由 wdhd 于 2016-9-20 13:52 编辑不会超过的, 当N=512, fs=21900时, 频率800,4Hz的最大值谱线位置在rr=round(f1*N/fs)=19处, 所以频率800.4赫在N=512,fs=21900的频谱中没有超过范围., 你画出a1, 频率800.4Hz的峰值在19处.
所以在最大值处 19 求p1和p2, 得到的频率校正值再*fs/N.为 fff=(rr+g)*fs/N.
这里有频率转换关系
[ 本帖最后由 zhwang554 于 2010-1-27 18:50 编辑 ]
回复 27楼 zhwang554 的帖子
我画出的a1和a2最大值谱线号分别为20和19，应该取哪个呢？还有就是f1是未知的，是不是可以直接让rr=A1（也就是第一个序列的最大值谱线号）？是g=g=mod((p2-p1)/pi/2,1)吗？这样算出来的误差较大回复 楼上 czk108的帖子
a1,a2 图最大值都是20,什么会19,一个20? 是噪音引起一个19,一个20?图中最大值是20, 就取p1(20), p2（20）
rr算出耒=19, 最大值在rr+1=20.
无噪运行结果, 很好呀,是噪音引起误差大?
[ 本帖最后由 zhwang554 于 2010-1-27 21:52 编辑 ]
回复 29楼 zhwang554 的帖子
可能是实测过程中存在噪声的问题吧，10组数据中都是存在最大值一个是19，一个20。我在程序中把rr设为19，p1，p2都取的是19处的值，令g=(p2-p1)/(2*pi),最后频率fff=(19+g)*fs/N
这样求出的结果跟我想要的误差比较小，但是不太清楚这样做是否是正确的