讨论：在matlab中如何实现由系统的阶跃响应求系统的传递函数

jimin

在matlab中如何实现由系统的阶跃响应求系统的传递函数
假设该系统为二阶系统

[ 本帖最后由 jimin 于 2006-11-20 19:13 编辑 ]

lxq

由阶跃响应求系统的传递函数?

jimin

原帖由 lxq 于 2006-11-9 12:19 发表
由阶跃响应求系统的传递函数?

是的，已知一个系统的阶跃响应曲线，要确定一个二阶系统的模型，使之输入阶跃信号后的阶跃响应与已知的曲线类似，应该要用到系统辨识的知识，我系统辨识了解的不多，那位朋友做过的类似的能不能给一个例子或者给我一些思路，先谢谢了

[ 本帖最后由 jimin 于 2006-11-9 12:39 编辑 ]

songzy41

我这里给出一个求阶跃响应的例子，不知是否是楼主所需要的：
a=[1 1 1/4];    %已知系统的参数a和b
b=[1];
t=0:15;           %设置时间
x=ones(1,length(t));   %设置阶跃函数
y=filter(b,a,x);           %计算阶跃响应
plot(t,y);
title('离散系统阶跃响应')
xlabel('t');
ylabel('Am(k)')

jimin

谢谢songzy41
我的刚好是逆过程
即知道系统的阶跃响应，要求系统模型的参数，可假设系统为二阶

songzy41

一般在已知系统的激冲响应下，能方便地求出系统的参数（b,a-也就能列出系统的传递函数）。但在已知系统的阶跃响应的条件下怎么来求？我是这样考虑的：激冲函数的积分（对时间）是阶跃函数，即阶跃函数的微分是激冲函数。故把阶跃响应微分，得激冲响应，从激冲响应求得系统参数。在数字系统中用差分代替微分。如果假设阶跃响应的时间函数是y，则求系统参数为：
y=[0 y];        
for k=1 : N
    y1(k)=y(k+1)-y(k);
end
[ma1,ar1] = prony (y1, p, p);
ma1为传递函数（Z变换的）中分子的系数，ar1为传递函数中分母的系数。

songzy41

为了证明这方法的正确性，以下有一个证明的程序，对于任意一个二阶系统假设了它的参数a和b。
a=[1 1 0.25];     %假设系统参数a和b
b=[1];
t=0:15;
N=length(t);
fs=200;           %假设采样频率，仅为计算系统激冲响应用
[x1,t1]=impz(b,a,t,fs);   %计算系统激冲响应
p=2;
[ma,ar] = prony (x1, p, p); %从系统激冲响应计算系统参数

x=ones(1,length(t));  %设置阶跃函数
y=filter(b,a,x);      %计算系统阶跃响应
plot(t1,x1,'r','linewidth',2); hold on;
plot(t1,y,'g');
title('离散系统阶跃响应')
xlabel('t');
ylabel('Am(k)');
y=[0 y];        %从系统阶跃响应求系统激冲响应
for k=1 : N
    y1(k)=y(k+1)-y(k);
end
%把系统激冲响应、阶跃响应和计算得的激冲响应曲线图
plot(t1,y1,'b'); hold off; grid;   
legend('激冲响应','阶跃响应','计算得激冲响应');
[ma1,ar1] = prony (y1, p, p);
fprintf('激冲响应\n');
fprintf('ma=%5.6f  %5.6f  %5.6f\n',ma);
fprintf('ma=%5.6f  %5.6f  %5.6f\n',ar);
fprintf('计算得激冲响应\n');
fprintf('ma=%5.6f  %5.6f  %5.6f\n',ma1);
fprintf('ma=%5.6f  %5.6f  %5.6f\n',ar1);
从试验来看，这方法不限于二阶系统，可适用于高阶系统。

songzy41

再从数学上来证明阶跃响应的微分是激冲响应：

jimin

再次感谢  
songzy41
你的方法我先试试

pigmin

谢谢，学习中

jimin

我用songzy41提供的方法做了下，可是求出来的阶跃曲线和原始的阶跃响应曲线相差很远，不知道问题出在哪里:@(
songzy41如果有空的话，麻烦您帮我看下，谢谢！！:@)
ttt=importdata('ttt.txt');
ttt=ttt';
p=ttt(1,:);
t=ttt(2,:);
plot(p,t)
q=2;
t1=[0 t];
for k=1:2500
    y(k)=t1(k+1)-t1(k);
end
[ma1,ar1] = prony (y, q, q);
x=ones(1,2500);
y1=filter(ma1,ar1,x);  
plot(p,t,p,y1)

songzy41

原帖由 jimin 于 2006-11-17 18:12 发表
我用songzy41提供的方法做了下，可是求出来的阶跃曲线和原始的阶跃响应曲线相差很远，不知道问题出在哪里:@(
songzy41如果有空的话，麻烦您帮我看下，谢谢！！:@)
ttt=importdata('ttt.txt');
ttt=ttt';
p=tt ...

