Skip to main content

 Bộ lọc Kalman - giành cho người mới bắt đầu


Ý tưởng: 

Để thực hiện việc loại bỏ nhiễu trong việc đo lường tín hiệu, ta sử dụng một hệ thống lọc nhiễu hay gọi vắn tắt là bộ lọc (Filter). Các cách lọc nhiễu trong hệ thống căn bản gồm các bộ lọc thông thấp (Lowpass Filters), bộ lọc thông dải (Bandpass Filters), bộ lọc chặn (Notch Filters) và lọc thông cao (Highpass Filters). Dựa trên quan sát sự khác biệt giữa tín hiệu và nhiễu để thực hiện việc tách tín hiệu khỏi nhiễu.



Fig1. Các kiểu bộ lọc cơ bản

Trong trường hợp cơ bản, để đo chất lượng tín hiệu ta dùng chỉ số SNR (Signal-to-Noise ratio) để đánh giá mức độ mạnh giữa tín hiệu so với nhiễu. Nếu giữa tín hiệu và nhiễu có tính độc lập ( Tương quan là 0) thì việc tách tín hiệu là khả thi. Chẳng hạn các tín hiệu ta đo lường thường có một số tần số nhất định, ví dụ các tín hiệu sinh học thường khu trú trong vùng tần số thấp < 25 Hz. Các tín hiệu nhiễu nền thường là nhiễu trắng có phân bố rộng, nhiễu điện lưới có tần số xác định 50/60 Hz. 

Trong trường hợp xuất hiện các loại nhiễu không xác định được, đây là nhiễu xảy ra trong hệ thống do quá trình vận hành, nhiễu từ một nguồn tác động bên ngoài. Khi đó việc dùng các bộ lọc cơ bản sẽ không đủ để ta có thể loại bỏ nhiễu. Lúc này một cách thức lọc thích nghi (adaptive filter) cho phép cập nhật các hệ số điều chỉnh của bộ lọc có thể được áp dụng. Kalman thuộc về nhóm các bộ lọc thích nghi này. 

Ý tưởng của Kalman dựa trên giả định rằng hệ thống là tuyến tính. Và căn cứ vào tính chất đã biết của hệ thống, ta đưa ra dự đoán về tín hiệu ở thời điểm tiếp theo. Sau đó, căn cứ vào sự khác biệt giữa dự đoán và loại tín hiệu đo được trên thực tế để hiệu chỉnh độ chính xác của dự đoán để cải thiện dự đoán cho các lần tiếp theo. 

Như vậy, bộ lọc Kalman mang hình thức lai ghép giữa mô hình học thống kê dựa trên Bayes, tuy nhiên đơn giản hơn Bayes dựa trên thống kê của toàn bộ tập dữ liệu. Mô hình học của Kalman dựa trên một hệ thống đã biết quy luật và nhiễu không xác định quy luật để từ đó đưa ra dự đoán bám theo sự thay đổi chưa định trước. Phương pháp học Kalman chỉ quan tâm tới trạng thái thời điểm trước đó và kiểm nghiệm vị trí hiện tại để điểu chỉnh giá trị của tín hiệu sau khi đo. 

Một cách toán học ta diễn ta mô hình học của Kalman như sau:

                                            x(t ) = F(t )x(t-1)  + B(t )u(t) + v(t)         (1)

Trong đó:

x(t) là vector trạng thái của hệ thống chứa các thành phần của hệ thống mà ta quan tâm. Ví dụ như để theo dõi chuyển động của một vật thì x sẽ bao gồm vị trí, vận tốc. Với một tín hiệu thì nó sẽ gồm giá trị trung bình \mu và độ tin cậy dựa trên giá trị độ lệch biến \sigma.

u(t) là vector chứa các yếu tố điều khiển tác động vào hệ thống và có tính xác định. Ví dụ như lực thắng, tốc độ đánh lái.

F(t) là ma trận mô tả quan hệ giữa trạng thái hiện tại và trạng thái kế tiếp gọi là ma trận chuyển tiếp trạng thái (state transitional matrix). Ví dụ như quan hệ giữa vị trí và tốc độ ở hai thời điểm kế tiếp nhau. 

B(t) là ma trận điều khiển mô tả tác động của yếu tố điều khiển xác định lên hệ thống ở hiện tại. 

