一封感谢信（欢迎大家讨论）

chengtianmmz

前不久，我所上的机械系统建模与动态分析课上，老师布置了一个作业大致如下：用matlab作系统辨识，要求分别使用最小二乘、递推最小二乘、广义最小二乘三种方法实现，还要求验证数据饱和，显示预报误差等等。当时把我愁坏了，因为我一点儿编程的基础都没有，对matlab更是一无所知，难度实在是太大了。于是开始上网疯狂地找相关资料，碰巧来到了这个论坛，下载了一些有用的东西。但是无意间看到了eight的一封信，希望大家自己动手动脑，不要动不动就到处求助源代码。我突然感到很不好意思，好像是自己偷偷摸摸做坏事被别人发现了一样，于是就下决心自己完成。
当然刚开始困难很大，不得不借鉴了别人的一些东西，比如信号源的产生就是在别人M序列基础上改成的逆M序列。后来我就觉得写起来顺手多了，虽然方法很笨很烦琐，但毕竟是独立完成的，心里感到很高兴。现在和大家分享第一阶段的成果，希望大家能多提意见，使这个程序更加完美（因为觉得自己写得太老土，太不专业了）。
首先运行input_output程序，得到输入信号数据（存放在input_m.txt中）和输出数据（存放在output1.txt中）,并得到信号波形图：
clear;close all;
a=input('a=')
b=input('b=')
c=20*log10(b);
v=wgn(1,a,c);
B1=0.2;
A=[1,-0.5,0.7];
x1=v;
y1=filter(B1,A,x1);
X1=1;X2=0;X3=1;X4=0; for i=1:a    %1#
   Y4=X4;  Y3=X3;  Y2=X2;   Y1=X1;S1=0;S2=1;S3=0;S4=1;
   X4=Y3;  X3=Y2;   X2=Y1;  
   X1=xor(Y3,Y4);
   Z1=xor(X1,S1);
   Z2=xor(X2,S2);
   Z3=xor(X3,S3);
   Z4=xor(X4,S4);
   if Z4==0   
       U(i)=0;
   else
   U(i)=Z4;
end
end
M=U;
B2=[0,0.2,-0.1];
A=[1,-0.5,0.7];
x2=U;
y2=filter(B2,A,x2);
y=y1+y2;
%绘图
k=1:a;
subplot(3,1,1);plot(k,y1)
title('滤波后的白噪声')
xlabel('n')
ylabel('幅值')
subplot(3,1,2);plot(k,y2)
title('滤波后的逆M序列')
xlabel('n')
ylabel('幅值')
subplot(3,1,3);plot(k,y)
title('Output1')
xlabel('n')
ylabel('幅值')
fid=fopen('input_m.txt','wt');
fprintf(fid,'%f\n',M);
fclose(fid);
fid1=fopen('output1.txt','wt');
fprintf(fid1,'%f\n',y);
fclose(fid1);
然后运行LS1程序，得到估计参数的值F=[-0.5001   0.7003   0.1998   -0.1001]'。
fid=fopen('input_m.txt','r+');
C13=fscanf(fid,'%f',a-1);
fclose(fid);
fid=fopen('input_m.txt','r+');
C14=fscanf(fid,'%f',a-2);
fclose(fid);
fid1=fopen('output1.txt','wt');
fprintf(fid1,'%12.8f\n',y);
fclose(fid1);
fid1=fopen('output1.txt','r+');
C=fscanf(fid1,'%f',a);
fclose(fid1);
fid1=fopen('output1.txt','r+');
C11=fscanf(fid1,'%f',a-1);
fclose(fid1);
fid1=fopen('output1.txt','r+');
C12=fscanf(fid1,'%f',a-2);
fclose(fid1);
D11=-1*[zeros(1,1);C11];
D12=-1*[zeros(2,1);C12];
D13=[zeros(1,1);C13];
D14=[zeros(2,1);C14];
D=[D11 D12 D13 D14];
E=[Rot90(D11);Rot90(D12);Rot90(D13);Rot90(D14)];
F=inv(E*D)*E*C

花如月

貌似感谢8兄，实则分享自己的成果和经验。支持一下，希望楼主能后边的成果也都拿来分享。给后来者一个不错的参考:handshake

[ 本帖最后由 eight 于 2007-6-12 22:49 编辑 ]

chengtianmmz

两者都有呵呵:lol
主要也想让大家提下意见，而且其实运行的时候还有一个警告，就是说Rot90的使用。
因为我没找到怎样求一个矩阵的转置，trsp命令用不了:@( 所以就用了一个土办法。
我会继续努力的，谢谢:@)

rocwoods

A'是A的共轭转置。A.'是A的普通转置。
当然A如果是实数矩阵，两者一样的

chengtianmmz

楼上是说直接写A.'就可以了吗？好像我试过不行。不知道是不是我机器上matlab安装的问题，因为simulink也用不了，本来想省事不编程了，可是唉:@( 卸了重装也不行，郁闷。

rocwoods

以下是我机子上运行结果。你试试，不行就是你机子的问题了
A=rand(5)
A.'
A =
    0.8147    0.0975    0.1576    0.1419    0.6557
    0.9058    0.2785    0.9706    0.4218    0.0357
    0.1270    0.5469    0.9572    0.9157    0.8491
    0.9134    0.9575    0.4854    0.7922    0.9340
    0.6324    0.9649    0.8003    0.9595    0.6787
ans =
    0.8147    0.9058    0.1270    0.9134    0.6324
    0.0975    0.2785    0.5469    0.9575    0.9649
    0.1576    0.9706    0.9572    0.4854    0.8003
    0.1419    0.4218    0.9157    0.7922    0.9595
    0.6557    0.0357    0.8491    0.9340    0.6787

VibrationMaster

A'是共扼转置,但A为复数时就不是简单的转置

chengtianmmz

试了:@) 好了谢谢。
现在正在搞递推算法......