21ic电子技术开发论坛 › 单片机与嵌入式系统 › 侃单片机论坛 › 授之以渔： 卡尔曼滤波器 ....大泻蜜 ................. ...
12345678910... 14下一页
返回列表 发新帖我要提问

授之以渔： 卡尔曼滤波器 ....大泻蜜 .........................

[复制链接]
115005|281
Cortex-M0 发表于 2011-12-8 18:15 | 显示全部楼层
继续听highgear老师上课~~~
回复

收藏85 举报

 楼主| highgear 发表于 2011-12-8 21:46 | 显示全部楼层
第一个例子中的代码


复制
  1.  
    

  2. Ts = 1;                    采样时间
    

  3. s = tf('s');
    

  4. sysc = 100/(60*s +1);      真实系统的传递函数
    

  5. sysd = c2d(sysc, Ts);
    


  6. 

  7. t = 0:TsTs*300);
    

  8. u = ones(1, length(t));     系统输入
    

  9. ys = lsim(sysd, u, t, 0);   系统输出
    

  10. ys = ys + 10*rand(size(ys));  测量量 = 系统输出 + 噪声
    


  11. 

  12. A = [1, 0; 0, 1];
    

  13. B = [1; 0];
    

  14. C = [1, 0];
    

  15. N = length(A);          系统维数
    

  16. Len = length(u);
    


  17. 

  18. yout = zeros(Len, 1);   滤波器输出
    

  19. Xh = zeros(N, 1);       状态变量
    

  20. P = eye(N);              
    

  21. Q = 2 * eye(N);         系统噪声
    

  22. R = 50;                 测量噪声
    

  23. for i = 1 : Len
    

  24.     Xh = A*Xh + B*u(i);
    

  25.     P = A*P*A' + Q;
    

  26.     K = P*C'* inv(C*P*C' + R);
    

  27.     Xh = Xh + K*(ys(i) - C*Xh);
    

  28.     P = P - K*C*P;
    

  29.     yout(i) = C*Xh;
    

  30. end
    

  31. Ts = 1;                    采样时间
    

  32. s = tf('s');
    

  33. sysc = 100/(60*s +1);      真实系统的传递函数
    

  34. sysd = c2d(sysc, Ts);
    

  35.  
    

  36. t = 0:Ts:(Ts*300);
    

  37. u = ones(1, length(t));     系统输入
    

  38. ys = lsim(sysd, u, t, 0);   系统输出
    

  39. ys = ys + 10*rand(size(ys));  测量量 = 系统输出 + 噪声
    

  40.  
    

  41. A = [1, 0; 0, 1];
    

  42. B = [1; 0];
    

  43. C = [1, 0];
    

  44. N = length(A);          系统维数
    

  45. Len = length(u);
    

  46.  
    

  47. yout = zeros(Len, 1);   滤波器输出
    

  48. Xh = zeros(N, 1);       状态变量
    

  49. P = eye(N);              
    

  50. Q = 2 * eye(N);         系统噪声
    

  51. R = 50;                 测量噪声
    

  52. for i = 1 : Len
    

  53.     Xh = A*Xh + B*u(i);
    

  54.     P = A*P*A' + Q;
    

  55.     K = P*C'* inv(C*P*C' + R);
    

  56.     Xh = Xh + K*(ys(i) - C*Xh);
    

  57.     P = P - K*C*P;
    

  58.     yout(i) = C*Xh;
    

  59. end
    


  60. 

  61. plot(t, ys, t, yout);
    


  62.

, t, yout);
回复

收藏85 举报

 楼主| highgear 发表于 2011-12-8 21:50 | 显示全部楼层
