已知输入输出 求传递函数 四种方法比较

sunyuxinhe

看了论坛上很多关于测量传递函数的帖子：
具有代表性的：
http://forum.vibunion.com/thread-110984-1-1.html

里面有H1,H2,Hv的程序
我自己找了个例子实现了一下，发现这三个方法均差与互谱除以自谱的方法。
互谱和自谱分别用的是matlab中的psd,csd来做。


有兴趣的可以运行一下：

1. %% 传递函数
   
2. clear
   
3. close all
   
4. num=1;
   
5. den=[1  6.2832e+001 9.8696e+004];
   
6. sys=tf(num,den);
   
7. %% 输入
   
8. nfft=1024;
   
9. x=randn(1,nfft);
   
10. fs=1000;
    
11. t=(1:nfft)/fs;
    
12. f1=(0:nfft/2-1)/(nfft/2)*fs/2;
    
13. %% 求输出
    
14. y=lsim(sys,x,t);
    
15. y=y';

复制代码

1. %%  GXF./GXX
   
2. X=fft(x);
   
3. Y=fft(y);
   
4. X=X(1:end/2);
   
5. Y=Y(1:end/2);
   
6. Sxx=X.*conj(X);
   
7. Syx=Y.*conj(X);
   
8. Txy=Syx./Sxx;
   
9. figure;plot(f1,abs(Txy)/max(abs(Txy)));
   
10. ang=angle(Txy);
    
11. for i=1:length(ang)
    
12.     if ang(i)<0
    
13.          ang(i)=ang(i)+2*pi;
    
14.     end
    
15. end
    
16. figure;plot(f1,ang);

复制代码

1. %%   GXF./GXX
   
2. Syy=Y.*conj(Y);
   
3. Sxy=X.*conj(Y);
   
4. Txy=Syy./Sxy;
   
5. figure;plot(f1,abs(Txy)/max(abs(Txy)));
   
6. ang=angle(Txy);
   
7. for i=1:length(ang)
   
8.     if ang(i)<0
   
9.          ang(i)=ang(i)+2*pi;
   
10.     end
    
11. end
    
12. figure;plot(f1,ang);
    

复制代码

1. %%  Hv
   
2. for ii=1:nfft/2;
   
3.         G=[Syy(1,ii),Sxy(1,ii);Syx(1,ii),Sxx(1,ii)];
   
4.         [xx,d]=eig(G);
   
5.         orig_lambda=diag(d);
   
6.         [Y,I]=sort(real(orig_lambda));
   
7.         lambda=orig_lambda(I);
   
8.         psi=xx(:,I);
   
9.         Hv(1,ii)=-psi(2,1)/psi(1,1);
   
10. end;
    
11. Txy=Hv;
    
12. figure;plot(f1,abs(Txy)/max(Txy));
    
13. ang=angle(Txy);
    
14. for i=1:length(ang)
    
15.     if ang(i)<0
    
16.          ang(i)=ang(i)+2*pi;
    
17.     end
    
18. end
    
19. figure;plot(f1,ang);

复制代码

1. %%  PSD 互谱/自谱
   
2. 
   
3. window=gausswin(nfft/4)';
   
4. [Sxx,f1]=psd(x,nfft/4,fs,window);
   
5. [Sxy,f1]=csd(x,y,nfft/4,fs,window);
   
6. Txy=Sxy./Sxx;
   
7. figure;plot(f1,abs(Txy)/max(abs(Txy)));
   
8. ang=angle(Txy);
   
9. for i=1:length(ang)
   
10.     if ang(i)<0
    
11.          ang(i)=ang(i)+2*pi;
    
12.     end
    
13. end
    
14. figure;plot(f1,ang);

复制代码

sunyuxinhe

自己是第一次涉足，不对的地方还请提宝贵意见

westrongmc

sunyuxinhe 发表于 2015-1-19 21:12
自己是第一次涉足，不对的地方还请提宝贵意见

所谓的利用matlab的 PSD 互谱/自谱中，进行了四次平均，加了窗函数；
而前面的H1，H2，Hv算子中的互谱、自谱没有进行平均，没有加窗。

呵呵，不在一个条件下，可以对比吗？

sunyuxinhe

westrongmc 发表于 2015-1-20 13:54
所谓的利用matlab的 PSD 互谱/自谱中，进行了四次平均，加了窗函数；
而前面的H1，H2，Hv算子的中的互谱 ...

多谢所提建议！我也是刚开始接触，看来平均和窗函数确实很重要！

star_cheney

正好在做，只是刚做到滤波截断。过来学习了