传递函数，相位响应计算的问题求助
本帖最后由 liushuangq05 于 2011-6-3 21:30 编辑各位老师：
我有个明确的传递函数（图4），用它们做幅频响应图和书上的一致！但是相位响应却和书上的差别很大。所以不知道是不是我没搞明白相位响应的计算问题。
我是这样计算相位响应的，就是将s替换为i*f(不是i*w,这个我对比了书上的结果，其中w=2*pi*f)，然后利用H(s)的虚部和实部的比值的反正切得到，就是利用matlab的angle函数。遗憾的是，这样的结果却不知道是否对（那书上的图可能也不对）。
所以想让老师们帮指点迷津，不甚感谢。程序如下：
clc;
clear;
clf;
n1=1;
D=;
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();
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 于 2011-6-3 21:24 编辑
上面三个传递函数，只是分子的s的幂次不同。其中1，3是的s是偶数次，2式是奇数次幂！
书上的三个相位响应图的形式跟图1，图3相似，都是由左到右的斜坡。其中图1的纵坐标0~180度；图2为-90~90度，图3为-180~0度。
由于还没有找到由传递函数具体计算相位的明确公式和步骤，所以非常希望大家指点！ function sysfunction2image
%Accomplished by Sqliu,2011-6-6
clear;
clc;
clf;
close all;
D=;
n1=1;
b=; %分子
%位移，速度，加速度==> b=; b=; b=;
for i=1:1:length(D)
D1=D(i);
%D1=0.707;
%D1=e1/n1;
e1=2*n1*D1;
a=; %分母
sys=tf(b,a);
pzmap(sys); %绘制零、极点
=bode(b,a,{0.01 100}); %计算幅频特性和相频特性
N=length(mag);
mags(1:N,i)=20*log10(mag(1:N));
phases(1:N,i)=phase(1:N);
W(1:N,i)=w(1:N);
end
figure('color','w','position',);
subplot(211);
semilogx(W(:,1),(mags(:,1)),'r','LineWidth',2);
hold on
semilogx(W(:,2),(mags(:,2)),'b','LineWidth',2);
semilogx(W(:,3),(mags(:,3)),'c','LineWidth',2);
semilogx(W(:,4),(mags(:,4)),'m','LineWidth',2);
semilogx(W(:,5),(mags(:,5)),'k','LineWidth',2);
grid on
xlabel('角频率/(rad/s)');
ylabel('幅度/dB');
legend('D=0.02','D=0.1','D=0.2','D=0.707','D=1.5','Location','NorthEastOutside')
axis tight;
subplot(212);
semilogx(W(:,1),(phases(:,1)),'r','LineWidth',2);
hold on
semilogx(W(:,2),(phases(:,2)),'b','LineWidth',2);
semilogx(W(:,3),(phases(:,3)),'c','LineWidth',2);
semilogx(W(:,4),(phases(:,4)),'m','LineWidth',2);
semilogx(W(:,5),(phases(:,5)),'k','LineWidth',2);
ylabel('相位(deg)');
xlabel('角频率(rad/sec)');
grid on
axis tight;
legend('D=0.02','D=0.1','D=0.2','D=0.707','D=1.5','Location','NorthEastOutside')
%figure;
%bode(b,a);
%系统特性：|H(w)|=H0*|jw-z1|*|jw-z2|*...*|jw-zm|/(|jw-s1|*|jw-s2|*...*|jw-sp|)
%phi(w)=arg(jw-z1)+arg(jw-z2)+...+arg(jw-zm)-arg(jw-s1)-arg(jw-s2)-...-arg(
%jw-sp)
%其中m为零点个数，p为极点个数。
问题基本解决，供参考！ 要是能描述一下你解决问题的思路，和大家共享一下，会更好！ 问题出在：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的地震仪器函数响应。
%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);
=bode(b,a,{0.001 100});
figure('color','w','position',);
subplot(211);
loglog(w,mag,'color','r','Linewidth',3);grid on
axis();
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();
ylabel('Phase');
title('BBVS-60: 数字地震观测系统相频特性');
xlabel('Frequency');
页:
[1]