footPDR.c 24 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <stdlib.h>
  4. #include <math.h>
  5. #include <stdint.h>
  6. #define ZUPT_threshold 0.81
  7. #define SIGMA 0.01
  8. #define SIGMA_V 0.01
  9. #define PI 3.1416
  10. #define RUN 1
  11. #define STAND 0
  12. extern void sen_step_to_host(uint8_t, uint16_t);
  13. //当地的重力加速度
  14. float g = 9.788f;
  15. float dt = 0.01f;
  16. float P[81], acc_n[3];
  17. float Temporary_array1[9], Temporary_array2[9];
  18. float K[27], P_prev[81], delta_x[9];
  19. float C[9], C_prev[9];
  20. float vel_n[3], pos_n[3];
  21. float last_pos_n[3];
  22. float pos_offset[3];
  23. int frame_index = 0;
  24. int stand_num = 0;
  25. float gyr_norm_window[10];
  26. float gyr_z_window[10];
  27. float gyr_extreme[6];
  28. float gyr_mean[3];
  29. float num_peak;
  30. float acc_mean[3];
  31. float gyrBias[3];
  32. float last_gyr[3];
  33. float last_acc[3];
  34. float last_vel_n[3];
  35. float accSum;
  36. float accSize;
  37. int press_data[10];
  38. int ZUPT_STATUS;
  39. int HAS_RESULT;
  40. int IS_DOWN;
  41. int press_index;
  42. float acc_shoes[3];
  43. float last_acc_shoes[3];
  44. int last_move_index;
  45. int down_pass1;
  46. int down_pass2;
  47. float theta;
  48. //last_stage:0 为 走路状态;
  49. //last_stage:1 为 静止状态
  50. int last_stage;
  51. uint8_t step_count;
  52. uint16_t step_distance;
  53. int step_time = 0;
  54. //
  55. short pos_offset_max = 0;
  56. short pos_offset_min = 0;
  57. void cal_step_data(void)
  58. {
  59. static int step_count_sum = 0;
  60. step_count = 1;
  61. step_count_sum += step_count;
  62. float step_length = sqrt((last_pos_n[0] - pos_n[0])*(last_pos_n[0] - pos_n[0]) + (last_pos_n[1] - pos_n[1])*(last_pos_n[1] - pos_n[1]));
  63. if(step_length > 5.0f)
  64. {
  65. step_length = 1.2f;
  66. }
  67. step_distance = (uint16_t)(step_length * 100.0f);
  68. sen_step_to_host(step_count, step_distance);
  69. }
  70. uint8_t get_step_count(void)
  71. {
  72. return step_count;
  73. }
  74. uint16_t get_step_length(void)
  75. {
  76. return step_distance;
  77. }
  78. void invert3x3(float * src, float * dst)
  79. {
  80. float det;
  81. /* Compute adjoint: */
  82. dst[0] = +src[4] * src[8] - src[5] * src[7];
  83. dst[1] = -src[1] * src[8] + src[2] * src[7];
  84. dst[2] = +src[1] * src[5] - src[2] * src[4];
  85. dst[3] = -src[3] * src[8] + src[5] * src[6];
  86. dst[4] = +src[0] * src[8] - src[2] * src[6];
  87. dst[5] = -src[0] * src[5] + src[2] * src[3];
  88. dst[6] = +src[3] * src[7] - src[4] * src[6];
  89. dst[7] = -src[0] * src[7] + src[1] * src[6];
  90. dst[8] = +src[0] * src[4] - src[1] * src[3];
  91. /* Compute determinant: */
  92. det = src[0] * dst[0] + src[1] * dst[3] + src[2] * dst[6];
  93. /* Multiply adjoint with reciprocal of determinant: */
  94. det = 1.0f / det;
  95. dst[0] *= det;
  96. dst[1] *= det;
  97. dst[2] *= det;
  98. dst[3] *= det;
  99. dst[4] *= det;
  100. dst[5] *= det;
  101. dst[6] *= det;
  102. dst[7] *= det;
  103. dst[8] *= det;
  104. }
  105. void multiply3x3(float *a, float *b, float *dst)
  106. {
  107. dst[0] = a[0] * b[0] + a[1] * b[3] + a[2] * b[6];
  108. dst[1] = a[0] * b[1] + a[1] * b[4] + a[2] * b[7];
  109. dst[2] = a[0] * b[2] + a[1] * b[5] + a[2] * b[8];
  110. dst[3] = a[3] * b[0] + a[4] * b[3] + a[5] * b[6];
  111. dst[4] = a[3] * b[1] + a[4] * b[4] + a[5] * b[7];
  112. dst[5] = a[3] * b[2] + a[4] * b[5] + a[5] * b[8];
  113. dst[6] = a[6] * b[0] + a[7] * b[3] + a[8] * b[6];
  114. dst[7] = a[6] * b[1] + a[7] * b[4] + a[8] * b[7];
  115. dst[8] = a[6] * b[2] + a[7] * b[5] + a[8] * b[8];
  116. }
  117. void multiply3x1(float *a, float *b, float *dst)
  118. {
  119. dst[0] = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
  120. dst[1] = a[3] * b[0] + a[4] * b[1] + a[5] * b[2];
  121. dst[2] = a[6] * b[0] + a[7] * b[1] + a[8] * b[2];
  122. }
  123. void init_attitude_matrix(float *C, float *acc)
  124. {
  125. float pitch = asin(-acc[0] / g);
  126. float roll = atan2(acc[1], acc[2]);
  127. float pitch_sin = sin(pitch);
  128. float pitch_cos = cos(pitch);
  129. float roll_sin = sin(roll);
  130. float roll_cos = cos(roll);
  131. C[0] = cos(pitch);
  132. C[1] = pitch_sin * roll_sin;
  133. C[2] = pitch_sin * roll_cos;
  134. C[4] = roll_cos;
  135. C[5] = -roll_sin;
  136. C[6] = -pitch_sin;
  137. C[7] = pitch_cos * roll_sin;
  138. C[8] = pitch_cos * roll_cos;
  139. }
  140. void reset_yaw_C(float *C)
  141. {
  142. float pitch = asin(-C[6]);
  143. float roll = atan2(C[7], C[8]);
  144. float pitch_sin = sin(pitch);
  145. float pitch_cos = cos(pitch);
  146. float roll_sin = sin(roll);
  147. float roll_cos = cos(roll);
  148. C[0] = pitch_cos;
  149. C[1] = pitch_sin * roll_sin;
  150. C[2] = pitch_sin * roll_cos;
  151. C[3] = 0.0;
  152. C[4] = roll_cos;
  153. C[5] = -roll_sin;
  154. }
  155. // F * P * F'
  156. void State_covariance_matrix_update(float *P, float *acc_n)
  157. {
  158. // P2 + P3*dt,
  159. P[3] = P[3] + P[6] * dt;
  160. P[4] = P[4] + P[7] * dt;
  161. P[5] = P[5] + P[8] * dt;
  162. P[12] = P[12] + P[15] * dt;
  163. P[13] = P[13] + P[16] * dt;
  164. P[14] = P[14] + P[17] * dt;
  165. P[21] = P[21] + P[24] * dt;
  166. P[22] = P[22] + P[25] * dt;
  167. P[23] = P[23] + P[26] * dt;
  168. //P4 + P7*dt,
  169. P[27] = P[27] + P[54] * dt;
  170. P[28] = P[28] + P[55] * dt;
  171. P[29] = P[29] + P[56] * dt;
  172. P[36] = P[36] + P[63] * dt;
  173. P[37] = P[37] + P[64] * dt;
  174. P[38] = P[38] + P[65] * dt;
  175. P[45] = P[45] + P[72] * dt;
  176. P[46] = P[46] + P[73] * dt;
  177. P[47] = P[47] + P[74] * dt;
  178. // P5 + P8*dt + dt*(P6 + P9*dt)
  179. P[30] = P[30] + P[57] * dt + dt * (P[33] + P[60] * dt);
  180. P[31] = P[31] + P[58] * dt + dt * (P[34] + P[61] * dt);
  181. P[32] = P[32] + P[59] * dt + dt * (P[35] + P[62] * dt);
  182. P[39] = P[39] + P[66] * dt + dt * (P[42] + P[69] * dt);
  183. P[40] = P[40] + P[67] * dt + dt * (P[43] + P[70] * dt);
  184. P[41] = P[41] + P[68] * dt + dt * (P[44] + P[71] * dt);
  185. P[48] = P[48] + P[75] * dt + dt * (P[51] + P[78] * dt);
  186. P[49] = P[49] + P[76] * dt + dt * (P[52] + P[79] * dt);
  187. P[50] = P[50] + P[77] * dt + dt * (P[53] + P[80] * dt);
  188. //P6 + P9*dt + matr*(P4 + P7*dt)
  189. P[33] = P[33] + P[60] * dt + acc_n[2] * dt * P[28] - acc_n[1] * dt * P[29];
  190. P[34] = P[34] + P[61] * dt - acc_n[2] * dt * P[27] + acc_n[0] * dt * P[29];
  191. P[35] = P[35] + P[62] * dt + acc_n[1] * dt * P[27] - acc_n[0] * dt * P[28];
  192. P[42] = P[42] + P[69] * dt + acc_n[2] * dt * P[37] - acc_n[1] * dt * P[38];
  193. P[43] = P[43] + P[70] * dt - acc_n[2] * dt * P[36] + acc_n[0] * dt * P[38];
  194. P[44] = P[44] + P[71] * dt + acc_n[1] * dt * P[36] - acc_n[0] * dt * P[37];
  195. P[51] = P[51] + P[78] * dt + acc_n[2] * dt * P[46] - acc_n[1] * dt * P[47];
  196. P[52] = P[52] + P[79] * dt - acc_n[2] * dt * P[45] + acc_n[0] * dt * P[47];
  197. P[53] = P[53] + P[80] * dt + acc_n[1] * dt * P[45] - acc_n[0] * dt * P[46];
  198. //P7 + P1*matr
  199. P[54] = P[54] + P[9] * acc_n[2] * dt - P[18] * acc_n[1] * dt;
  200. P[55] = P[55] + P[10] * acc_n[2] * dt - P[19] * acc_n[1] * dt;
  201. P[56] = P[56] + P[11] * acc_n[2] * dt - P[20] * acc_n[1] * dt;
  202. P[63] = P[63] - P[0] * acc_n[2] * dt + P[18] * acc_n[0] * dt;
  203. P[64] = P[64] - P[1] * acc_n[2] * dt + P[19] * acc_n[0] * dt;
  204. P[65] = P[65] - P[2] * acc_n[2] * dt + P[20] * acc_n[0] * dt;
  205. P[72] = P[72] + P[0] * acc_n[1] * dt - P[9] * acc_n[0] * dt;
  206. P[73] = P[73] + P[1] * acc_n[1] * dt - P[10] * acc_n[0] * dt;
  207. P[74] = P[74] + P[2] * acc_n[1] * dt - P[11] * acc_n[0] * dt;
  208. //P8 + P2*matr + dt*(P9 + P3*matr),
  209. P[57] = P[57] + dt * P[60] + P[12] * acc_n[2] * dt - P[21] * acc_n[1] * dt;
  210. P[58] = P[58] + dt * P[61] + P[13] * acc_n[2] * dt - P[22] * acc_n[1] * dt;
  211. P[59] = P[59] + dt * P[62] + P[14] * acc_n[2] * dt - P[23] * acc_n[1] * dt;
  212. P[66] = P[66] + dt * P[69] - P[3] * acc_n[2] * dt + P[21] * acc_n[0] * dt;
  213. P[67] = P[67] + dt * P[70] - P[4] * acc_n[2] * dt + P[22] * acc_n[0] * dt;
  214. P[68] = P[68] + dt * P[71] - P[5] * acc_n[2] * dt + P[23] * acc_n[0] * dt;
  215. P[75] = P[75] + dt * P[78] + P[3] * acc_n[1] * dt - P[12] * acc_n[0] * dt;
  216. P[76] = P[76] + dt * P[79] + P[4] * acc_n[1] * dt - P[13] * acc_n[0] * dt;
  217. P[77] = P[77] + dt * P[80] + P[5] * acc_n[1] * dt - P[14] * acc_n[0] * dt;
  218. // P9 + P3*matr + matr*(P7 + P1*matr)
  219. P[60] = P[60] + P[15] * acc_n[2] * dt - P[24] * acc_n[1] * dt + acc_n[2] * dt * P[55] - acc_n[1] * dt * P[56];
  220. P[61] = P[61] + P[16] * acc_n[2] * dt - P[25] * acc_n[1] * dt - acc_n[2] * dt * P[54] + acc_n[0] * dt * P[56];
  221. P[62] = P[62] + P[17] * acc_n[2] * dt - P[26] * acc_n[1] * dt + acc_n[1] * dt * P[54] - acc_n[0] * dt * P[55];
  222. P[69] = P[69] - P[6] * acc_n[2] * dt + P[24] * acc_n[0] * dt + acc_n[2] * dt * P[64] - acc_n[1] * dt * P[65];
  223. P[70] = P[70] - P[7] * acc_n[2] * dt + P[25] * acc_n[0] * dt - acc_n[2] * dt * P[63] + acc_n[0] * dt * P[65];
  224. P[71] = P[71] - P[8] * acc_n[2] * dt + P[26] * acc_n[0] * dt + acc_n[1] * dt * P[63] - acc_n[0] * dt * P[64];
  225. P[78] = P[78] + P[6] * acc_n[1] * dt - P[15] * acc_n[0] * dt + acc_n[2] * dt * P[73] - acc_n[1] * dt * P[74];
  226. P[79] = P[79] + P[7] * acc_n[1] * dt - P[16] * acc_n[0] * dt - acc_n[2] * dt * P[72] + acc_n[0] * dt * P[74];
  227. P[80] = P[80] + P[8] * acc_n[1] * dt - P[17] * acc_n[0] * dt + acc_n[1] * dt * P[72] - acc_n[0] * dt * P[73];
  228. //P3 + P1 * matr
  229. P[6] = P[6] + P[1] * acc_n[2] * dt - P[2] * acc_n[1] * dt;
  230. P[7] = P[7] - P[0] * acc_n[2] * dt + P[2] * acc_n[0] * dt;
  231. P[8] = P[8] + P[0] * acc_n[1] * dt - P[1] * acc_n[0] * dt;
  232. P[15] = P[15] + P[10] * acc_n[2] * dt - P[11] * acc_n[1] * dt;
  233. P[16] = P[16] - P[9] * acc_n[2] * dt + P[11] * acc_n[0] * dt;
  234. P[17] = P[17] + P[9] * acc_n[1] * dt - P[10] * acc_n[0] * dt;
  235. P[24] = P[24] + P[19] * acc_n[2] * dt - P[20] * acc_n[1] * dt;
  236. P[25] = P[25] - P[18] * acc_n[2] * dt + P[20] * acc_n[0] * dt;
  237. P[26] = P[26] + P[18] * acc_n[1] * dt - P[19] * acc_n[0] * dt;
  238. float noise = SIGMA * SIGMA * dt *dt;
  239. for (int i = 0; i < 9; i++)
  240. {
  241. P[i * 9 + i] += noise;
  242. }
  243. }
  244. void Kalfman_gain(float *P, float *Temporary_array, float *Temporary_array1, float *K)
  245. {
  246. Temporary_array[0] = P[60] + SIGMA_V * SIGMA_V * 0.01f;
  247. Temporary_array[1] = P[61];
  248. Temporary_array[2] = P[62];
  249. Temporary_array[3] = P[69];
  250. Temporary_array[4] = P[70] + SIGMA_V * SIGMA_V * 0.01f;
  251. Temporary_array[5] = P[71];
  252. Temporary_array[6] = P[78];
  253. Temporary_array[7] = P[79];
  254. Temporary_array[8] = P[80] + SIGMA_V * SIGMA_V * 0.01f;
  255. invert3x3(Temporary_array, Temporary_array1);
  256. memcpy(Temporary_array, Temporary_array1, 9 * sizeof(float));
  257. K[0] = P[6] * Temporary_array[0] + P[7] * Temporary_array[3] + P[8] * Temporary_array[6];
  258. K[1] = P[6] * Temporary_array[1] + P[7] * Temporary_array[4] + P[8] * Temporary_array[7];
  259. K[2] = P[6] * Temporary_array[2] + P[7] * Temporary_array[5] + P[8] * Temporary_array[8];
  260. K[3] = P[15] * Temporary_array[0] + P[16] * Temporary_array[3] + P[17] * Temporary_array[6];
  261. K[4] = P[15] * Temporary_array[1] + P[16] * Temporary_array[4] + P[17] * Temporary_array[7];
  262. K[5] = P[15] * Temporary_array[2] + P[16] * Temporary_array[5] + P[17] * Temporary_array[8];
  263. K[6] = P[24] * Temporary_array[0] + P[25] * Temporary_array[3] + P[26] * Temporary_array[6];
  264. K[7] = P[24] * Temporary_array[1] + P[25] * Temporary_array[4] + P[26] * Temporary_array[7];
  265. K[8] = P[24] * Temporary_array[2] + P[25] * Temporary_array[5] + P[26] * Temporary_array[8];
  266. K[9] = P[33] * Temporary_array[0] + P[34] * Temporary_array[3] + P[35] * Temporary_array[6];
  267. K[10] = P[33] * Temporary_array[1] + P[34] * Temporary_array[4] + P[35] * Temporary_array[7];
  268. K[11] = P[33] * Temporary_array[2] + P[34] * Temporary_array[5] + P[35] * Temporary_array[8];
  269. K[12] = P[42] * Temporary_array[0] + P[43] * Temporary_array[3] + P[44] * Temporary_array[6];
  270. K[13] = P[42] * Temporary_array[1] + P[43] * Temporary_array[4] + P[44] * Temporary_array[7];
  271. K[14] = P[42] * Temporary_array[2] + P[43] * Temporary_array[5] + P[44] * Temporary_array[8];
  272. K[15] = P[51] * Temporary_array[0] + P[52] * Temporary_array[3] + P[53] * Temporary_array[6];
  273. K[16] = P[51] * Temporary_array[1] + P[52] * Temporary_array[4] + P[53] * Temporary_array[7];
  274. K[17] = P[51] * Temporary_array[2] + P[52] * Temporary_array[5] + P[53] * Temporary_array[8];
  275. K[18] = P[60] * Temporary_array[0] + P[61] * Temporary_array[3] + P[62] * Temporary_array[6];
  276. K[19] = P[60] * Temporary_array[1] + P[61] * Temporary_array[4] + P[62] * Temporary_array[7];
  277. K[20] = P[60] * Temporary_array[2] + P[61] * Temporary_array[5] + P[62] * Temporary_array[8];
  278. K[21] = P[69] * Temporary_array[0] + P[70] * Temporary_array[3] + P[71] * Temporary_array[6];
  279. K[22] = P[69] * Temporary_array[1] + P[70] * Temporary_array[4] + P[71] * Temporary_array[7];
  280. K[23] = P[69] * Temporary_array[2] + P[70] * Temporary_array[5] + P[71] * Temporary_array[8];
  281. K[24] = P[78] * Temporary_array[0] + P[79] * Temporary_array[3] + P[80] * Temporary_array[6];
  282. K[25] = P[78] * Temporary_array[1] + P[79] * Temporary_array[4] + P[80] * Temporary_array[7];
  283. K[26] = P[78] * Temporary_array[2] + P[79] * Temporary_array[5] + P[80] * Temporary_array[8];
  284. }
  285. void multiply9x3(float *K, float *vel_n, float* delta_x)
  286. {
  287. int i = 0;
  288. for (i = 0; i < 9; i++)
  289. {
  290. delta_x[i] = K[i * 3] * vel_n[0] + K[i * 3 + 1] * vel_n[1] + K[i * 3 + 2] * vel_n[2];
  291. }
  292. }
  293. void State_covariance_matrix_corr(float *P, float *P_tmp, float *K)
  294. {
  295. int i = 0;
  296. int j = 0;
  297. for (i = 0; i < 9; i++) {
  298. for (j = 0; j < 9; j++) {
  299. P_tmp[i * 9 + j] = P[i * 9 + j] - K[3 * i] * P[54 + j] - K[3 * i + 1] * P[63 + j] - K[3 * i + 2] * P[72 + j];
  300. }
  301. }
  302. memcpy(P, P_tmp, 81 * sizeof(float));
  303. }
  304. void Att_matrix_corr(float *C, float *C_prev, float *Temporary_array, float *Temporary_array1, float *delta_x)
  305. {
  306. Temporary_array[0] = 2.0;
  307. Temporary_array[1] = -delta_x[2];
  308. Temporary_array[2] = delta_x[1];
  309. Temporary_array[3] = delta_x[2];
  310. Temporary_array[4] = 2.0;
  311. Temporary_array[5] = -delta_x[0];
  312. Temporary_array[6] = -delta_x[1];
  313. Temporary_array[7] = delta_x[0];
  314. Temporary_array[8] = 2.0;
  315. invert3x3(Temporary_array, Temporary_array1);
  316. Temporary_array[0] = 2.0;
  317. Temporary_array[1] = delta_x[2];
  318. Temporary_array[2] = -delta_x[1];
  319. Temporary_array[3] = -delta_x[2];
  320. Temporary_array[4] = 2.0;
  321. Temporary_array[5] = delta_x[0];
  322. Temporary_array[6] = delta_x[1];
  323. Temporary_array[7] = -delta_x[0];
  324. Temporary_array[8] = 2.0;
  325. multiply3x3(Temporary_array, Temporary_array1, C_prev);
  326. memcpy(Temporary_array, C_prev, 9 * sizeof(float));
  327. multiply3x3(Temporary_array, C, C_prev);
  328. memcpy(C, C_prev, 9 * sizeof(float));
  329. }
  330. void pos_n_corr(float *pos_n, float *delta_x)
  331. {
  332. pos_n[0] -= delta_x[3];
  333. pos_n[1] -= delta_x[4];
  334. pos_n[2] -= delta_x[5];
  335. }
  336. void vel_n_corr(float *vel_n, float *delta_x)
  337. {
  338. vel_n[0] -= delta_x[6];
  339. vel_n[1] -= delta_x[7];
  340. vel_n[2] -= delta_x[8];
  341. }
  342. void State_covariance_matrix_orthogonalization(float *P)
  343. {
  344. int i = 0;
  345. int j = 0;
  346. float temp;
  347. for (i = 0; i < 9; i++)
  348. for (j = i + 1; j < 9; j++)
  349. {
  350. temp = 0.5f*(P[i * 9 + j] + P[j * 9 + i]);
  351. P[i * 9 + j] = temp;
  352. P[j * 9 + i] = temp;
  353. }
  354. }
  355. void Initialize(float *gyr, float *acc)
  356. {
  357. frame_index = 1;
  358. stand_num = 0;
  359. accSize = 1.0f;
  360. accSum = 0.0f;
  361. ZUPT_STATUS = 0;
  362. HAS_RESULT = 0;
  363. memset(last_pos_n, 0, 3 * sizeof(float));
  364. memset(pos_offset, 0, 3 * sizeof(float));
  365. memset(gyr_norm_window, 0, 10 * sizeof(float));
  366. memset(gyr_z_window, 0, 10 * sizeof(float));
  367. memset(P, 0, 81 * sizeof(float));
  368. memset(acc_n, 0, 3 * sizeof(float));
  369. memset(vel_n, 0, 3 * sizeof(float));
  370. memset(pos_n, 0, 3 * sizeof(float));
  371. memset(Temporary_array1, 0, 9 * sizeof(float));
  372. memset(Temporary_array2, 0, 9 * sizeof(float));
  373. memset(K, 0, 27 * sizeof(float));
  374. memset(P_prev, 0, 81 * sizeof(float));
  375. memset(delta_x, 0, 9 * sizeof(float));
  376. memset(C, 0, 9 * sizeof(float));
  377. memset(Temporary_array1, 0, 9 * sizeof(float));
  378. memset(Temporary_array2, 0, 9 * sizeof(float));
  379. memset(press_data, 0, 10 * sizeof(int));
  380. init_attitude_matrix(C, acc);
  381. memcpy(C_prev, C, 9 * sizeof(float));
  382. }
  383. void attitude_matrix_update(float *C, float *Temporary_array1, float *Temporary_array2, float *gyr, float dt)
  384. {
  385. Temporary_array1[0] = 2.0f;
  386. Temporary_array1[1] = dt * gyr[2];
  387. Temporary_array1[2] = -dt * gyr[1];
  388. Temporary_array1[3] = -dt * gyr[2];
  389. Temporary_array1[4] = 2.0f;
  390. Temporary_array1[5] = dt * gyr[0];
  391. Temporary_array1[6] = dt * gyr[1];
  392. Temporary_array1[7] = -dt * gyr[0];
  393. Temporary_array1[8] = 2.0f;
  394. invert3x3(Temporary_array1, Temporary_array2);
  395. memset(Temporary_array1, 0, 9 * sizeof(float));
  396. Temporary_array1[0] = 2 * C[0] + C[1] * dt * gyr[2] - C[2] * dt * gyr[1];
  397. Temporary_array1[1] = 2 * C[1] - C[0] * dt * gyr[2] + C[2] * dt * gyr[0];
  398. Temporary_array1[2] = 2 * C[2] + C[0] * dt * gyr[1] - C[1] * dt * gyr[0];
  399. Temporary_array1[3] = 2 * C[3] + C[4] * dt * gyr[2] - C[5] * dt * gyr[1];
  400. Temporary_array1[4] = 2 * C[4] - C[3] * dt * gyr[2] + C[5] * dt * gyr[0];
  401. Temporary_array1[5] = 2 * C[5] + C[3] * dt * gyr[1] - C[4] * dt * gyr[0];
  402. Temporary_array1[6] = 2 * C[6] + C[7] * dt * gyr[2] - C[8] * dt * gyr[1];
  403. Temporary_array1[7] = 2 * C[7] - C[6] * dt * gyr[2] + C[8] * dt * gyr[0];
  404. Temporary_array1[8] = 2 * C[8] + C[6] * dt * gyr[1] - C[7] * dt * gyr[0];
  405. multiply3x3(Temporary_array1, Temporary_array2, C);
  406. }
  407. float max_window_val(float *window, int window_size)
  408. {
  409. float val = window[0];
  410. for (int i = 0; i < window_size; i++)
  411. {
  412. if (window[i] > val)
  413. val = window[i];
  414. }
  415. return val;
  416. }
  417. int max_window_int(int *window, int window_size)
  418. {
  419. int val = window[0];
  420. for (int i = 0; i < window_size; i++)
  421. {
  422. if (window[i] > val)
  423. val = window[i];
  424. }
  425. return val;
  426. }
  427. float min_window_val(float *window, int window_size)
  428. {
  429. float val = window[0];
  430. for (int i = 0; i < window_size; i++)
  431. {
  432. if (window[i] < val)
  433. val = window[i];
  434. }
  435. return val;
  436. }
  437. int min_window_int(int *window, int window_size)
  438. {
  439. int val = window[0];
  440. for (int i = 0; i < window_size; i++)
  441. {
  442. if (window[i] < val)
  443. val = window[i];
  444. }
  445. return val;
  446. }
  447. //press_tren 函数功能:提供走路过程中上升沿,下降沿
  448. //1 为上升 2 为 下降 0为不需要得状态
  449. int press_trend(int index ,int *window, int window_size)
  450. {
  451. int i;
  452. int max_val = window[(index - 1) % window_size];
  453. int max_index = index;
  454. int min_val = max_val;
  455. int min_index = max_index;
  456. for(i = 1; i < window_size + 1; i++)
  457. {
  458. if(max_val < window[(index - i) % window_size])
  459. {
  460. max_index = index - i + 1;
  461. max_val = window[(index - i) % window_size];
  462. }
  463. if(min_val > window[(index - i) % window_size])
  464. {
  465. min_index = index - i + 1;
  466. min_val = window[(index - i) % window_size];
  467. }
  468. }
  469. if(max_index > min_index && max_val > min_val + 50000)
  470. {
  471. return 1;
  472. }
  473. if(max_index < min_index && max_val > min_val + 50000)
  474. {
  475. return 2;
  476. }
  477. return 0;
  478. }
  479. unsigned char footPDR(int num, float *gyr, float *acc, int press, int16_t* pos_res, int16_t* att, int16_t* zupt)
  480. {
  481. unsigned char movement_e = 0;
  482. for (int i = 0; i < 3; i++)
  483. {
  484. gyr[i] *= (PI / 180);
  485. acc[i] *= g;
  486. }
  487. if(num_peak == 0)
  488. {
  489. for(int i = 0; i < 3; i++)
  490. {
  491. gyr_extreme[2 * i] = gyr[i];
  492. gyr_extreme[2 * i + 1] = gyr[i];
  493. }
  494. }
  495. for (int i = 0; i < 3; i++)
  496. {
  497. if (gyr[i] < gyr_extreme[2 * i])
  498. {
  499. gyr_extreme[2 * i] = gyr[i];
  500. }
  501. if (gyr[i] > gyr_extreme[2 * i + 1])
  502. {
  503. gyr_extreme[2 * i + 1] = gyr[i];
  504. }
  505. }
  506. accSum += sqrt(acc[0] * acc[0] + acc[1] * acc[1] + acc[2] * acc[2]);
  507. for (int i = 0; i < 3; i++)
  508. {
  509. gyr_mean[i] += gyr[i];
  510. }
  511. num_peak++;
  512. if (num_peak > 500)
  513. {
  514. if (gyr_extreme[1] - gyr_extreme[0] < 0.005f && gyr_extreme[3] - gyr_extreme[2] < 0.005f && gyr_extreme[5] - gyr_extreme[4] < 0.005f)
  515. {
  516. gyrBias[0] = gyr_mean[0] / num_peak;
  517. gyrBias[1] = gyr_mean[1] / num_peak;
  518. gyrBias[2] = gyr_mean[2] / num_peak;
  519. }
  520. num_peak = 0;
  521. accSum = 0.0f;
  522. memset(gyr_mean, 0, 3 * sizeof(float));
  523. }
  524. gyr[0] -= gyrBias[0];
  525. gyr[1] -= gyrBias[1];
  526. gyr[2] -= (gyrBias[2]);
  527. float gyr_norm_xyz = sqrt(gyr[0]*gyr[0] + gyr[1]*gyr[1] + gyr[2]*gyr[2]);
  528. //需要一个滑动窗口来判断脚步是否在地上
  529. frame_index++;
  530. //下面为惯导解算
  531. if (num == 1 || frame_index < 0)
  532. {
  533. Initialize(gyr, acc);
  534. return movement_e;
  535. }
  536. //惯导解算拆分为5次迭代,缓解高量程下由于采样率过导致惯导解算有错
  537. //插值法
  538. float gyr_temp[3];
  539. float acc_temp[3];
  540. gyr_temp[0] = gyr[0] - last_gyr[0];
  541. gyr_temp[1] = gyr[1] - last_gyr[1];
  542. gyr_temp[2] = gyr[2] - last_gyr[2];
  543. acc_temp[0] = acc[0] - last_acc[0];
  544. acc_temp[1] = acc[1] - last_acc[1];
  545. acc_temp[2] = acc[2] - last_acc[2];
  546. for(int i = 1; i < 6; i++)
  547. {
  548. last_gyr[0] += 0.2f * gyr_temp[0];
  549. last_gyr[1] += 0.2f * gyr_temp[1];
  550. last_gyr[2] += 0.2f * gyr_temp[2];
  551. last_acc[0] += 0.2f * acc_temp[0];
  552. last_acc[1] += 0.2f * acc_temp[1];
  553. last_acc[2] += 0.2f * acc_temp[2];
  554. attitude_matrix_update(C, Temporary_array1, Temporary_array2, last_gyr, 0.2f*dt);
  555. multiply3x1(C, last_acc, acc_n);
  556. vel_n[0] = last_vel_n[0] + acc_n[0] * dt * 0.2f;
  557. vel_n[1] = last_vel_n[1] + acc_n[1] * dt * 0.2f;
  558. vel_n[2] = last_vel_n[2] + (acc_n[2] - g) * dt * 0.2f;
  559. pos_n[0] = pos_n[0] + (vel_n[0] + last_vel_n[0]) * dt * 0.1f;
  560. pos_n[1] = pos_n[1] + (vel_n[1] + last_vel_n[1]) * dt * 0.1f;
  561. pos_n[2] = pos_n[2] + (vel_n[2] + last_vel_n[2]) * dt * 0.1f;
  562. memcpy(last_vel_n, vel_n, 3 * sizeof(float));
  563. }
  564. // attitude_matrix_update(C, Temporary_array1, Temporary_array2, gyr, dt);
  565. // multiply3x1(C, acc, acc_n);
  566. // vel_n[0] = vel_n[0] + acc_n[0] * dt;
  567. // vel_n[1] = vel_n[1] + acc_n[1] * dt;
  568. // vel_n[2] = vel_n[2] + (acc_n[2] - g) * dt;
  569. // pos_n[0] = pos_n[0] + vel_n[0] * dt;
  570. // pos_n[1] = pos_n[1] + vel_n[1] * dt;
  571. // pos_n[2] = pos_n[2] + vel_n[2] * dt;
  572. memcpy(last_gyr, gyr, 3 * sizeof(float));
  573. memcpy(last_acc, acc, 3 * sizeof(float));
  574. //P = F*P*F' + Q;
  575. State_covariance_matrix_update(P, acc_n);
  576. int window_index = (frame_index - 1) % 10;
  577. float gyr_norm_xz = sqrt(gyr[0] * gyr[0] + gyr[1] * gyr[1] +gyr[2] * gyr[2]);
  578. gyr_norm_window[window_index] = gyr_norm_xz;
  579. press_data[window_index] = press;
  580. //当press_trend函数返回是1,判断为踩地上
  581. // 返回2 的时候,判断为离地
  582. // 返回0 的时候,需要保持状态
  583. int press_trend_val = press_trend(frame_index, press_data, 10);
  584. if(press_trend_val == 1)
  585. {
  586. ZUPT_STATUS = 1;
  587. }
  588. else if(press_trend_val == 2)
  589. {
  590. ZUPT_STATUS = 2;
  591. }
  592. //RUN_ZUPT mean detect on floor when running
  593. int RUN_ZUPT = 0;
  594. if((frame_index > 10 && ZUPT_STATUS == 1))
  595. RUN_ZUPT = 1;
  596. //STAND_ZUPT mean detect on floor when no any moving
  597. int STAND_ZUPT = 0;
  598. if((frame_index > 15 && gyr_norm_window[window_index] < 0.35f && fabs(min_window_val(gyr_norm_window, 10) - max_window_val(gyr_norm_window, 10)) < 0.1f))
  599. STAND_ZUPT = 1;
  600. //zupt
  601. if ((STAND_ZUPT || RUN_ZUPT))
  602. {
  603. //计算一步的距离及步数+1
  604. //做一个计算时间时间,防止非正常走路的情况出来
  605. if(last_stage == 2 && frame_index - step_time > 30)
  606. {
  607. cal_step_data();
  608. step_time = frame_index;
  609. }
  610. stand_num = stand_num + 1;
  611. //K = P*H'/(H*P*H' + R);
  612. Kalfman_gain(P, Temporary_array1, Temporary_array2, K);
  613. //delta_x = K * [vel_n(:,i);];
  614. multiply9x3(K, vel_n, delta_x);
  615. State_covariance_matrix_corr(P, P_prev, K);
  616. //这里先从设置 delta_x(3) = atan2(C(2,1),C(1,1));
  617. //意味着每一步的参考方向是IMU X轴方向
  618. // delta_x[2] = atan2(C[3], C[0]);
  619. //// theta = -1.7801 + atan2(C[3], C[0]);
  620. // theta = 0.0f;
  621. Att_matrix_corr(C, C_prev, Temporary_array1, Temporary_array2, delta_x);
  622. pos_n_corr(pos_n, delta_x);
  623. vel_n_corr(vel_n, delta_x);
  624. memset(vel_n, 0, 3*sizeof(float));
  625. last_pos_n[0] = pos_n[0];
  626. last_pos_n[1] = pos_n[1];
  627. last_pos_n[2] = pos_n[2];
  628. HAS_RESULT = 1;
  629. last_stage = 1;
  630. pos_offset_min = 0;
  631. *zupt = 1;
  632. }
  633. else
  634. {
  635. stand_num = 0;
  636. *zupt = 0;
  637. }
  638. /*
  639. if(RUN_ZUPT)
  640. {
  641. //计算一步的距离及步数+1
  642. //做一个计算时间时间,防止非正常走路的情况出来
  643. if(last_stage == 2 && frame_index - step_time > 30)
  644. {
  645. cal_step_data();
  646. step_time = frame_index;
  647. }
  648. }
  649. */
  650. State_covariance_matrix_orthogonalization(P);
  651. memcpy(last_vel_n, vel_n, 3 * sizeof(float));
  652. /*theta = -0.61;
  653. temp = [cos(theta), sin(theta); -sin(theta), cos(theta)];
  654. pos_temp = temp * [pos_result(1); pos_result(2)];
  655. pos_result(1) = pos_temp(1);
  656. pos_result(2) = pos_temp(2);*/
  657. float pos_offset_temp0 = pos_n[0] - last_pos_n[0];
  658. float pos_offset_temp1 = pos_n[1] - last_pos_n[1];
  659. float pos_offset_temp2 = pos_n[2] - last_pos_n[2];
  660. // pos_offset[0] = cos(theta) * pos_offset_temp0 + sin(theta) * pos_offset_temp1;
  661. // pos_offset[1] = -sin(theta) * pos_offset_temp0+ cos(theta) * pos_offset_temp1;
  662. pos_offset[0] = pos_offset_temp0;
  663. pos_offset[1] = pos_offset_temp1;
  664. pos_offset[2] = pos_offset_temp2;
  665. // memcpy(pos_res, acc_n, 3 * sizeof(float));
  666. //修改位置判断
  667. // if(HAS_RESULT == 1 && pos_offset[0] > 0.2f)
  668. // {
  669. // movement_e = 1;
  670. // HAS_RESULT = 0;
  671. // }
  672. att[0] = (short) (atan2(C[3], C[0]) * 10000.f); //yaw
  673. att[1] = (short) (asin(-C[6]) * 10000.f); //pitch
  674. att[2] = (short) (atan2(C[7], C[8]) * 10000.f); //roll
  675. pos_res[0] = (short) (pos_offset[0] * 100.0f);
  676. pos_res[1] = (short) (pos_offset[1] * 100.0f);
  677. pos_res[2] = (short) (pos_offset[2] * 100.0f);
  678. // if(pos_res[0] < pos_offset_min)
  679. // {
  680. // pos_offset_min = pos_res[0];
  681. // }
  682. //
  683. //
  684. // if(HAS_RESULT == 1 && pos_res[0] - pos_offset_min > 20)
  685. // {
  686. // movement_e = 1;
  687. // HAS_RESULT = 0;
  688. // }
  689. if((pos_res[0] * pos_res[0] + pos_res[1] * pos_res[1] + pos_res[2] * pos_res[2]) > 10)
  690. {
  691. last_stage = 2;
  692. }
  693. return movement_e;
  694. }