---
# System prepended metadata

title: Kalman filter for angular velocity estimation

---

卡爾曼濾波器（Kalman Filter, KF）是一種最優估計技術，可用於根據角度測量來估算速度。在這裡，我們假設有一個系統，其中角度是可測量的，而速度是我們要估算的隱藏狀態。

---

## **問題定義**
我們考慮一個一維系統：
- 狀態向量：$x = \begin{bmatrix} \theta \\ \dot{\theta} \end{bmatrix}$（角度和角速度）
- 狀態轉移模型：
$$
  x_{k+1} = A x_k + w_k
$$
  其中：
  $$
  A = \begin{bmatrix} 1 & T \\ 0 & 1 \end{bmatrix}
  $$
  $T$ 是離散時間間隔，$w_k$ 是過程噪聲。
- 測量方程：
  $$
  z_k = H x_k + v_k
  $$
  其中：
  $$
  H = \begin{bmatrix} 1 & 0 \end{bmatrix}
  $$
  $z_k$ 是測得的角度，$v_k$ 是測量噪聲。

---

## **Matlab Code**
以下是使用卡爾曼濾波器來估計角速度的 Matlab 代碼：

```matlab
clc; clear; close all;

% 參數設定
T = 0.1;  % 取樣時間
A = [1 T; 0 1];  % 狀態轉移矩陣
H = [1 0];  % 觀測矩陣
Q = [1e-4 0; 0 1e-2];  % 過程噪聲協方差
R = 1e-2;  % 測量噪聲協方差
P = eye(2);  % 初始協方差矩陣
x_est = [0; 0];  % 初始狀態 [theta; theta_dot]

% 生成模擬數據
N = 100;  % 時間步數
theta_true = zeros(1, N);
theta_dot_true = zeros(1, N);
z = zeros(1, N);

for k = 2:N
    theta_dot_true(k) = 0.1 * sin(0.1 * k);  % 真實速度變化
    theta_true(k) = theta_true(k-1) + T * theta_dot_true(k);
    z(k) = theta_true(k) + sqrt(R) * randn;  % 加入測量噪聲
end

% 卡爾曼濾波估計
x_est_hist = zeros(2, N);

for k = 1:N
    % 預測步驟
    x_pred = A * x_est;
    P_pred = A * P * A' + Q;
    
    % 更新步驟
    K = P_pred * H' / (H * P_pred * H' + R);
    x_est = x_pred + K * (z(k) - H * x_pred);
    P = (eye(2) - K * H) * P_pred;
    
    % 儲存結果
    x_est_hist(:, k) = x_est;
end

% 畫圖
figure;
subplot(2,1,1);
plot(1:N, theta_true, 'g', 1:N, z, 'b.', 1:N, x_est_hist(1,:), 'r');
legend('True Angle', 'Measured Angle', 'Estimated Angle');
xlabel('Time Step'); ylabel('Angle (rad)');

subplot(2,1,2);
plot(1:N, theta_dot_true, 'g', 1:N, x_est_hist(2,:), 'r');
legend('True Velocity', 'Estimated Velocity');
xlabel('Time Step'); ylabel('Angular Velocity (rad/s)');
```

---

## **代碼解釋**
1. **系統初始化**：
   - 設定狀態轉移矩陣 $A$、觀測矩陣 $H$ 及協方差矩陣 $Q, R$。
   - 初始化狀態 $x_{est}$ 和協方差 $P$。

2. **模擬真實數據**：
   - 使用 $\theta_{dot}(k) = 0.1 \sin(0.1 k)$ 產生角速度變化。
   - 通過數值積分獲得角度，並添加測量噪聲。

3. **卡爾曼濾波迴圈**：
   - **預測步驟**：
     - 計算預測狀態  $x_{pred} = A x_{est}$。
     - 更新誤差協方差矩陣 $P_{pred}$。
   - **更新步驟**：
     - 計算卡爾曼增益 $K$。
     - 更新狀態估計 $x_{est}$。
     - 更新誤差協方差矩陣 $P$。

4. **繪圖**：
   - 比較真實角度、測量角度與估計角度。
   - 比較真實角速度與估計角速度。

這樣就能使用卡爾曼濾波器來估算角速度，並在有噪聲的情況下得到較準確的結果！

## 限制式
### 連續form
$$I\ddot{\alpha}<0.6$$$$I\ddot{\beta}<0.6$$
### Discrete form
$$I \frac{\alpha_{k+1} - 2\alpha_k + \alpha_{k-1}}{T^2} < 0.6$$$$I \frac{\beta_{k+1} - 2\beta_k + \beta_{k-1}}{T^2} < 0.6$$