v(t) là yếu tố nhiễu phát sinh không xác định được hay còn gọi là nhiễu hệ thống. 

Như vậy ở góc nhìn của điều khiển hệ thống tuyến tính, nếu ta biết quan hệ giữa trạng thái hiện tại và trạng thái kế F, tác động của các lực tác động vào hệ thống B, thì trong điều kiện lý tưởng không có nhiễu v = 0, ta  sẽ tính được tình trạng kế tiếp của hệ thống mà không cần bất cứ đo đạc nào thêm. Hệ thống như vậy được xem là xác định và ta chỉ cần đo một lần duy nhất thời điểm ban đầu, thì ta có thể xác định hệ thống ở các thời điểm tiếp theo. Đây là bài toán tìm quỹ đạo của vật mà ta đã từng biết trong chương trình vật lý phổ thông.

Tuy nhiên, trong trường hợp xuất hiện của nhiễu, việc xác định x(t) ta chỉ có thể căn cứ vào F B đã biết. Ta gọi dự đoán này là dự đoán trước khi hiệu chỉnh, hay dự đoán dựa trên hệ thống. 

                                                x_hat(t|t-1 ) = F(t )x(t-1|t-1)  + B(t )u(t)     (2)

So (2) và (1) rõ ràng ta cần tìm cách ước lượng nhiễu v(t) để có thể mô tả chính xác trạng thái thật của hệ thống. Để làm việc này ta sử dụng thủ thuật dò vết hệ thống. Nghĩa là ta so sánh sự khác biệt giữa dự đoán so với giá trị đo thật để điều chỉnh giá trị nhiễu hệ thống thêm vào cho phương trình (2)

Và do đó giá trị cuối cùng của hệ thống sẽ có dạng:

                    x_hat(t|t ) =  x_hat(t|t-1) + Ky(t) = F(t )x(t-1|t-1)  + B(t )u(t) + Kw(t)       (3)

ở đây K là ma trận Kalman, và w vector độ lệch giữa giá trị đo được so với dự đoán.