使用最小二乘法例子的代码：
复制
  1.  
    

  2.  
    

  3. Ts = 1;               采样时间
    

  4. s = tf('s');
    

  5. sysc = 100/(60*s +1);   真实系统的传递函数
    

  6. sysd = c2d(sysc, Ts);
    

  7.  
    

  8. t = 0:Ts:(Ts*300);
    

  9. u = ones(1, length(t));   系统输入
    

  10. ys = lsim(sysd, u, t, 0); 系统输出
    

  11. ys = ys + 10*rand(size(ys)); 测量量 = 系统输出 + 噪声
    

  12.  
    

  13. [numd, dend] = LeastSquare(u, ys, 3);   最小二乘法获取预测系统模型
    

  14. [A, B, C, D] = tf2ss(numd, dend);       变换到状态空间形式
    

  15. Len = length(u);
    

  16. N = length(A);         系统维数 
    

  17.  
    

  18. yout = zeros(Len, 1);  滤波器输出
    

  19. Xh = zeros(N, 1);      状态变量
    

  20. P = eye(N);
    

  21. Q = 0.02 * eye(N);     系统噪声
    

  22. R = 50;                测量噪声
    

  23. for i = 1 : Len
    

  24.     Xh = A*Xh + B*u(i);
    

  25.     P = A*P*A' + Q;
    

  26.  
    

  27.     K = P*C'* inv(C*P*C' + R);
    

  28.     Xh = Xh + K*(ys(i) - C*Xh);
    

  29.     P = P - K*C*P;
    

  30.     yout(i) = C*Xh;
    

  31. end
    

  32.  
    

  33. plot(t, ys, t, yout);


超过 100 楼则公布如何通过输入输出获取系统模型的最小二乘法程序
回复

收藏85 举报

Cortex-M0 发表于 2011-12-9 08:27 | 显示全部楼层
顶~~~

继续听highgear老师上课~~~
回复

收藏85 举报

Cortex-M0 发表于 2011-12-9 08:28 | 显示全部楼层
有钱出钱
回复

收藏85 举报

Cortex-M0 发表于 2011-12-9 08:28 | 显示全部楼层
有力出力
回复

收藏85 举报

Cortex-M0 发表于 2011-12-9 08:29 | 显示全部楼层
没本事的多灌灌水
回复

收藏85 举报

Cortex-M0 发表于 2011-12-9 08:29 | 显示全部楼层
多顶顶highgear老师之贴
回复

收藏85 举报

Cortex-M0 发表于 2011-12-9 08:30 | 显示全部楼层
多听听highgear老师上课
回复

收藏85 举报

lxyppc 发表于 2011-12-9 08:30 | 显示全部楼层
100楼太容易了
回复

收藏85 举报

Cortex-M0 发表于 2011-12-9 10:18 | 显示全部楼层
LS盆友正解~~~

这该上1000楼才对~~~ :lol
回复

收藏85 举报

DownCloud 发表于 2011-12-9 10:50 | 显示全部楼层
100楼，必须的。
回复

收藏85 举报

skyfight 发表于 2011-12-9 11:01 | 显示全部楼层
不错，支持一下
回复

收藏85 举报

DreamNk 发表于 2011-12-9 11:14 | 显示全部楼层
:lol
回复

收藏85 举报

pidcontrol 发表于 2011-12-9 12:44 | 显示全部楼层
精彩精彩
回复

收藏85 举报

xxmmtt 发表于 2011-12-9 13:46 | 显示全部楼层
:D
回复

收藏85 举报

Periodic 发表于 2011-12-9 15:02 | 显示全部楼层
高深 目前还是听不大 懂  功底不太好
回复

收藏85 举报

Periodic 发表于 2011-12-9 15:03 | 显示全部楼层
呵呵 没有钱出力
回复

收藏85 举报

Periodic 发表于 2011-12-9 15:04 | 显示全部楼层
21家 上 谈技术的 热情 大不如前了啊 !!现在 才 79楼
回复

收藏85 举报

Cortex-M0 发表于 2011-12-9 15:28 | 显示全部楼层
80楼 :victory:
回复

收藏85 举报

下一页 »
12345678910... 14下一页
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 注册

本版积分规则

快速回复 在线客服 返回列表 返回顶部