全相位FFT问题——全相位时移相位差法 请教！！！

wangze2008scu

       最近学习王兆华教授的全相位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）的突然跳变会对全相位结果无太大的影响，求解释？！！
实验代码如下：

1. 
   
2. %附3 apFFT时移相位差法校正程序
   
3. close all;clear all;clc;
   
4. N=1024;
   
5. FSX=5;
   
6. fs=FSX*N;
   
7. t1=-N+1:0;
   
8. t2=1:N-1;
   
9. t3=N:2*N-1;
   
10. f=50;
    
11. ph=90;
    
12. A1=120;A2=60;A3=10000;%试图改变A1、A2、A3的值...
    
13. s1=A1*cos(2*pi*t1*f/fs+ph*pi/180);
    
14. s2=A2*cos(2*pi*t2*f/fs+ph*pi/180);
    
15. s3=A3*cos(2*pi*t3*f/fs+ph*pi/180);
    
16. ss=[s1,s2,s3];
    
17. 
    
18. win=hanning(N)';win1=hann(N)';
    
19. win2=conv(win,win1);
    
20. win2=win2/sum(win2);
    
21. w=pi*2;
    
22. s1=ss(1:2*N-1); %第1组（2N-1）个数据
    
23. y1=s1.*win2;
    
24. y1a=y1(N:end)+[0 y1(1:N-1)];
    
25. out1=fft(y1a,N);
    
26. a1=abs(out1);
    
27. p1=mod(phase(out1),2*pi);
    
28. s2=ss(1+N:3*N-1); %第2组（2N-1）个数据
    
29. y2=s2.*win2;
    
30. y2a=y2(N:end)+[0 y2(1:N-1)];
    
31. out2=fft(y2a,N);
    
32. a2=abs(out2);
    
33. p2=mod(phase(out2),2*pi);
    
34. 
    
35. rr=round(f/FSX);
    
36. dp=p2(rr+1)-p1(rr+1);
    
37. dph=mod(dp,2*pi); %由于频率较小，为减小噪声信号的干扰，作如下改进频率和幅值
    
38. if dph>pi %dph即为两序列谱峰线间的相位差，（0~2*pi）
    
39. dph=dph-2*pi;
    
40. elseif dph<-pi
    
41. dph=dph+2*pi;
    
42. end
    
43. g0=dph/pi/2;
    
44. if g0==0
    
45. h=2;
    
46. else h=2*pi*g0.*(1-g0.*g0)./sin(pi*g0);
    
47. end
    
48. 
    
49. AA=abs((h.^2).*a1)/2;
    
50. fff=(rr+g0)*FSX
    
51. aaa=AA(rr+1)
    
52. ppp=p1(rr+1)*180/pi
    
53. 
    

复制代码

结果：fff =  50.0251        aaa =  90.0058           ppp =  89.5443
即幅值aaa为前两段数据幅值（120、60）的平均值，改变其他数据同理。我想知道，第三段数据干什么去了？？？

zhwang554

回复 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时移相位差法也一样, 若前段数据在采样时出现干扰（如幅值跳变），则对幅值结果影响较大；而若幅值跳变出现（不影响相位测量值）在后段，则对幅值影响较小。

wangze2008scu

回复 2 # zhwang554 的帖子

首先非常感谢王教授的解答，您说的方式之前我也试过，的确是您说的那样只与（2N-1）个数据有关，无论是前(2N-1)个还是后（2N-1）个。
那么是否可以这样理解，以之前源程序为例，全相位时移相位差法的结果(主要指幅值)很大程度上依靠前两段数据的值，而与第三段数据关系不大。换句话说，在实际采样过程中，若前两段数据在采样时出现干扰（如幅值跳变），则对结果影响较大；而若干扰出现在第三段，则影响较小？？?

zhwang554

回复 3 # wangze2008scu 的帖子

可以这样说：若前两段数据在采样时出现干扰（如幅值跳变），则对幅值结果影响较大；而若幅值跳变出现（不影响相位测量值）在第三段，则对幅值影响较小。
但不能说全相位时移相位差法的结果(主要指幅值)很大程度上依靠前两段数据的值，而与第三段数据关系不大。要点是第三段中的干扰是否影响相位测量值。第三段数据对幅值测量的影响反映在p2相位测量上。

wangze2008scu

回复 4 # zhwang554 的帖子

感谢教授细心的解答，以后肯定会有其他关于全相位的问题向您请教。
向老师致敬！