Thông thường trong các hệ thống đo lường, ta không thể đo được x(t) trực tiếp. Bản thân x(t) ta đang đo dựa trên mô hình ở phương trình (1) đã chữa nhân tố bất định. Cho nên ta sẽ dùng một hế thống đo thứ hai để ước lượng giá trị x(t). Một cách tổng quát giả sử hệ thống hai đo lường dựa trên một cơ sở khác, ví dụ như ta đang muốn ước lượng vị trí của vật thể. Giờ ta dùng hệ thống thứ hai ước lượng dựa trên sự dịch thời gian giữa sóng phát ra từ anten trên đối tượng so với vị trí quan sát. Rõ ràng hệ thống đo thời gian lệch sẽ dùng đơn vị giây (s) trong khi hệ thống đo ban đầu xác định vị trí vật thể dựa trên đơn vị met (m). Do đó trong trường hợp tổng quát. Hệ thống thứ hai đo lường sẽ có dạng:

                                                        z(t) = Hx(t|t-1) + w(t)     (3')

Ở đây H là ma trận chuyển cơ sở để chuyển đổi từ hệ thống ta muốn quan sát sang hệ thống ta đo lường được. Do đó (3) có thể được viết lại là :

                x_hat(t|t ) =  F(t )x(t-1|t-1)  + B(t )u(t)  + K(z(t) - Hx(t|t-1))         (4)

Làm thế nào để tính được ma trận K?

Lưu ý ở (1) tồn tại một loại nhiễu mô hình hệ thống Q. Ở  (3') tồn tại một loại nhiễu do đo lường R. Do đó K phải hiệu chỉnh để điều chỉnh đáp ứng đối với sự tác động của hai yếu tố nhiễu này tạo ra. Trong các hệ thống tuyến tính ta giả sử các nhiễu này tuân theo phân bố Gaussian có trung bình là 0. Q, R  là hai ma trận hiệp biến để mô tả mức tác động nhiễu. Về mặt thống kê thì 

        Q = cov(vv') =  E[(v - v_bar)(v- v_bar)') và R = cov (ww') = E[(w - w_bar)(ww_bar)')         

Trong đó v_bar, w_bar là giá trị trung bình của nhiễu hệ thống và nhiễu đo lường.   

Để đánh giá sự sai lệch giữa giá trị dự đoán so với mong muốn, Kalman sử dụng một ma trận mới để tích lũy sự sai lệch giữa giá trị dự đoán và giá trị đo lường. 

Gọi P là ma trận dự đoán, thể hiện tương quan giữa giá trị của hệ thống ở thời điểm trước đó và thời điểm hiện tại. 

                                        P(t|t-1) = cov([x(t) - x(t|t-1)][x(t) - x(t|t-1)]') 

                                                    = F(t)P(t-1|t-1)F(t)' + Q(t)

P(t|t-1) là tương quan giữa giá trị dự đoán so với thực tế. Ta không biết thực tế nên giá trị tính toán của nó phụ thuộc vào giá trị khởi tạo ban đầu. P(t-1|t-1) giá trị tương quan đã hiệu chỉnh ở thời điểm trước. Một lần nữa ta không thể biết Q(t) nên ta sẽ dùng giá trị đo được để hiệu chỉnh.

                                    P(t|t) = P(t|t-1) - K(t)H(t)P(t|t-1)          (5)

Ở đây ma trận Kalman ở (4) và (5) được xác định dựa trên giá trị P ở thời điểm trước đó như sau:

                        K(t) = P(t|t-1)H'(t)/(H(t)P(t|t-1)H'(t) + R(t))                    (6)

Ta thấy mẫu số S(t) = H(t)P(t|t-1)H'(t) + R(t)  chính là ma trận dự đoán sai lệch so với giá trị đo được ở thang đo s. 

Công thức (6) cho thấy tỉ số giữa dự đoán trong không gian m so với không gian s. Cho nên ma trận Kalman được gọi độ độ lợi của bộ lọc




Comments

Popular posts from this blog

Arduino Code for test Heart Rate 7 Click

Heart Rate 7 Click is the newest module from MikroE which uses VEM8080 photodetector has a wide range spectrum from 300nm - 1000nm. To control and acquisition data, AFE4404 from TI inc. is adopted. This chip permits control 3 LED channels, and sample heart rate default 100 SPS.  A 22-bit ADC permit collecting very small changed voltage from a PD sensor. In this example we config Arduino Mega  2650 as below: Pin 4 for RST PIN 5 CLK PIN 6 INT PIN 20 SDA PIN 21 SCL Config registers follow the default of AFE4404 datasheet Page 27 with some minor changes. 1. Config Internal Clock through  AFE_CONTROL2 register addr.: 0x23 value: 0x104218  // setup all dynamic power to save energy 2. Control LED2 current through AFE_LEDCNTRL register addr: 0x22 value: 0x000100 3. Read data using PROG_TG_EN signal through AFE_TIA_GAIN register addr: 0x21 value: 0x000125 Time to start and end of PROG_TG setup through two registers: AFE_PROG_TG_STC register (...
Serial Port Profile with CC2650 and CC2640R2 (part 2) Modified to transfer ADC values up to PC Step 1:  Include ACD.h library #include <ti/drivers/ADC.h>  Step 2: Insert period event in order to update ADC read values from the potentiometer #define SSSS_ICALL_EVT                         ICALL_MSG_EVENT_ID // Event_Id_31 #define SSSS_QUEUE_EVT                         UTIL_QUEUE_EVENT_ID // Event_Id_30 #define SSSS_PERIODIC_EVT                      Event_Id_00 //SOLUTION // Bitwise OR of all events to pend on #define SSSS_ALL_EVENTS                        (SSSS_ICALL_EVT        | \                                ...
Programing on Launchpad MSP432 Interrupt button Prerequisite: Hardware: Launchpad MSP432P401 IDE software: IAR/CCS/KEIL Tool: Simplelink SDK MSP432 v3.20.00.06 Description: Simple program using two buttons connected to P1.0  and P1.4 to turn on or turn off LED on P1.0 Step1. Open example of Simplelink SDK without RTOS Step2. modify the code as below Step3. Compile and debug the program /* DriverLib Includes */ #include <ti/devices/msp432p4xx/driverlib/driverlib.h> /* Standard Includes */ #include <stdint.h> #include <stdbool.h> volatile uint32_t SW1, SW2; //semaphores void VectorButtons_Init(void){   __disable_interrupt;   //DisableInterrupts();          // set the I bit during initialization   SW1 = 0; SW2 = 0;             // clear semaphores   P1->SEL0 &= ~0x12;   P1->SEL1 &= ~0x12;      P1->D...