传递函数，相位响应计算的问题求助

liushuangq05

各位老师：
       我有个明确的传递函数（图4），用它们做幅频响应图和书上的一致！但是相位响应却和书上的差别很大。所以不知道是不是我没搞明白相位响应的计算问题。
      我是这样计算相位响应的，就是将s替换为i*f(不是i*w,这个我对比了书上的结果，其中w=2*pi*f)，然后利用H(s)的虚部和实部的比值的反正切得到，就是利用matlab的angle函数。遗憾的是，这样的结果却不知道是否对（那书上的图可能也不对）。
      所以想让老师们帮指点迷津，不甚感谢。程序如下：
clc;
clear;
clf;
n1=1;
D=[0.02 0.1 0.2 0.707 1.5];
e1=D*n1;
s1=log(0.01);
s2=log(100);
s3=s1:(s2-s1)/100:s2;
s4=exp(s3);
s5=i*s4;
pi2angle=180/pi;

for k=1:1:length(D)
    for m=1:1:length(s5)
         H(m)=-s5(m)^0/(s5(m)^2+2*e1(k)*s5(m)+n1^2);
   
    end
    HH(:,k)=abs(H);         %取模
   
    Hangle(:,k)=pi2angle*angle(H);
  end
HH=20*log10(HH);
subplot(211);
semilogx(s4,(HH(:,1)),'r','LineWidth',2);
hold on
semilogx(s4,(HH(:,2)),'b','LineWidth',2);
semilogx(s4,(HH(:,3)),'c','LineWidth',2);
semilogx(s4,(HH(:,4)),'m','LineWidth',2);
semilogx(s4,(HH(:,5)),'k','LineWidth',2);
grid on
%axis tight
xlabel('角频率/(rad/s)');
ylabel('幅度/dB');
legend('D=0.02','D=0.1','D=0.2','D=0.707','D=1.5','Location','NorthEastOutside')
axis([0.01 100 -45 40]);

subplot(212);
semilogx(s4,Hangle(:,1),'r','LineWidth',2);
hold on
semilogx(s4,Hangle(:,2),'b','LineWidth',2);
semilogx(s4,Hangle(:,3),'c','LineWidth',2);
semilogx(s4,Hangle(:,4),'m','LineWidth',2);
semilogx(s4,Hangle(:,5),'k','LineWidth',2);
grid on
axis tight
xlabel('角频率/(rad/s)');
ylabel('相位/degree');
legend('D=0.02','D=0.1','D=0.2','D=0.707','D=1.5','Location','NorthEastOutside')

figure(2)
semilogx(s4,real(H))
hold on
semilogx(s4,imag(H),':')
hold off

liushuangq05

上面三个传递函数，只是分子的s的幂次不同。其中1，3是的s是偶数次，2式是奇数次幂！
书上的三个相位响应图的形式跟图1，图3相似，都是由左到右的斜坡。其中图1的纵坐标0~180度；图2为-90~90度，图3为-180~0度。

由于还没有找到由传递函数具体计算相位的明确公式和步骤，所以非常希望大家指点！

liushuangq05

1. function sysfunction2image
   
2. %Accomplished by Sqliu,2011-6-6
   
3. clear;
   
4. clc;
   
5. clf;
   
6. close all;
   
7. 
   
8. D=[0.02 0.1 0.2 0.707 1.5];
   
9. n1=1;
   
10. b=[0 1 0]; %分子
    
11. %位移，速度，加速度==> b=[1 0 0]; b=[0 1 0]; b=[0 0 1];
    
12. 
    
13. for i=1:1:length(D)
    
14. D1=D(i);
    
15. %D1=0.707;
    
16. %D1=e1/n1;
    
17. e1=2*n1*D1;
    
18. a=[1 e1 n1]; %分母
    
19. sys=tf(b,a);
    
20. pzmap(sys); %绘制零、极点
    
21. [mag,phase,w]=bode(b,a,{0.01 100}); %计算幅频特性和相频特性
    
22. 
    
23. N=length(mag);
    
24. mags(1:N,i)=20*log10(mag(1:N));
    
25. phases(1:N,i)=phase(1:N);
    
26. W(1:N,i)=w(1:N);
    
27. end
    
28. 
    
29. figure('color','w','position',[60 60 700 600]);
    
30. subplot(211);
    
31. semilogx(W(:,1),(mags(:,1)),'r','LineWidth',2);
    
32. hold on
    
33. semilogx(W(:,2),(mags(:,2)),'b','LineWidth',2);
    
34. semilogx(W(:,3),(mags(:,3)),'c','LineWidth',2);
    
35. semilogx(W(:,4),(mags(:,4)),'m','LineWidth',2);
    
36. semilogx(W(:,5),(mags(:,5)),'k','LineWidth',2);
    
