系统辨识作业给两组数据求频谱分析.docx
《系统辨识作业给两组数据求频谱分析.docx》由会员分享，可在线阅读，更多相关《系统辨识作业给两组数据求频谱分析.docx（20页珍藏版）》请在冰豆网上搜索。
系统辨识作业给两组数据求频谱分析
模型一数据一：
clc
fs=500;
t=[1/fs:
1/fs:
1];
a=importdata（'C:
\Documents and Settings\Administrator\桌面\系统辨识数据（两组）\model1\model1_1_1a.txt'）;
N=length（a）;
p1=a（:
1）; %得到采样输入数据p，下面对p进行频谱分析
figure
（1）
plot（t,p1）;
grid on
title（'输入数据p1'）; xlabel（'t'）;ylabel（'p1'）
p2=a（:
2）; %得到采样输入数据p，下面对p进行频谱分析
figure
（2）
plot（t,p2）;
grid on
title（'输入数据p2'）; xlabel（'t'）;ylabel（'p2'）
Y1=fft（p1）;
Y2=fft（p2）;
Y=Y2./Y1;
magY=abs（Y（1:
1:
N/2））*2/N;
f=（0:
N/2-1）'*fs/N;
figure（3）
plot（f,magY）;
h=stem（f,magY,'fill','--'）;
set（h,'MarkerEdgeColor','red','Marker','*'）;
grid on
title（'频谱图'）; xlabel（'f （Hz）'）; ylabel（'幅值'）
去掉
h=stem（f,magY,'fill','--'）;set（h,'MarkerEdgeColor','red','Marker','*'）;
模型一数据二：
clc
fs=500;
t=[1/fs:
1/fs:
1];
a=importdata（'C:
\Documents and Settings\Administrator\桌面\系统辨识数据（两组）\model1\model1_1_2a.txt'）;
N=length（a）;
p1=a（:
1）; %得到采样输入数据p，下面对p进行频谱分析
figure
（1）
plot（t,p1）;
grid on
title（'输入数据p1'）; xlabel（'t'）;ylabel（'p1'）
p2=a（:
2）; %得到采样输入数据p，下面对p进行频谱分析
figure
（2）
plot（t,p2）;
grid on
title（'输入数据p2'）; xlabel（'t'）;ylabel（'p2'）
Y1=fft（p1）;
Y2=fft（p2）;
Y=Y2./Y1;
magY=abs（Y（1:
1:
N/2））*2/N;
f=（0:
N/2-1）'*fs/N;
figure（3）
plot（f,magY）;
grid on
title（'频谱图'）; xlabel（'f （Hz）'）; ylabel（'幅值'）
模型一数据三：
clc
fs=500;
t=[1/fs:
1/fs:
1];
a=importdata（'C:
\Documents and Settings\Administrator\桌面\系统辨识数据（两组）\model1\model1_1_3a.txt'）;
N=length（a）;
p1=a（:
1）; %得到采样输入数据p，下面对p进行频谱分析
figure
（1）
plot（t,p1）;
grid on
title（'输入数据p1'）; xlabel（'t'）;ylabel（'p1'）
p2=a（:
2）; %得到采样输入数据p，下面对p进行频谱分析
figure
（2）
plot（t,p2）;
grid on
title（'输入数据p2'）; xlabel（'t'）;ylabel（'p2'）
Y1=fft（p1）;
Y2=fft（p2）;
Y=Y2./Y1;
magY=abs（Y（1:
1:
N/2））*2/N;
f=（0:
N/2-1）'*fs/N;
figure（3）
plot（f,magY）;
grid on
title（'频谱图'）; xlabel（'f （Hz）'）; ylabel（'幅值'）
模型二数据一：
clc
fs=500;
t=[1/fs:
1/fs:
1];
a=importdata（'C:
\Documents and Settings\Administrator\桌面\系统辨识数据（两组）\model2\model2_1_1a.txt'）;
N=length（a）;
p1=a（:
1）; %得到采样输入数据p，下面对p进行频谱分析
figure
（1）
plot（t,p1）;
grid on
title（'输入数据p1'）; xlabel（'t'）;ylabel（'p1'）
p2=a（:
2）; %得到采样输入数据p，下面对p进行频谱分析
figure
（2）
plot（t,p2）;
grid on
title（'输入数据p2'）; xlabel（'t'）;ylabel（'p2'）
Y1=fft（p1）;
Y2=fft（p2）;
Y=Y2./Y1;
magY=abs（Y（1:
1:
N/2））*2/N;
f=（0:
N/2-1）'*fs/N;
figure（3）
plot（f,magY）;
grid on
title（'频谱图'）; xlabel（'f （Hz）'）; ylabel（'幅值'）
模型二数据二：
clc
fs=500;
t=[1/fs:
1/fs:
1];
a=importdata（'C:
\Documents and Settings\Administrator\桌面\系统辨识数据（两组）\model2\model2_1_2a.txt'）;
N=length（a）;
p1=a（:
1）; %得到采样输入数据p，下面对p进行频谱分析
figure
（1）
plot（t,p1）;
grid on
title（'输入数据p1'）; xlabel（'t'）;ylabel（'p1'）
p2=a（:
2）; %得到采样输入数据p，下面对p进行频谱分析
figure
（2）
plot（t,p2）;
grid on
title（'输入数据p2'）; xlabel（'t'）;ylabel（'p2'）
Y1=fft（p1）;
Y2=fft（p2）;
Y=Y2./Y1;
magY=abs（Y（1:
1:
N/2））*2/N;
f=（0:
N/2-1）'*fs/N;
figure（3）
plot（f,magY）;
grid on
title（'频谱图'）; xlabel（'f （Hz）'）; ylabel（'幅值'）
模型二数据三：
clc
fs=500;
t=[1/fs:
1/fs:
1];
a=importdata（'C:
\Documents and Settings\Administrator\桌面\系统辨识数据（两组）\model2\model2_1_3a.txt'）;
N=length（a）;
p1=a（:
1）; %得到采样输入数据p，下面对p进行频谱分析
figure
（1）
plot（t,p1）;
grid on
title（'输入数据p1'）; xlabel（'t'）;ylabel（'p1'）
p2=a（:
2）; %得到采样输入数据p，下面对p进行频谱分析
figure
（2）
plot（t,p2）;
grid on
title（'输入数据p2'）; xlabel（'t'）;ylabel（'p2'）
Y1=fft（p1）;
Y2=fft（p2）;
Y=Y2./Y1;
magY=abs（Y（1:
1:
N/2））*2/N;
f=（0:
N/2-1）'*fs/N;
figure（3）
plot（f,magY）;
grid on
title（'频谱图'）; xlabel（'f （Hz）'）; ylabel（'幅值'）
判别系统的阶次;
%模型仿真与参数估计
%产生输入输出数据
clearall;
fs=500;
t=[1/fs:
1/fs:
1];
a=importdata（'C:
\DocumentsandSettings\Administrator\桌面\系统辨识数据（两组）\model1\model1_1_1a.txt'）;
N=length（a）;%数据长度
p1=a（:
1）;%得到采样输入数据p，下面对p进行频谱分析
figure
（1）
plot（t,p1）;
gridon
title（'输入数据p1'）;xlabel（'t'）;ylabel（'p1'）
p2=a（:
2）;%得到采样输出数据p，下面对p进行频谱分析
figure
（2）
plot（t,p2）;
gridon
title（'输入数据p2'）;xlabel（'t'）;ylabel（'p2'）
A=[1-1.50.7];B=[010.5];C=[1-10.2];
th0=idpoly（A,B,C）;%创建一个armax模型
y=sim（th0,[p1p2]）;%在线模型仿真
z=iddata（y,p1）;%取得输入输出值
figure（3）
plot（z）%绘制输入输出曲线
%armax模型阶次的估计
NN=struc（1:
2,1:
2,1:
4）;
Loss_fun=arxstruc（z,z,NN）;
order=selstruc（Loss_fun,'aic'）;
order=[order
（1）,order
（2）,1,order（3）];
order;%模型阶次
%模型参数的估计
Model_para=armax（z,order）;%估计armax模型参数
present（Model_para）;%显示辨识结果
figure（4）
compare（z,Model_para）%实际参数与估计参数比较
程序结果
分析：
z=iddata（y,p1）
Timedomaindatasetwith500samples.
Samplinginterval:
1
OutputsUnit（ifspecified）
y1
InputsUnit（ifspecified）
u1
NN=struc（1:
2,1:
2,1:
4）
NN=
111
112
113
114
121
122
123
124
211
212
213
214
221
222
223
224
Loss_fun=arxstruc（z,z,NN）
Loss_fun=
1.0e+003*
Columns1through11
0.00020.00030.00030.00030.00020.00030.00030.00030.00010.00010.0001
0.00100.00100.00100.00100.00100.00100.00100.00100.00200.00200.0020
0.00100.00100.00100.00100.00200.00200.00200.00200.00100.00100.0010
0.00100.00200.00300.00400.00100.00200.00300.00400.00100.00200.0030
Columns12through17
0.00010.00000.00010.00010.00010.4950
0.00200.00200.00200.00200.00201.4925
0.00100.00200.00200.00200.00200
0.00400.00100.00200.00300.00400
order=selstruc（Loss_fun,'aic'）
order=
221
order=[order
（1）,order
（2）,1,order（3）]
order=
2211
order
order=
2211
Model_para=armax（z,order）
Discrete-timeIDPOLYmodel:
A（q）y（t）=B（q）u（t）+C（q）e（t）
A（q）=1-1.507q^-1+0.7033q^-2
B（q）=1.14q^-1+0.548q^-2
C（q）=1-0.9697q^-1
EstimatedusingARMAXfromdatasetz
Lossfunction0.0215163andFPE0.0219484
Samplinginterval:
1
present（Model_para）
Discrete-timeIDPOLYmodel:
A（q）y（t）=B（q）u（t）+C（q）e（t）
A（q）=1-1.507（+-0.003833）q^-1+0.7033（+-0.003341）q^-2
B（q）=1.14（+-0.02723）q^-1+0.548（+-0.02639）q^-2
C（q）=1-0.9697（+-0.01052）q^-1
EstimatedusingARMAXfromdatasetz
Lossfunction0.0215163andFPE0.0219484
Samplinginterval:
1
Created:
18-Jun-201419:
51:
42
Lastmodified:
18-Jun-201419:
53:
39
升级会员 