声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 5621|回复: 34

[综合] 请教apfft时移相位差法的问题

[复制链接]
发表于 2010-1-5 10:58 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
我想用apfft时移相位差法求实测数据序列的频率值,也就是信号的频率是未知的,按照您书中的原理,我是不是找到两段apfft后主谱线的相位后求其相位差,然后除以两序列的时延n0就可以了?但是n0该怎么求,不应该是对应的时间上的延迟吗,怎么书中又说到当n0=N时,就不用进行相位补偿,这个n0到底应该怎么计算啊?
谢谢了!

[ 本帖最后由 czk108 于 2010-1-5 15:43 编辑 ]

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2010-1-5 13:27 | 显示全部楼层

回复 楼主czk108的帖子

你想用apfft时移相位差法求实测数据序列的频率值.设频率值为k+dk,k为整数,dk为小於1的频偏值.

时移相位差法的两序列的时延一股都选为N,  所以频偏dk=(p2-p1)/(2*pi). 十分简单. < 数字信号全相位谱分析与滤波技术>附录3就给出了程序, 你一试就知道了.

如果你非要知道为什么n0=N时,频偏dk=(p2-p1)/(2*pi), 你可参见时移相位差法的论文

[ 本帖最后由 zhwang554 于 2010-1-5 15:30 编辑 ]
 楼主| 发表于 2010-1-5 14:55 | 显示全部楼层

回复 沙发 zhwang554 的帖子

您好,您所说的把频率值设为k+dk中的k是否是序列apfft谱的振幅最大值所对应的谱线号?
发表于 2010-1-5 15:02 | 显示全部楼层

回复 板凳czk108的帖子

对, k就是序列apfft谱的振幅最大值所对应的谱线号, p1和p2是频率k处的相位值
在程序中 k=rr=round(f1)

[ 本帖最后由 zhwang554 于 2010-1-5 16:56 编辑 ]
 楼主| 发表于 2010-1-5 16:12 | 显示全部楼层

回复 地板 zhwang554 的帖子

现在的问题是k值可以求取,而f1的值未知。我把程序中的f1改为f1=k+dk,dk=(p2-p1)/(2*pi)。把已知频率的正弦信号改为我的实测数据,引入采样率,其他程序不变。但是运行结果显示disp('频率校正值')
fff=floor(f1)+g(rr+1)  中g的个数超过了g本身的范围,不知道什么原因。  程序中的k值我是把最大值谱线号换算成对应的频率值的
发表于 2010-1-5 17:47 | 显示全部楼层

回复 6楼 czk108 的帖子

对实测数据, 最好自动找峰值,
你可以参考
whoru 发表于 2009-6-23 16:33   的贴子"求教:apFFT频谱校正问题!'的沙发中的程序
http://forum.vibunion.com/thread-83661-1-1.html
其中有判断dk>0.5或dk<0.5的语句,


你说:"fff=floor(f1)+g(rr+1)  中g的个数超过了g本身的范围,不知道什么原因。  程序中的k值我是把最大值谱线号换算成对应的频率值的"
你先从振幅谱中求出k值, 式中的g(rr+1)变为g(k), g的个数只有一个.
若有几个峰值k1,k2,...,分别代入, g的个数也是有限个.
 楼主| 发表于 2010-1-8 16:09 | 显示全部楼层

回复 6楼 zhwang554 的帖子