wangze2008scu

回复 4 # zhwang554 的帖子

王老师，您好！再问一个问题哈：
       我看你们著作《数字信号全相位谱分析与滤波技术》中附录2和附录3中程序均涉及到一个问题，在最终的频率、幅值、相位校正中都用了rr=round(f1)，即找出谱峰线的位置（rr+1），然后在此结果上进行校正，这样在做仿真研究的时候显然是可行的。但是若提前未知 f1数据时，难道方法会失效么？？
       我的理解是，两个程序均建立在已知频率估计值时，若将仿真信号s替换为一组等效的实际数据，即本身不知道f1、A1、ph1参数时，如何使用全相位来进行三个参数的准确校正呢？
       请老师不吝赐教，谢谢！

zhwang554

回复 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
这个程序可分析实际数据，也可对已知频率(已包恬振幅,相位)做估计

wangze2008scu

回复 7 # zhwang554 的帖子

万分感谢，好老师~~~:loveliness:
学习了。。

wangze2008scu

还有问题请教您！

wangze2008scu

回复 7 # zhwang554 的帖子

王教授，您好。今天再借助此贴专程向您请教一下，关于全相位时移相位差法求取谐波相位的问题，望不吝赐教！具体阐述如下：
        一般在工程应用中，多数是从已经采集好的数据样本中截取一段（长度满足3*N-1）进行数据分析。现在问题就来了，通过截取怎么求取谐波的相位（或者相对相位）？
       现在举例如下——数据中含两个频率（50Hz、100Hz），相位分别为（60度、100度），模拟产生数据的程序如下：

1. close all;clear all;clc;
   
2. fs=512;
   
3. N=512;
   
4. FSX=fs/N;
   
5. %% ***模拟数据产生***
   
6. t=-N+1:10*N-1;
   
7. A1=1.0; f1=50; ph1=60;
   
8. A2=0.3; f2=100; ph2=100;
   
9. y=A1*cos(2*pi*t*f1/fs+ph1*pi/180)+A2*cos(2*pi*t*f2/fs+ph2*pi/180);
   
10. s=y(1:3*N-1); %假设设置的情况满足校正程序的要求，即中心样点过零
    
11. % s=y(1+117:3*N-1+117); % %假设另一种截取数据情况。
    
12. 
    
13. fid=fopen('s.txt','wt');
    
14. fprintf(fid, '%f\n',s);
    
15. fclose(fid);
    
16. disp('数据产生成功！');

复制代码

利用全相位时移相位差法，求取参数程序如下：

1. close all;clear all;clc;
   
2. fs=512;
   
3. N=512;
   
4. FSX=fs/N;
   
5. %% ***数据导入***
   
6. s=load('s.txt');
   
7. s=s';
   
8. 
   
9. %% ***全相位时移相位差校正法***
   
10. win=hanning(N)';win1=hann(N)';
    
11. win2=conv(win,win1);
    
12. win2=win2/sum(win2);
    
13. w=pi*2;
    
14. s1=s(1:2*N-1); %第1组（2N-1）个数据
    
15. y1=s1.*win2;
    
16. y1a=y1(N:end)+[0 y1(1:N-1)];
    
17. out1=fft(y1a,N);
    
18. a1=abs(out1);
    
19. p1=mod(angle(out1),2*pi);
    
20. s2=s(1+N:3*N-1); %第2组（2N-1）个数据
    
21. y2=s2.*win2;
    
22. y2a=y2(N:end)+[0 y2(1:N-1)];
    
23. out2=fft(y2a,N);
    
24. a2=abs(out2);
    
25. p2=mod(angle(out2),2*pi);
    
26. 
    
27. %% ***其中频率f可以搜索峰值得到，此处直接设定为了简便！***
    
28. rr=[round(50/FSX) round(100/FSX)];
    
29. 
    
30. 
    
31. %% ***校正程序***
    
32. for k=1:length(rr)
    
33. dp=p2(rr(k)+1)-p1(rr(k)+1);
    
34. dph=mod(dp,2*pi); %由于频率较小，为减小噪声信号的干扰，作如下改进频率和幅值
    
35. if dph>pi %dph即为两序列谱峰线间的相位差，（0~2*pi）
    
36. dph=dph-2*pi;
    
37. elseif dph<-pi
    
38. dph=dph+2*pi;
    
39. end
    
40. g0=dph/pi/2;
    
41. if g0==0
    
42. h=2;
    
