高手帮我看看程序,看有什么问题(纵轴为相位差,横轴为频率图)?

matlab-vib

例图

如何根据已知两信号怎么得到如图所示的例图?信号的频率随时间呈非线性变化,信号为2ms采样,长为8s,我的程序:

1. a1=load('1.txt');
   
2. a2=load('2.txt');
   
3. b=length(a1);
   
4. f=0:500/b:500-500/b;
   
5. p1=unwrap(angle(fft(a1)));
   
6. p2=unwrap(angle(fft(a2)));
   
7. plot(f,p1-p2)

复制代码

和示例图结果相差很远,...
高手帮我看看程序,看看问题出在哪儿?

[ 本帖最后由 matlab-vib 于 2009-1-18 17:30 编辑 ]

songzy41

从楼主是问题来看，不是要求每个频率的相位差，而是求两信号的相位差。在以前的帖子中已提到FFT不适用于非稳态信号，这里用hilbert变换求出信号的相位，再求两信号的相位差，和信号的频率。从做的图看出信号的频率是非线性的。

1. x=load('1.txt');
   
2. y=load('2.txt');
   
3. subplot 211; plot(x);title('信号1');
   
4. subplot 212; plot(y);title('信号2');
   
5. fs=500;
   
6. xx=hilbert(x);
   
7. yy=hilbert(y);
   
8. AX=unwrap(angle(xx));
   
9. AY=unwrap(angle(yy));
   
10. N=length(xx);
    
11. m=1:N;
    
12. DA=AX-AY;
    
13. figure
    
14. subplot 211; plot(m,AX,'r',m,AY,'b'); grid;
    
15. title('两信号的相位角');
    
16. subplot 212; plot(m,DA); grid;
    
17. title('两信号的相位差');
    
18. ffx=diff(AX);
    
19. ff=ffx*fs;
    
20. N1=length(ff);
    
21. m1=1:N1;
    
22. figure
    
23. plot(m1,ff,'r'); grid;
    
24. title('测算的频率');

复制代码

matlab-vib

SONG老师,你好,
首先很感谢你的帮助;
还有几个问题,
1,我的例图是这两个信号的实际的相位差与频率图(国外的软件做的),不是示意图,比较了一下你的结果,我把你的图放到频率为10-64Hz的范围,还是和例图有十几度的相位误差,不知道是什么原因?
2.你做的那个频率测算用的:ffx=diff(AX);是相位角的导数,这里和频率的关系是?不太明白还请指点

songzy41

对于笫1个问题我也无法回答，我只是根据对相位角和频率的计算得到的结果。笫2个问题是计算瞬时频率的基本关系，求出的ffx是圆频率，实际频率还要除2pi。可以参看一下程乾生 “数字信号处理” 笫110页 “实连续信号的包络、瞬时相位、瞬时频率”一节。很多文献也给出过由hilbert变换求瞬时频率的关系。

matlab-vib

感谢SONG老师:victory:

hitman1990

为什么看不了图