zhwang554 发表于 2010-1-14 17:09

回复 楼上 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 编辑 ]

czk108 发表于 2010-1-14 20:38

回复 16楼 zhwang554 的帖子

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

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

zhwang554 发表于 2010-1-15 05:41

回复 楼上 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 编辑 ]

czk108 发表于 2010-1-15 09:46

回复 18楼 zhwang554 的帖子

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

zhwang554 发表于 2010-1-15 15:39

回复 楼上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 发表于 2010-1-21 11:35

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

zhwang554 发表于 2010-1-21 15:04

回复 楼上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 编辑 ]

dsfire 发表于 2010-1-22 11:36

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

czk108 发表于 2010-1-27 15:40

回复 20楼 zhwang554 的帖子

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

zhwang554 发表于 2010-1-27 16:48

回复 楼上 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 编辑 ]

czk108 发表于 2010-1-27 17:15

回复 25楼 zhwang554 的帖子

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

zhwang554 发表于 2010-1-27 18:00

回复 楼上 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 编辑 ]

czk108 发表于 2010-1-27 20:55

回复 27楼 zhwang554 的帖子

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

zhwang554 发表于 2010-1-27 21:17

回复 楼上 czk108的帖子

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

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

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

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

czk108 发表于 2010-1-27 23:54

回复 29楼 zhwang554 的帖子

可能是实测过程中存在噪声的问题吧,10组数据中都是存在最大值一个是19,一个20。
我在程序中把rr设为19,p1,p2都取的是19处的值,令g=(p2-p1)/(2*pi),最后频率fff=(19+g)*fs/N
这样求出的结果跟我想要的误差比较小,但是不太清楚这样做是否是正确的
页: 1 [2] 3
查看完整版本: 请教apfft时移相位差法的问题