声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1686|回复: 7

[综合讨论] 一封感谢信(欢迎大家讨论)

[复制链接]
发表于 2007-6-12 22:01 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
前不久,我所上的机械系统建模与动态分析课上,老师布置了一个作业大致如下:用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


output1.JPG

评分

1

查看全部评分

回复
分享到:

使用道具 举报

发表于 2007-6-12 22:07 | 显示全部楼层
貌似感谢8兄,实则分享自己的成果和经验。支持一下,希望楼主能后边的成果也都拿来分享。给后来者一个不错的参考:handshake

[ 本帖最后由 eight 于 2007-6-12 22:49 编辑 ]
 楼主| 发表于 2007-6-12 22:21 | 显示全部楼层
两者都有呵呵:lol
主要也想让大家提下意见,而且其实运行的时候还有一个警告,就是说Rot90的使用。
因为我没找到怎样求一个矩阵的转置,trsp命令用不了:@( 所以就用了一个土办法。
我会继续努力的,谢谢:@)
发表于 2007-6-12 22:41 | 显示全部楼层
A'是A的共轭转置。A.'是A的普通转置。
当然A如果是实数矩阵,两者一样的

评分

1

查看全部评分

 楼主| 发表于 2007-6-12 22:53 | 显示全部楼层
楼上是说直接写A.'就可以了吗?好像我试过不行。不知道是不是我机器上matlab安装的问题,因为simulink也用不了,本来想省事不编程了,可是唉:@( 卸了重装也不行,郁闷。
发表于 2007-6-12 23:02 | 显示全部楼层
以下是我机子上运行结果。你试试,不行就是你机子的问题了
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
发表于 2007-6-13 07:46 | 显示全部楼层
A'是共扼转置,但A为复数时就不是简单的转置
 楼主| 发表于 2007-6-13 23:42 | 显示全部楼层
试了:@) 好了谢谢。
现在正在搞递推算法......
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 03:07 , Processed in 0.085480 second(s), 26 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表