全相位FFT问题——全相位时移相位差法 请教！！！
最近学习王兆华教授的全相位FFT方面的知识，对于全相位时移相位差法有点疑问，望高手们解答。全相位时移相位差法的采样长度为3N-1，倘若将其分为三段：1~N、N+1~2N-1、2N~3N-1，来模拟一下实际采样中噪声干扰可能带来的偏差。
设定三段的频率、相位相同，但是幅值取不同值（保证了频谱的连续性），然后按照附录3的原理性程序做全相位FFT，并对最后的结果分析：
结果的频率和相位没什么好说的，主要是结果中幅值一项基本就等于前两段数据（1~N、N+1~2N-1）理论幅值的平均值，而第三段数据（2N~3N-1）幅值的变化对结果基本没什么大的影响。换句话说，全相位时移相位差法所得的结果基本等于前两段数据FFT的平均值，而几乎与第三段无关，这是为什么？？？而平常FFT做例如10个周波的FFT，那么结果一定等于各个周波FFT的平均值（已证明）。那么全相位第三段数据（2N~3N-1）的突然跳变会对全相位结果无太大的影响，求解释？！！
实验代码如下：
%附3 apFFT时移相位差法校正程序
close all;clear all;clc;
N=1024;
FSX=5;
fs=FSX*N;
t1=-N+1:0;
t2=1:N-1;
t3=N:2*N-1;
f=50;
ph=90;
A1=120;A2=60;A3=10000;%试图改变A1、A2、A3的值...
s1=A1*cos(2*pi*t1*f/fs+ph*pi/180);
s2=A2*cos(2*pi*t2*f/fs+ph*pi/180);
s3=A3*cos(2*pi*t3*f/fs+ph*pi/180);
ss=;
win=hanning(N)';win1=hann(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
w=pi*2;
s1=ss(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=ss(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);
rr=round(f/FSX);
dp=p2(rr+1)-p1(rr+1);
dph=mod(dp,2*pi); %由于频率较小，为减小噪声信号的干扰，作如下改进频率和幅值
if dph>pi %dph即为两序列谱峰线间的相位差，（0~2*pi）
dph=dph-2*pi;
elseif dph<-pi
dph=dph+2*pi;
end
g0=dph/pi/2;
if g0==0
h=2;
else h=2*pi*g0.*(1-g0.*g0)./sin(pi*g0);
end
AA=abs((h.^2).*a1)/2;
fff=(rr+g0)*FSX
aaa=AA(rr+1)
ppp=p1(rr+1)*180/pi
结果：fff =50.0251 aaa =90.0058 ppp =89.5443
即幅值aaa为前两段数据幅值（120、60）的平均值，改变其他数据同理。我想知道，第三段数据干什么去了？？？
本帖最后由 zhwang554 于 2011-10-8 06:15 编辑
回复 1 # wangze2008scu 的帖子
由桯序中第48句
AA=abs((h.^2).*a1)/2
可知,振幅校正值AA由a1算得;a1是由s1=ss(1:2*N-1)的apfft1,它只和1:2N-1个数据有关, 即只与A1,A2有关.
你的问题是A3干什么了?A3是第2个apfft所必须的,从apfft2测出相位p2用於计算频偏g0, 从而由g0插值出AA
若振幅校正值AA由a2算得,即桯序中第48句改为
AA=abs((h.^2).*a2)/2
结果：aaa =5031.56142755636
即幅值aaa近似为后两段数据幅值（60、10000）的平均值，
同样问题是A1干什么了?A1是第1个apfft所必须的,从apfft1测出相位p1用於计算频偏g0, 从而由g0插值出AA
FFT时移相位差法也一样, 若前段数据在采样时出现干扰（如幅值跳变），则对幅值结果影响较大；而若幅值跳变出现（不影响相位测量值）在后段，则对幅值影响较小。
回复 2 # zhwang554 的帖子
首先非常感谢王教授的解答，您说的方式之前我也试过，的确是您说的那样只与（2N-1）个数据有关，无论是前(2N-1)个还是后（2N-1）个。
那么是否可以这样理解，以之前源程序为例，全相位时移相位差法的结果(主要指幅值)很大程度上依靠前两段数据的值，而与第三段数据关系不大。换句话说，在实际采样过程中，若前两段数据在采样时出现干扰（如幅值跳变），则对结果影响较大；而若干扰出现在第三段，则影响较小？？? 本帖最后由 zhwang554 于 2011-10-7 22:29 编辑
回复 3 # wangze2008scu 的帖子
可以这样说：若前两段数据在采样时出现干扰（如幅值跳变），则对幅值结果影响较大；而若幅值跳变出现（不影响相位测量值）在第三段，则对幅值影响较小。
但不能说全相位时移相位差法的结果(主要指幅值)很大程度上依靠前两段数据的值，而与第三段数据关系不大。要点是第三段中的干扰是否影响相位测量值。第三段数据对幅值测量的影响反映在p2相位测量上。
回复 4 # zhwang554 的帖子
感谢教授细心的解答，以后肯定会有其他关于全相位的问题向您请教。
向老师致敬！ 回复 4 # zhwang554 的帖子
王老师，您好！再问一个问题哈：
我看你们著作《数字信号全相位谱分析与滤波技术》中附录2和附录3中程序均涉及到一个问题，在最终的频率、幅值、相位校正中都用了rr=round(f1)，即找出谱峰线的位置（rr+1），然后在此结果上进行校正，这样在做仿真研究的时候显然是可行的。但是若提前未知 f1数据时，难道方法会失效么？？
我的理解是，两个程序均建立在已知频率估计值时，若将仿真信号s替换为一组等效的实际数据，即本身不知道f1、A1、ph1参数时，如何使用全相位来进行三个参数的准确校正呢？
请老师不吝赐教，谢谢！
本帖最后由 zhwang554 于 2011-10-12 20:55 编辑
回复 6 # wangze2008scu 的帖子
附录2和附录3的2个程序都是用於理解fft/apfft和 apfft/apfft校正原理的程序, 不是实用程序, 所以
1)采样频率选为N, 即fft的频率分辨率为1Hz/格, 不作频率转换使初学者搞不清楚
2)不用相位差判断语句, 所以出现振幅校正公式中有floor语句
3)不加噪音
4)不加自动搜索振幅峰值语句, rr=round(f1)是为了找出谱峰线的位置（rr+1），
5)加已知余弦信号, 这样可以知道校正是否正确
采样频率选为fs, 加自动搜索振幅峰值和加相位差判断语句的 apfft/apfft 和 fft/apfft 校正程序,
见本论坛zhwang554日志
”自动搜索峰值的 apfft/apfft 和 fft/apfft 校正程序”
http://forum.vibunion.com/home-spa ... -blog-id-18060.html
这个程序可分析实际数据，也可对已知频率(已包恬振幅,相位)做估计
回复 7 # zhwang554 的帖子
万分感谢，好老师~~~:loveliness:
学习了。。
谢谢王老师！
本帖最后由 wangze2008scu 于 2012-4-22 12:54 编辑还有问题请教您！ 回复 7 # zhwang554 的帖子
王教授，您好。今天再借助此贴专程向您请教一下，关于全相位时移相位差法求取谐波相位的问题，望不吝赐教！具体阐述如下：
一般在工程应用中，多数是从已经采集好的数据样本中截取一段（长度满足3*N-1）进行数据分析。现在问题就来了，通过截取怎么求取谐波的相位（或者相对相位）？
现在举例如下——数据中含两个频率（50Hz、100Hz），相位分别为（60度、100度），模拟产生数据的程序如下：
close all;clear all;clc;
fs=512;
N=512;
FSX=fs/N;
%% ***模拟数据产生***
t=-N+1:10*N-1;
A1=1.0; f1=50; ph1=60;
A2=0.3; f2=100; ph2=100;
y=A1*cos(2*pi*t*f1/fs+ph1*pi/180)+A2*cos(2*pi*t*f2/fs+ph2*pi/180);
s=y(1:3*N-1); %假设设置的情况满足校正程序的要求，即中心样点过零
% s=y(1+117:3*N-1+117); % %假设另一种截取数据情况。
fid=fopen('s.txt','wt');
fprintf(fid, '%f\n',s);
fclose(fid);
disp('数据产生成功！');利用全相位时移相位差法，求取参数程序如下：
close all;clear all;clc;
fs=512;
N=512;
FSX=fs/N;
%% ***数据导入***
s=load('s.txt');
s=s';
%% ***全相位时移相位差校正法***
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(angle(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(angle(out2),2*pi);
%% ***其中频率f可以搜索峰值得到，此处直接设定为了简便！***
rr=;
%% ***校正程序***
for k=1:length(rr)
dp=p2(rr(k)+1)-p1(rr(k)+1);
dph=mod(dp,2*pi); %由于频率较小，为减小噪声信号的干扰，作如下改进频率和幅值
if dph>pi %dph即为两序列谱峰线间的相位差，（0~2*pi）
dph=dph-2*pi;
elseif dph<-pi
dph=dph+2*pi;
end
g0=dph/pi/2;
if g0==0
h=2;
else h=2*pi*g0.*(1-g0.*g0)./sin(pi*g0);
end
AA=abs((h.^2).*a1)/2;
fff(k)=(rr(k)+g0)*FSX;
aaa(k)=AA(rr(k)+1);
ppp(k)=p1(rr(k)+1)*180/pi;
fprintf('%f\t%f\t%f\n',fff(k),aaa(k),ppp(k));
end第一种情况：
s=y(1:3*N-1); %假设设置的情况满足校正程序的要求，即中心样点过零
结果：
50.000000 1.000000 60.000000
100.000000 0.300000 99.999994第二种情况： s=y(1+117:3*N-1+117); % %假设另一种截取数据情况。
结果： 50.000000 1.000000 213.281250
100.000000 0.300000 46.562494这样子，相位就不准确了。不过没关系，我们可以通过补偿校正相位，在ppp(k)=p1(rr(k)+1)*180/pi; 的下面添加如下语句：
ppp(k)=ppp(k)-mod(fff(k)/FSX*(117)/N*360,360);
ppp(k)=mod(ppp(k),360); 即可得到第一种情况的结果。
【问题来了】：
如此相位校正，是假设我们已知截取的数据段距离117点（相对t=-N+1:3N-1的情况）。在工程中，我们不知道截取的一段（3N-1长度）到底是哪一段，如何处理得到正确的相位（或者正确的相位差[如此处：40度]）？
希望老师解答一下，不胜感激！
本帖最后由 zhwang554 于 2012-4-22 21:48 编辑
回复 10 # wangze2008scu 的帖子
碥定一个正弦波,要3个参数:振幅,频率和相位,
这相位是任一己知时刻的相位,可以是初相是t=0的相位,或是距离117点的相位.
Apfft测的是2N-1数据点中间点的相位,即时间为取样点N的时刻的相位,这个相位的时刻知道,相位值也知道,这个正弦波就全确定了,它的波形全知道了,即任一点的相位都可求出,这有什么问题
你不是从117点算出初相(笫一个取样的相位).你也可以从其它已知时刻的相位(如apfft的笫N个取样)算出笫一个取样的相位,并算出任一其它对刻取样点的相位确定一个正弦波不一定要知初相,可以是任一己知时刻的相位, 回复 11 # zhwang554 的帖子
王老师，您好！首先感谢您的回复。
前面我可能描述得较细，我是想问如何从任意一段采样数据中求出这两个频率分量的相位ph1=60度、ph2=100度，或者求两者的相位差为40度？？ 回复 11 # zhwang554 的帖子
举个实例：
我这儿有一堆现在已经采样好的数据s.txt（见附件），已知采样频率fs=512Hz，余弦多次谐波信号。
怎么分离出个两个频率分量的相位ph1=60、ph2=100,或者两者的相位差为40度？
该如何处理呢？这个问题困扰我好久了，请王老解答一下，谢谢！
本帖最后由 zhwang554 于 2012-4-23 00:41 编辑
回复 13 # wangze2008scu 的帖子
s.txt中笫一个采样相位分别是60度和100度吗? 你如何做到的? 我算一下好样不是. 那一个采样相位是60度和100度.
相位差一般片指两个相同频率的任一时刻的相位差值,它们在任一时刻的相位差是相同的. 你说的两个不同频率(50Hz,100Hz)的40度初相相位差只在第一个取样存在,其它点相位差值是不同的
回复 14 # zhwang554 的帖子
首先非常感谢王老这么晚了还回复！
是这样子的：我就是按照10楼中“模拟数据产生”程序，其中s=y(1+151:3*N-1+151);产生数据文件s.txt，各频率参数和程序中给出的不变。
实际工程中就是这样子的：只能已知采样频率fs(如：10240Hz)，然后就是一串数据（类似s.txt中的数据），然后我用全相位时移相位差的方法进行求解，频率和幅值还好说，相位真心不知道应该怎么计算？一般算谐波相位（此处的2次谐波100Hz）都是要求相对于基波(50Hz)的相位，无论算出来基波相位是多少，此处的2次谐波相位=基波相位+40。
我就比较不太清楚，工程中的一段数据拿来分析时，往往就只知道采样频率，相位计算时到底以哪个为参照？