jimin，仔细研究了你提供的数据，作了如下的处理：
1，因为原始的数据起伏很大，用滑动窗求平均的方法，求出原始阶跃响应的光滑数据；
2，把数据只取变化较大的一部分，即201-1200的一段，并下采样(10:1)；
3，求差分，得到激冲响应曲线（核心的方法还是按以前贴子所说的）；
4，对激冲响应用ARMA方法分解，并求阶跃响应。
程序清单如下：
clear all; clc;
ttt=importdata('ttt.txt');
ttt=ttt';
p=ttt(1,:);
t=ttt(2,:);
N=length(p);
%plot(p,t)
%用滑动窗进行平滑处理
M=10;
n1=1:N-M;
xx=zeros(1,N-M);
for i=1 : N-M
    for k=1 : M
        xx(i)=xx(i)+t(i-1+k);
    end
end
xx=xx/M;
%只取201~1200之间的数据;
x=xx(201:1200);
plot(n1,t(n1),'b',n1,xx(n1),'r'); grid;
legend('原始数据','平滑数据');
%下采样
y1=x(1:10:end);
%求差分和ARMA分析
y2=[0 y1];
for k=1 : 100
    z(k)=y2(k+1)-y2(k);
end
figure(2)
n2=0:99;
q=40;
[ma1,ar1] = prony (z, q, q);
%求估算的阶跃响应
x2=ones(1,100);
y3=filter(ma1,ar1,x2);
plot(n2,y1,'b','linewidth',2); hold on;
plot(n2,y3,'r'); grid;
legend('原始数据','估算数据');
其中阶数用到40阶才能较好地重合，远远不像所说的2-3阶。得到的估算阶跃响应和原始的阶跃响应比较如下，如果谁有更好的方法不妨也能贴出来，供大家学习。

[ 本帖最后由 songzy41 于 2006-11-20 18:50 编辑 ]

jimin

已知一个系统的阶跃响应曲线，若系统的阶数为三阶，求系统的模型，若阶数未知，求一个系统的模型。
程序在matlab6.5下可以运行，在7.0的版本可能会有些问题
附件中有数据文件，程序及word文档
ttt=importdata('ttt.txt');
ttt=ttt';
global zlast;
p=ttt(1,:);
t=ttt(2,:);
z=t;
u=ones(1,2500);
h=[-z(3:2499)',-z(2:2498)',-z(1:2497)', u(4:2500)',u(3:2499)',u(2:2498)',u(1:2497)'];
beta=[0.001 0.001 0.001 0.001 0.001 0.001 0.001]';
zold=z(4:2500);
betafit = nlinfit(h,zold,@myjieyue,beta)
a1 = betafit(1);
a2 = betafit(2);
a3 = betafit(3);
b1 = betafit(4);
b2 = betafit(5);
b3=betafit(6);
b4=betafit(7);
zlast=zlast';
zlast=[0.128 0.136 0.128 zlast];
plot(p,t,p,zlast)
legend('实际系统阶跃响应曲线','三阶系统阶跃响应曲线');
xlabel('2500相当于时间1秒钟');
ylabel('电压（V)');
disp('误差为')
err=z-zlast;
errsum=sum(abs(err))
figure
plot(p,err)
xlabel('2500相当于时间1秒钟');
ylabel('误差');
num=[b1 b2 b3 b4];
den=[1 a1 a2 a3];
disp('三阶系统的传递函数为')
h=tf(num,den)
function znew =myjieyue(beta,x)
global zlast;
a1 = beta(1);
a2 = beta(2);
a3 = beta(3);
b1 = beta(4);
b2 = beta(5);
b3=beta(6);
b4=beta(7);
x1 = x(:,1);
x2 = x(:,2);
x3 = x(:,3);
x4 = x(:,4);
x5 = x(:,5);
x6 = x(:,6);
x7 = x(:,7);
znew=-a1*x1-a2*x2-a3*x3+b1*x4+b2*x5+b3*x6+b4*x7;
zlast=znew;

误差为
errsum =
5.4474
三阶系统的传递函数为
Transfer function:
-0.005596 s^3 - 0.004943 s^2 + 0.01193 s + 3.915e-010
------------------------------------------------------------------------

         s^3 + 0.2029 s^2 + 0.923 s - 0.1272

vib

正在学习你们的成果，我的感觉既然阶跃响应已知了，用laplace得到他的拉氏变换，再除以阶跃函数的拉氏变换不就行了？s=jw,我是不是没弄明白你们在讲什莫？

jimin

原帖由 vib 于 2006-12-7 19:59 发表
正在学习你们的成果，我的感觉既然阶跃响应已知了，用laplace得到他的拉氏变换，再除以阶跃函数的拉氏变换不就行了？s=jw,我是不是没弄明白你们在讲什莫？

阶跃响应曲线是一些数据，你觉得能通过你所说的方法做呢？
数据都在上面，如果能做的，做出来，让大家学习一下