4545#include <fcntl.h>
4646#include <poll.h>
4747#include <errno.h>
48+ #include <float.h>
49+ #include <math.h>
4850
4951#include <nuttx/sensors/cxd5602pwbimu.h>
5052#include <MadgwickAHRS.h>
5961#define DEFAULT_GYROSCOPERANGE (1000)
6062#define DEFAULT_ACCELRANGE (4)
6163
64+ #define RAD2DEG (x ) ((x) * 180.0 / M_PI)
65+
6266/****************************************************************************
6367 * Private Functions
6468 ****************************************************************************/
6569
70+ /****************************************************************************
71+ * Private Data
72+ ****************************************************************************/
73+
74+ static uint32_t s_prev_time ;
75+ static uint64_t s_prev_unroll_time ;
76+
6677/** start_sensing()
6778 * @rate [in]: 15, 30, 60, 120, 240, 480, 960, 1920
6879 * @adrange [in]: 2, 4, 8, 16
@@ -118,21 +129,6 @@ static int start_sensing(int rate, int adrange, int gdrange, int nfifos)
118129 return fd ;
119130}
120131
121- static int drop_50msdata (int fd , int samprate , cxd5602pwbimu_data_t * imu )
122- {
123- int cnt = samprate / 20 ; /* data size of 50ms */
124-
125- if (cnt == 0 ) cnt = 1 ;
126-
127- while (cnt )
128- {
129- read (fd , imu , sizeof (imu [0 ]));
130- cnt -- ;
131- }
132-
133- return 0 ;
134- }
135-
136132static int read_imudata (int fd , cxd5602pwbimu_data_t * imudata )
137133{
138134 char c ;
@@ -181,6 +177,107 @@ static int read_imudata(int fd, cxd5602pwbimu_data_t *imudata)
181177 return ret ;
182178}
183179
180+ static int drop_50msdata (int fd , int samprate , cxd5602pwbimu_data_t * imu )
181+ {
182+ int cnt = samprate / 20 ; /* data size of 50ms */
183+
184+ if (cnt == 0 ) cnt = 1 ;
185+
186+ while (cnt )
187+ {
188+ read_imudata (fd , imu );
189+ cnt -- ;
190+ }
191+
192+ return 0 ;
193+ }
194+
195+ #ifdef CONFIG_EXAMPLES_AHRS_PWBIMU_EXEC_GYROBIAS_ESTIMATION
196+ static void max_min_check (float val , float * min_val , float * max_val )
197+ {
198+ * min_val = val < * min_val ? val : * min_val ;
199+ * max_val = val > * max_val ? val : * max_val ;
200+ }
201+
202+ static int gyrobias_estimation (int fd , int samprate , int estimation_time ,
203+ double * gyrobias , cxd5602pwbimu_data_t * imu )
204+ {
205+ int cnt = estimation_time * samprate ; /* data size of estimation time */
206+ int cnt_dwn = cnt ;
207+ int elp_cnt = 0 ;
208+ float min_gx = FLT_MAX , max_gx = - FLT_MAX ;
209+ float min_gy = FLT_MAX , max_gy = - FLT_MAX ;
210+ float min_gz = FLT_MAX , max_gz = - FLT_MAX ;
211+ float min_ax = FLT_MAX , max_ax = - FLT_MAX ;
212+ float min_ay = FLT_MAX , max_ay = - FLT_MAX ;
213+ float min_az = FLT_MAX , max_az = - FLT_MAX ;
214+ double accelave [3 ] = {0.0 };
215+ int i = 0 ;
216+
217+ printf ("Start gyro-bias estimation.\n" );
218+
219+ while (cnt_dwn )
220+ {
221+ read_imudata (fd , imu );
222+
223+ gyrobias [0 ] += imu -> gx ;
224+ gyrobias [1 ] += imu -> gy ;
225+ gyrobias [2 ] += imu -> gz ;
226+
227+ accelave [0 ] += imu -> ax ;
228+ accelave [1 ] += imu -> ay ;
229+ accelave [2 ] += imu -> az ;
230+
231+ max_min_check (imu -> gx , & min_gx , & max_gx );
232+ max_min_check (imu -> gy , & min_gy , & max_gy );
233+ max_min_check (imu -> gz , & min_gz , & max_gz );
234+
235+ max_min_check (imu -> ax , & min_ax , & max_ax );
236+ max_min_check (imu -> ay , & min_ay , & max_ay );
237+ max_min_check (imu -> az , & min_az , & max_az );
238+
239+ cnt_dwn -- ;
240+
241+ if (!(cnt_dwn % samprate ))
242+ {
243+ printf ("Elapsed time : %03d seconds.\n" , ++ elp_cnt );
244+ }
245+ }
246+
247+ printf ("Finish gyro-bias estimation.\n" );
248+
249+ for (i = 0 ; i < 3 ; i ++ )
250+ {
251+ gyrobias [i ] /= (double )cnt ;
252+ accelave [i ] /= (double )cnt ;
253+ }
254+
255+ printf ("gx: ave=%12.3e, min=%12.3e, max=%12.3e\n" , gyrobias [0 ], min_gx , max_gx );
256+ printf ("gy: ave=%12.3e, min=%12.3e, max=%12.3e\n" , gyrobias [1 ], min_gy , max_gy );
257+ printf ("gz: ave=%12.3e, min=%12.3e, max=%12.3e\n" , gyrobias [2 ], min_gz , max_gz );
258+
259+ printf ("ax: ave=%12.3e, min=%12.3e, max=%12.3e\n" , accelave [0 ], min_ax , max_ax );
260+ printf ("ay: ave=%12.3e, min=%12.3e, max=%12.3e\n" , accelave [1 ], min_ay , max_ay );
261+ printf ("az: ave=%12.3e, min=%12.3e, max=%12.3e\n" , accelave [2 ], min_az , max_az );
262+
263+ return 0 ;
264+ }
265+ #endif
266+
267+ static uint64_t generate_unroll_timestamp (uint32_t curr_time )
268+ {
269+ s_prev_unroll_time += (uint64_t )(curr_time - s_prev_time );
270+ s_prev_time = curr_time ;
271+ return s_prev_unroll_time ;
272+ }
273+
274+ static void set_initial_posture (int fd , struct ahrs_out_s * inst ,
275+ cxd5602pwbimu_data_t * imu )
276+ {
277+ read_imudata (fd , imu );
278+ setPostureByAccel (inst , imu -> ax , imu -> ay , imu -> az , 0.f );
279+ }
280+
184281/****************************************************************************
185282 * Public Functions
186283 ****************************************************************************/
@@ -191,9 +288,13 @@ int main(int argc, FAR char *argv[])
191288 struct ahrs_out_s ahrs ;
192289 cxd5602pwbimu_data_t imu ;
193290 float e [3 ];
291+ double gyrobias [3 ] = {0.0 };
194292 unsigned int cnt = 0 ;
195293 int print_hex = 0 ;
196294
295+ s_prev_time = 0 ;
296+ s_prev_unroll_time = 0 ;
297+
197298 if (argc == 2 && argv [1 ][0 ] == 'h' )
198299 {
199300 print_hex = 1 ;
@@ -213,25 +314,43 @@ int main(int argc, FAR char *argv[])
213314 * stationary for several seconds, use that value as the BIAS
214315 * value, and subtract that value from each gyro sample data.
215316 */
317+ #ifdef CONFIG_EXAMPLES_AHRS_PWBIMU_EXEC_GYROBIAS_ESTIMATION
318+ /* This is a example of gyro bias estimation.
319+ * This function estimate the gyro bias and the angular velocity
320+ * caused by the Earth's rotation together as offset values.
321+ * Please maintain a stationary state as much as possible
322+ * during the estimation.
323+ */
324+ gyrobias_estimation (fd , DEFAULT_SAMPLERATE ,
325+ CONFIG_EXAMPLES_AHRS_PWBIMU_GYROBIAS_ESTIMATION_TIME ,
326+ gyrobias , & imu );
327+ #endif
328+
329+ set_initial_posture (fd , & ahrs , & imu );
216330
217331 while (read_imudata (fd , & imu ))
218332 {
219- MadgwickAHRSupdateIMU (& ahrs , imu .gx , imu .gy , imu .gz ,
220- imu .ax , imu .ay , imu .az );
333+ MadgwickAHRSupdateIMU (& ahrs ,
334+ imu .gx - (float )gyrobias [0 ],
335+ imu .gy - (float )gyrobias [1 ],
336+ imu .gz - (float )gyrobias [2 ],
337+ imu .ax , imu .ay , imu .az );
221338
222339 quaternion2euler (ahrs .q , e );
223340
224341 if (++ cnt >= 20 ) /* Decimate data 1920Hz / 20 = 96Hz */
225342 {
226343 if (print_hex )
227344 {
228- printf ("%08x,%08x,%08x\n" , * (unsigned int * )& e [0 ],
229- * (unsigned int * )& e [1 ],
230- * (unsigned int * )& e [2 ]);
345+ printf ("%08x,%08x,%08x,%08x\n" , * (unsigned int * )& ahrs .q [0 ],
346+ * (unsigned int * )& ahrs .q [1 ],
347+ * (unsigned int * )& ahrs .q [2 ],
348+ * (unsigned int * )& ahrs .q [3 ]);
231349 }
232350 else
233351 {
234- printf ("R:%0.1f, P:%0.1f, Y:%0.1f\n" , e [0 ], e [1 ], e [2 ]);
352+ uint64_t unroll_time = generate_unroll_timestamp (imu .timestamp );
353+ printf ("T:%0.2f, R:%0.2f, P:%0.2f, Y:%0.2f\n" , (float )(unroll_time /19200000.0f ), RAD2DEG (e [0 ]), RAD2DEG (e [1 ]), RAD2DEG (e [2 ]));
235354 }
236355
237356 cnt = 0 ;
0 commit comments