请教apfft时移相位差法的问题

zhwang554

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)+[0 y22(1:N-1)];%构成长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 编辑 ]

czk108

您好  感谢您的指点
您上面的程序我试了下   发现当小数部分落在0.5到0.7之间时误差会达到2Hz，其他区间的误差都很小达到10的-4次幂的级别，这是什么原因呢

[ 本帖最后由 czk108 于 2010-1-14 22:42 编辑 ]

zhwang554

很抱歉, 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 编辑 ]

czk108

您好   程序中p1，p2改完之后所有点校正结果都不准，误差大小为零点几到2点多，

zhwang554

喔!找到你在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)+[0 y22(1:N-1)];%构成长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

这些程序看了好几遍，感觉还是有点晕呼呼的！

zhwang554

　　最好先弄清<书>中附录2的fft/apfft谱分析程序, 弄清fft/apfft校正原理,
　　附录3的apfft/apfft谱分析程序,弄清apfft/apfft校正原理
　　它不加判断语句,循环语句,噪声,不搜索峰值,不会弄得晕呼呼.
　　若无<书>,可见下网址
　　http://forum.vibunion.com/UChome ... -blog-id-17741.html
　　然后
　　10楼程序是apfft/apfft 加判断语句程序, 它在峰值读值,振幅校正效果好,无频偏时频率校正正确
　　12楼程序是apfft/apfft 加判断语句,加噪音及T=100次仿真结果
　　20楼程序是fft/apfft 加判断语句,加噪音及T=100次仿真结果
　　[ 本帖最后由 zhwang554 于 2010-1-21 22:43 编辑 ]

dsfire

谢谢王老师的回复，您的书我已经拿到了，正在看呢！

czk108

您好，我在用时移相位差法对实测数据处理时，只是按照书中原理介绍的求出两序列apFFT谱后主谱线的相位值，然后做差再除以时延N得到信号的数字角频率w,然后根据模拟角频率和数字角频率的关系计算出信号频率。不知道这样做正确与否，但是这样做的结果不正确。
另外，由apFFT谱得到的频率和直接由FFT谱得到的频率是否一样呢？

zhwang554

　　由於相位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 编辑 ]

czk108

由于受我的实测数据的影响，实际采样1536点，采样率为21900Hz，所以作谱点数最大只能取N=512了，所要估计的频率在800HZ左右，这样用您上面求dk的方法应该会超出它的范围，可不可以求频偏时就用两个谱最大值处的相位求，即p1=phase(A1)
p2=phase(A2)
dk=mod(p2-p1)/(2*pi),1)
其中A1，A2为两个谱最大值的谱线号。

zhwang554

　　不会超过的, 当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 编辑 ]

czk108

我画出的a1和a2最大值谱线号分别为20和19，应该取哪个呢？还有就是f1是未知的，是不是可以直接让rr=A1（也就是第一个序列的最大值谱线号）？是g=g=mod((p2-p1)/pi/2,1)吗？这样算出来的误差较大

zhwang554

a1,a2 图最大值都是20,什么会19,一个20? 是噪音引起一个19,一个20?
图中最大值是20, 就取p1(20), p2（20）

rr算出耒=19, 最大值在rr+1=20.

无噪运行结果, 很好呀,是噪音引起误差大?

[ 本帖最后由 zhwang554 于 2010-1-27 21:52 编辑 ]

czk108

可能是实测过程中存在噪声的问题吧，10组数据中都是存在最大值一个是19，一个20。
我在程序中把rr设为19，p1，p2都取的是19处的值，令g=(p2-p1)/(2*pi),最后频率fff=(19+g)*fs/N
这样求出的结果跟我想要的误差比较小，但是不太清楚这样做是否是正确的