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

czk108

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

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

zhwang554

你想用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 编辑 ]

czk108

您好，您所说的把频率值设为k+dk中的k是否是序列apfft谱的振幅最大值所对应的谱线号？

zhwang554

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

[ 本帖最后由 zhwang554 于 2010-1-5 16:56 编辑 ]

czk108

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

zhwang554

对实测数据, 最好自动找峰值,
你可以参考
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的个数也是有限个.

czk108

我在仿真的过程中发现一个问题：当我取信号频率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 编辑 ]

zhwang554

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

12楼中的程序
http://forum.vibunion.com/thread-77615-1-1.html

czk108

我的程序中N=1024，当f1=150时远小于N/2，我试了下f1=50时，校正后为51，如果让f1取小数的话，哪怕是取f1=50.001，校正结果都很准确，就是不能取整数，很奇怪   这是为什么呢？

zhwang554

频率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 编辑 ]

czk108

我想在信号加噪声的情况下用蒙特卡罗模拟进行多次测量求均值，程序如下：
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次的仿真结果求均值呢，或者是蒙特卡罗模拟具体怎么用？谢谢

zhwang554

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 编辑 ]

czk108

您好！
我仿照您上面的程序，改了一下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
谢谢！

zhwang554

dph=mod(p2-p1,w);      改为   dph=mod(p1-p2,w);
p22(k)=p2(ff+1);           改为  p22(k)=p2(ff+1)*180/pi;
即可

czk108

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