43. else h=2*pi*g0.*(1-g0.*g0)./sin(pi*g0);
    
44. end
    
45. 
    
46. AA=abs((h.^2).*a1)/2;
    
47. fff(k)=(rr(k)+g0)*FSX;
    
48. aaa(k)=AA(rr(k)+1);
    
49. ppp(k)=p1(rr(k)+1)*180/pi;
    
50. 
    
51. fprintf('%f\t%f\t%f\n',fff(k),aaa(k),ppp(k));
    
52. end

复制代码

第一种情况：

1. s=y(1:3*N-1); %假设设置的情况满足校正程序的要求，即中心样点过零
   

复制代码

      结果：

1. 50.000000 1.000000 60.000000
   
2. 100.000000 0.300000 99.999994

复制代码

第二种情况：

1. s=y(1+117:3*N-1+117); % %假设另一种截取数据情况。
   

复制代码

      结果：

1. 50.000000 1.000000 213.281250
   
2. 100.000000 0.300000 46.562494

复制代码

这样子，相位就不准确了。不过没关系，我们可以通过补偿校正相位，在ppp(k)=p1(rr(k)+1)*180/pi; 的下面添加如下语句：

1. ppp(k)=ppp(k)-mod(fff(k)/FSX*(117)/N*360,360);
   
2. ppp(k)=mod(ppp(k),360);

复制代码

即可得到第一种情况的结果。
【问题来了】：
如此相位校正，是假设我们已知截取的数据段距离117点（相对t=-N+1:3N-1的情况）。在工程中，我们不知道截取的一段（3N-1长度）到底是哪一段，如何处理得到正确的相位（或者正确的相位差[如此处：40度]）？

希望老师解答一下，不胜感激！

zhwang554

回复 10 # wangze2008scu 的帖子

碥定一个正弦波,要3个参数:振幅,频率和相位,
这相位是任一己知时刻的相位,可以是初相是t=0的相位,或是距离117点的相位.
Apfft测的是2N-1数据点中间点的相位,即时间为取样点N的时刻的相位,这个相位的时刻知道,相位值也知道,这个正弦波就全确定了,它的波形全知道了,即任一点的相位都可求出,这有什么问题

你不是从117点算出初相(笫一个取样的相位).你也可以从其它已知时刻的相位(如apfft的笫N个取样)算出笫一个取样的相位,并算出任一其它对刻取样点的相位

确定一个正弦波不一定要知初相,可以是任一己知时刻的相位,

wangze2008scu

回复 11 # zhwang554 的帖子

王老师，您好！首先感谢您的回复。
前面我可能描述得较细，我是想问如何从任意一段采样数据中求出这两个频率分量的相位ph1=60度、ph2=100度，或者求两者的相位差为40度？？

wangze2008scu

回复 11 # zhwang554 的帖子

举个实例：
        我这儿有一堆现在已经采样好的数据s.txt（见附件），已知采样频率fs=512Hz，余弦多次谐波信号。
怎么分离出个两个频率分量的相位ph1=60、ph2=100,或者两者的相位差为40度？
该如何处理呢？这个问题困扰我好久了，请王老解答一下，谢谢！

zhwang554

回复 13 # wangze2008scu 的帖子

s.txt中笫一个采样相位分别是60度和100度吗? 你如何做到的? 我算一下好样不是. 那一个采样相位是60度和100度.

相位差一般片指两个相同频率的任一时刻的相位差值,它们在任一时刻的相位差是相同的.

你说的两个不同频率(50Hz,100Hz)的40度初相相位差只在第一个取样存在,其它点相位差值是不同的

wangze2008scu

回复 14 # zhwang554 的帖子

首先非常感谢王老这么晚了还回复！
    是这样子的：我就是按照10楼中“模拟数据产生”程序，其中s=y(1+151:3*N-1+151);产生数据文件s.txt，各频率参数和程序中给出的不变。
    实际工程中就是这样子的：只能已知采样频率fs(如：10240Hz)，然后就是一串数据（类似s.txt中的数据），然后我用全相位时移相位差的方法进行求解，频率和幅值还好说，相位真心不知道应该怎么计算？  一般算谐波相位（此处的2次谐波100Hz）都是要求相对于基波(50Hz)的相位，无论算出来基波相位是多少，此处的2次谐波相位=基波相位+40。
    我就比较不太清楚，工程中的一段数据拿来分析时，往往就只知道采样频率，相位计算时到底以哪个为参照？