请教apfft时移相位差法的问题
我想用apfft时移相位差法求实测数据序列的频率值，也就是信号的频率是未知的，按照您书中的原理，我是不是找到两段apfft后主谱线的相位后求其相位差，然后除以两序列的时延n0就可以了？但是n0该怎么求，不应该是对应的时间上的延迟吗，怎么书中又说到当n0=N时，就不用进行相位补偿，这个n0到底应该怎么计算啊？谢谢了！
[ 本帖最后由 czk108 于 2010-1-5 15:43 编辑 ]
回复 楼主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 编辑 ]
回复 沙发 zhwang554 的帖子
您好，您所说的把频率值设为k+dk中的k是否是序列apfft谱的振幅最大值所对应的谱线号？回复 板凳czk108的帖子
对, k就是序列apfft谱的振幅最大值所对应的谱线号, p1和p2是频率k处的相位值在程序中 k=rr=round(f1)
[ 本帖最后由 zhwang554 于 2010-1-5 16:56 编辑 ]
回复 地板 zhwang554 的帖子
现在的问题是k值可以求取，而f1的值未知。我把程序中的f1改为f1=k+dk，dk=(p2-p1)/(2*pi)。把已知频率的正弦信号改为我的实测数据，引入采样率，其他程序不变。但是运行结果显示disp('频率校正值')fff=floor(f1)+g(rr+1)中g的个数超过了g本身的范围，不知道什么原因。程序中的k值我是把最大值谱线号换算成对应的频率值的
回复 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的个数也是有限个.
回复 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)+;%经全相预处理的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)+;
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 编辑 ]
回复 楼上 czk108 的帖子
是振幅校正有问题加一个判断语句可解决
参见 下网址
12楼中的程序
http://forum.vibunion.com/thread-77615-1-1.html
回复 8楼 zhwang554 的帖子
我的程序中N=1024，当f1=150时远小于N/2，我试了下f1=50时，校正后为51，如果让f1取小数的话，哪怕是取f1=50.001，校正结果都很准确，就是不能取整数，很奇怪 这是为什么呢？回复 楼上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)+;
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)+;
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 编辑 ]
回复 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)+;%经全相预处理的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)+;
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次的仿真结果求均值呢，或者是蒙特卡罗模拟具体怎么用？谢谢
回复 楼上 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)+;%经全相预处理的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)+;
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 编辑 ]
回复 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)+;%构成长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
谢谢！
回复 楼上 czk108 的帖子
dph=mod(p2-p1,w); 改为 dph=mod(p1-p2,w);p22(k)=p2(ff+1); 改为p22(k)=p2(ff+1)*180/pi;
即可
回复 14楼 zhwang554 的帖子
改完之后发现当f1取小数的时候精度很高，整数的时候频率误差大概有0.5Hz吧，由于以前不是学这个专业的无从分析啊，不知道是不是判断语句中的dph的范围取的不对， p22(k)=p2(ff+1); 跟 p22(k)=p2(ff+1)*180/pi;有什么区别？您在求相位差时用的是比真实频率大的最小整数，为什么要这样选，是不是因为频率分辨率为1的原因呢，但是如果是的话，不明白为什么会是1.谢谢您的解答