我在仿真的过程中发现一个问题:当我取信号频率f1为整数时,频率和振幅校正都会出现很大的偏差,然而取小数时校正效果非常理想,这是什么原因呢?有没有通用的程序对频率的取值没有限制的?
还有就是频率的取值是不是不能超过N的大小呢?
我的程序如下:
close all;clc;clear all;
N=1024;
t=-N+1:N*2-1;
f1=150;
A1=1;
ph1=60;
s=A1*cos(1*(2*pi*t*f1/N+ph1*pi/180));
% s=awgn(s,30);
win=hann(N)';
win1=hanning(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
w=pi*2;
s1=s(1:2*N-1);
y1=s1.*win2;%加窗后的信号
y1a=y1(N:end)+[0 y1(1:N-1)];%经全相预处理的N点序列

Out1=fft(y1a,N);
a1=abs(Out1);
p1=mod(angle(Out1),2*pi);
s2=s(1+N:3*N-1);
y2=s2.*win2;
y2a=y2(N:end)+[0 y2(1:N-1)];
Out2=fft(y2a,N);
a2=abs(Out2);
p2=mod(angle(Out2),2*pi);
g=mod((p2-p1)/pi/2,1);
h=2*pi*g.*(1-g.*g)./sin(pi*g);
aa1=abs((h.^2).*a2)/2;
rr=round(f1);
disp('频率校正值')
fff=floor(f1)+g(rr+1)
fid1=fopen('F:\仿真\apfftshiyi\频率.txt','a');  %存储路径
fprintf(fid1,'%.12f ',fff);
fclose(fid1);
disp('振幅校正值')
aaa=aa1(floor(f1)+1)
fid2=fopen('F:\仿真\apfftshiyi\振幅.txt','a');  %存储路径
fprintf(fid2,'%.12f ',aaa);
fclose(fid2);
disp('初相位校正值')
ppp=p1(rr+1)*180/pi
fid3=fopen('F:\仿真\apfftshiyi\初相位.txt','a');  %存储路径
fprintf(fid3,'%.12f ',ppp);
fclose(fid3);

[ 本帖最后由 czk108 于 2010-1-8 16:12 编辑 ]
发表于 2010-1-8 19:31 | 显示全部楼层

回复 楼上 czk108 的帖子

是振幅校正有问题
加一个判断语句可解决
参见 下网址

12楼中的程序
http://forum.vibunion.com/thread-77615-1-1.html
 楼主| 发表于 2010-1-8 23:47 | 显示全部楼层

回复 8楼 zhwang554 的帖子

我的程序中N=1024,当f1=150时远小于N/2,我试了下f1=50时,校正后为51,如果让f1取小数的话,哪怕是取f1=50.001,校正结果都很准确,就是不能取整数,很奇怪   这是为什么呢?
发表于 2010-1-9 04:30 | 显示全部楼层

回复 楼上czk108的帖子

频率f1<N即可,
问题出在当f1=150时
运行中g(rr+1)是1, 不是0
而g是对1取模的  g=mod((p2-p1)/pi/2,1) , 应是0
所以f1为整数时, 校正频率都多了1
为避开这个问题, 频偏取值也应在判断语句之后, 如下程序

close all;clc;clear all;
N=1024;
t=-N+1:N*2-1; f1=150.2;A1=1;ph1=30;
s=A1*cos(1*(2*pi*t*f1/N+ph1*pi/180));
win=hanning(N)';win1=hann(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
w=pi*2;
s1=s(1:2*N-1);%第1组(2N-1)个数据
y1=s1.*win2;
      y1a=y1(N:end)+[0 y1(1:N-1)];
Out1=fft(y1a,N);
a1=abs(Out1);
p1=mod(phase(Out1),2*pi);
s2=s(1+N:3*N-1); %第2组(2N-1)个数据
y2=s2.*win2;
y2a=y2(N:end)+[0 y2(1:N-1)];
Out2=fft(y2a,N);
a2=abs(Out2);
p2=mod(phase(Out2),2*pi);
dph=mod(p2-p1,2*pi);
rr=round(f1)
dph=dph(rr+1);
if dph>pi
dph=dph-2*pi;
elseif dph<-pi
dph=dph+2*pi;
end
g=dph/pi/2;
h=2*pi*g.*(1-g.*g)./sin(pi*g);
aa1=abs((h.^2).*a1)/2;
disp('频率校正值');
fff=rr+g
disp('振幅校正值');
aaa=aa1(rr+1)
disp('初相位校正值');
ppp=p1(rr+1)*180/pi

谢谢你发现了原程序中这个问题. 当时为了避免初学者弄不清楚, 在程序中不用判断语句

[ 本帖最后由 zhwang554 于 2010-1-9 09:24 编辑 ]
 楼主| 发表于 2010-1-12 10:38 | 显示全部楼层

回复 10楼 zhwang554 的帖子

我想在信号加噪声的情况下用蒙特卡罗模拟进行多次测量求均值,程序如下:
close all;clc;clear all;
for T=1:100
N=1024;
t=-N+1:N*2-1;
f1=155.325;
A1=1;
ph1=60;
s=A1*cos(1*(2*pi*t*f1/N+ph1*pi/180));
% s=awgn(s,15);
win=hann(N)';
win1=hanning(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
w=pi*2;
s1=s(1:2*N-1);
y1=s1.*win2;%加窗后的信号
y1a=y1(N:end)+[0 y1(1:N-1)];%经全相预处理的N点序列

Out1=fft(y1a,N);
a1=abs(Out1);
p1=mod(angle(Out1),2*pi);
s2=s(1+N:3*N-1);
y2=s2.*win2;
y2a=y2(N:end)+[0 y2(1:N-1)];
Out2=fft(y2a,N);
a2=abs(Out2);
p2=mod(angle(Out2),2*pi);
dph=mod(p2-p1,2*pi);
rr=round(f1);
dph=dph(rr+1);
if dph>pi
dph=dph-2*pi;
elseif dph<-pi
dph=dph+2*pi;
end
g=dph/pi/2;
h=2*pi*g.*(1-g.*g)./sin(pi*g);
aa1=abs((h.^2).*a2)/2;
end
disp('频率校正值')
fff=sum(rr+g)/100
disp('振幅校正值')
aaa=mean(aa1(floor(f1)+1))
disp('初相位校正值')
ppp=mean(p1(rr+1)*180/pi)
运行结果频率校正值为1.5532,在求频率校正值时应该是采用的T=100时的值,怎么才能使校正值为100次的仿真结果求均值呢,或者是蒙特卡罗模拟具体怎么用?谢谢
发表于 2010-1-13 07:03 | 显示全部楼层

回复 楼上 dzk108 的帖子

close all;clc;clear all;
T=100;k=1;
while T>0;
    T=T-1;
N=1024;
t=-N+1:N*2-1;
f1=155.325;
A1=1;
ph1=60;
s=A1*cos(1*(2*pi*t*f1/N+ph1*pi/180));
s=awgn(s,15);
win=hann(N)';
win1=hanning(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
w=pi*2;
s1=s(1:2*N-1);
y1=s1.*win2;%加窗后的信号
y1a=y1(N:end)+[0 y1(1:N-1)];%经全相预处理的N点序列
Out1=fft(y1a,N);
a1=abs(Out1);
p1=mod(angle(Out1),2*pi);
s2=s(1+N:3*N-1);
y2=s2.*win2;
y2a=y2(N:end)+[0 y2(1:N-1)];
Out2=fft(y2a,N);
a2=abs(Out2);
p2=mod(angle(Out2),2*pi);
dph=mod(p2-p1,2*pi);
rr=round(f1);
dph=dph(rr+1);
if dph>pi
dph=dph-2*pi;
elseif dph<-pi
dph=dph+2*pi;
end
g=dph/pi/2;
h=2*pi*g.*(1-g.*g)./sin(pi*g);
aa1=abs((h.^2).*a2)/2;
fff(k)=rr+g;
aaa(k)=aa1(floor(f1)+1);
ppp(k)=p1(rr+1)*180/pi;
k=k+1;
end
disp('频率校正值')
fff=sqrt(mean(fff.^2))
disp('振幅校正值')
aaa=sqrt(mean(aaa.^2))
disp('初相位校正值')
ppp=sqrt(mean(ppp.^2))


运行结果
频率校正值
fff =155.3247421237747
振幅校正值
aaa =1.00093021860210
初相位校正值
ppp =60.06201714402110

[ 本帖最后由 zhwang554 于 2010-1-13 08:08 编辑 ]
 楼主| 发表于 2010-1-14 10:50 | 显示全部楼层

回复 12楼 zhwang554 的帖子

您好!
我仿照您上面的程序,改了一下fft/apfft综合相位差法100次蒙特卡洛模拟的程序,如下,但是运行结果误差校大,不知道什么原因?
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=153.8;
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=mod(angle(y11_fft),w);
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=mod(angle(y2_fft),w);
ff=round(f1);
dph=mod(p2-p1,w);
dph=dph(ff+1);
if dph>pi
dph=dph-w;
elseif dph<-pi
dph=dph+w;
end
        ee=mod(dph/pi/(1-1/N),1);
        aa=(a1.^2)./a2*2;
        p22(k)=p2(ff+1);
        f22(k)=ee+floor(f1);
        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 =

    1.0465

frequency correction

f22 =

  153.1992

amplitude correction

aaa =

    1.0006

频率误差

df =

   -0.6008
谢谢!
发表于 2010-1-14 12:56 | 显示全部楼层

回复 楼上 czk108 的帖子

dph=mod(p2-p1,w);      改为   dph=mod(p1-p2,w);
p22(k)=p2(ff+1);           改为  p22(k)=p2(ff+1)*180/pi;
即可
 楼主| 发表于 2010-1-14 15:28 | 显示全部楼层

回复 14楼 zhwang554 的帖子

改完之后发现当f1取小数的时候精度很高,整数的时候频率误差大概有0.5Hz吧,由于以前不是学这个专业的无从分析啊,不知道是不是判断语句中的dph的范围取的不对, p22(k)=p2(ff+1);     跟 p22(k)=p2(ff+1)*180/pi;有什么区别?
您在求相位差时用的是比真实频率大的最小整数,为什么要这样选,是不是因为频率分辨率为1的原因呢,但是如果是的话,不明白为什么会是1.谢谢您的解答
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-26 03:58 , Processed in 0.106400 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表