37. grid on
    
38. xlabel('角频率/(rad/s)');
    
39. ylabel('幅度/dB');
    
40. legend('D=0.02','D=0.1','D=0.2','D=0.707','D=1.5','Location','NorthEastOutside')
    
41. axis tight;
    
42. 
    
43. subplot(212);
    
44. semilogx(W(:,1),(phases(:,1)),'r','LineWidth',2);
    
45. hold on
    
46. semilogx(W(:,2),(phases(:,2)),'b','LineWidth',2);
    
47. semilogx(W(:,3),(phases(:,3)),'c','LineWidth',2);
    
48. semilogx(W(:,4),(phases(:,4)),'m','LineWidth',2);
    
49. semilogx(W(:,5),(phases(:,5)),'k','LineWidth',2);
    
50. ylabel('相位(deg)');
    
51. xlabel('角频率(rad/sec)');
    
52. grid on
    
53. axis tight;
    
54. legend('D=0.02','D=0.1','D=0.2','D=0.707','D=1.5','Location','NorthEastOutside')
    
55. 
    
56. 
    
57. %figure;
    
58. %bode(b,a);
    
59. %系统特性：|H(w)|=H0*|jw-z1|*|jw-z2|*...*|jw-zm|/(|jw-s1|*|jw-s2|*...*|jw-sp|)
    
60. %phi(w)=arg(jw-z1)+arg(jw-z2)+...+arg(jw-zm)-arg(jw-s1)-arg(jw-s2)-...-arg(
    
61. %jw-sp)
    
62. %其中m为零点个数，p为极点个数。
    
63. 
    
64. 
    
65. 
    
66. 
    
67. 
    

复制代码

问题基本解决，供参考！

ljelly

要是能描述一下你解决问题的思路，和大家共享一下，会更好！

liushuangq05

问题出在：phi(w)=arg(jw-z1)+arg(jw-z2)+...+arg(jw-zm)-arg(jw-s1)-arg(jw-s2)-...-arg(%jw-sp)这个计算的方程，不能直接对H(s)反正切得到，即expand,sym2poly无法用angle函数直接取代！可参看下面的例子，这是BBVS-60s的地震仪器函数响应。

1. 
   
2. 
   

复制代码

%function sysfunc_example
%Accomplished by Sqliu,2011-6-6
clear;
clc;
clf;
close all;
syms s;
k=1.488009*10^19;   %归一化因子
r=2*pi;
%----零点----------
z1=0/r;
z2=0/r;
z3=(149.067+219.11*i)/r;
z4=(149.067-219.11*i)/r;
z5=(33.4425+264.982*i)/r;
z6=(33.4425-264.982*i)/r;
%----极点----------
p1=(-0.0740369+0.0740592*i)/r;
p2=(-0.0740369-0.0740592*i)/r;
p3=(-7.11489+266.09*i)/r;
p4=(-7.11489-266.09*i)/r;
p5=(-12.3067+248.881*i)/r;
p6=(-12.3067-248.881*i)/r;
p7=(-23.0299-223.483*i)/r;
p8=(-23.0299+223.483*i)/r;
p9=(-34.8891-183.589*i)/r;
p10=(-34.8891+183.589*i)/r;
p11=(-47.2687-130.253*i)/r;
p12=(-47.2687+130.253*i)/r;
p13=(-58.258-66.4236*i)/r;
p14=(-58.258+66.4236*i)/r;
p15=(-64.1115)/r;
%================================
Hs1=k*(s-z1)*(s-z2)*(s-z3)*(s-z4)*(s-z5)*(s-z6)*(2*pi)^6;
Hs2=(s-p1)*(s-p2)*(s-p3)*(s-p4)*(s-p5)*(s-p6)*(s-p7)*(s-p8)*...
    (s-p9)*(s-p10)*(s-p11)*(s-p12)*(s-p13)*(s-p14)*(s-p15)*(2*pi)^15;
Hs1=expand(Hs1);
Hs2=expand(Hs2);
b=sym2poly(Hs1);    %分子
a=sym2poly(Hs2);    %分母
sys=tf(b,a);
pzmap(sys);
[mag,phase,w]=bode(b,a,{0.001 100});
figure('color','w','position',[100 160 700 550]);
subplot(211);
loglog(w,mag,'color','r','Linewidth',3);grid on
axis([0.001 100 0.01 1.5]);
ylabel('Velocity Response');
title('BBVS-60: 数字地震观测系统幅频特性');
subplot(212);
phase=phase-round(max(phase)/360)*360;
semilogx(w,phase,'color','b','Linewidth',3);grid on
axis([0.001 100 -1000 200]);
ylabel('Phase[Degree]');
title('BBVS-60: 数字地震观测系统相频特性');
xlabel('Frequency[Hz]');