math_utils.c
8.7 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
/*
* 常用数学函数实现
*/
#include "math_utils.h"
#include "isoCalculate.h"
// normalizing radian angle to [0, 2pi]
float _normalizeAngle(float angle)
{
float a = fmod(angle, _2PI); // fmod()函数用于浮点数的取余运算
return a >= 0.0f ? a : (a + _2PI);
}
//======================================== LowPassFilter SimpleVer ==============================
// 低通滤波
void LowPassFilter(float *value)
{
static float valueLast = 0.0f;
*value = 0.9f * valueLast + 0.1f * (*value);
valueLast = *value;
}
// 低通滤波
void LowPassFilter1(float *value)
{
static float valueLast = 0.0f;
*value = 0.9f * valueLast + 0.1f * (*value);
valueLast = *value;
}
//======================================== LowPassFilter STD Ver ==============================
void LowPassFilterSTD_Init(LowPassFSTD_s *config, float Tf)
{
config->Tf = Tf;
TimeFlash(&(config->timeS));
}
float LowPassFilterSTD_Process(LowPassFSTD_s *config, float x)
{
float timePass = 0;
float Tf = config->Tf;
float y_prev = config->y_prev;
timePass = TimeFlash(&(config->timeS));
// unsigned long timestamp = _micros();
// if (dt < 0.0f ) dt = 1e-3f;
// else if(dt > 0.3f) {
// y_prev = x;
// timestamp_prev = timestamp;
// return x;
// }
float alpha = Tf / (Tf + timePass);
y_prev = alpha * y_prev + (1.0f - alpha) * x;
config->y_prev = y_prev;
return y_prev;
}
float LowPassFilterSTD_Process_EX(LowPassFSTD_s *config, float x, float y_prev)
{
float timePass = 0;
float Tf = config->Tf;
// float y_prev = config->y_prev;
timePass = TimeFlash(&(config->timeS));
// unsigned long timestamp = _micros();
// if (dt < 0.0f ) dt = 1e-3f;
// else if(dt > 0.3f) {
// y_prev = x;
// timestamp_prev = timestamp;
// return x;
// }
float alpha = Tf / (Tf + timePass);
y_prev = alpha * y_prev + (1.0f - alpha) * x;
config->y_prev = y_prev;
return y_prev;
}
//======================================== LowPassFilter Ring Ver ==============================
void LowPassFilterSTDRing_Init(LowPassFSTDRing_s *config, float Tf, float ringRange)
{
config->Tf = Tf;
config->ringSet = fabsf(ringRange);
TimeFlash(&(config->timeS));
}
float LowPassFilterSTDRing_Process_EX(LowPassFSTDRing_s *config, float x, float y_prev)
{
float timePass = 0;
float Tf = config->Tf;
// float y_prev = config->y_prev;
timePass = TimeFlash(&(config->timeS));
float err = shortARCCalculate(y_prev, x);
float x_New = y_prev + err;
float alpha = Tf / (Tf + timePass);
y_prev = alpha * y_prev + (1.0f - alpha) * x_New;
y_prev = keepIn180UpDown(y_prev);
// config->y_prev = y_prev;
return y_prev;
}
//======================================== Interval Filter Utils ==============================
#define INTERFILTER_OFFSET 3
int32_t interF_PreData = 0;
uint16_t intervalFilter(uint16_t input)
{
if ((input > interF_PreData) && (input < (uint16_t)(interF_PreData + INTERFILTER_OFFSET))) // Overflow UP
{
interF_PreData = input;
}
else
{
goto FLASH_INTERF_PREDATA;
}
if ((input < interF_PreData) && (input > (uint16_t)(interF_PreData + INTERFILTER_OFFSET)))
{
interF_PreData = input;
}
else
{
goto FLASH_INTERF_PREDATA;
}
return interF_PreData;
FLASH_INTERF_PREDATA:
interF_PreData = input;
return interF_PreData;
}
uint16_t Inter1_Rec = 0; // 区间滤波的中值
uint16_t Inter1_Range = 2;
uint16_t inter_filter(uint16_t x)
{
if (x > (Inter1_Rec + Inter1_Range))
{
Inter1_Rec = x - Inter1_Range; // 把区间顶上去
}
else if (x < (Inter1_Rec - Inter1_Range))
{
Inter1_Rec = x + Inter1_Range;
}
return Inter1_Rec;
}
uint16_t Inter2_Rec = 0; // 区间滤波的中值
uint16_t Inter2_Range = 3;
uint16_t inter_filterC1(uint16_t x)
{
if (x > (Inter2_Rec + Inter2_Range))
{
Inter2_Rec = x - Inter2_Range; // 把区间顶上去
}
else if (x < (Inter2_Rec - Inter2_Range))
{
Inter2_Rec = x + Inter2_Range;
}
return Inter2_Rec;
}
uint16_t inter_filter_GetRange(){
return Inter1_Range;
}
int inter_filter_SetRange(uint16_t range){
Inter1_Range = range;
return 0;
}
uint16_t inter_filterC1_GetRange(){
return Inter2_Range;
}
int inter_filterC1_SetRange(uint16_t range){
Inter2_Range = range;
return 0;
}
//======================================== Other Utils ==============================
float calculateArcMidpointAngle(float angle1, float angle2)
{
float a, b, arc1, arc2, c;
if (angle1 > angle2)
{
a = angle1;
b = angle2;
}
else
{
b = angle1;
a = angle2;
}
arc1 = _2PI - a + b;
arc2 = a - b;
if (arc1 < arc2)
{
c = a + arc1 * 0.5f;
if (c > _2PI)
{
return c - _2PI;
}
else
{
return c;
}
}
else
{
return arc2 * 0.5f + b;
}
}
//======================================== KalmanFilter Utils ==============================
int kalmanFilter_Init(KFP *kfp, KFP_CONFIG config)
{
memset(kfp, 0x00, sizeof(KFP));
kfp->config = config;
return 0;
}
/**
*卡尔曼滤波器
*@param KFP *kfp 卡尔曼结构体参数
* float input 需要滤波的参数的测量值(即传感器的采集值)
*@return 滤波后的参数(最优值)
*/
float kalmanFilter(KFP *kfp, float input)
{
// 预测协方差方程:k时刻系统估算协方差 = k-1时刻的系统协方差 + 过程噪声协方差
kfp->Now_P = kfp->LastP + kfp->config.Q;
// 卡尔曼增益方程:卡尔曼增益 = k时刻系统估算协方差 / (k时刻系统估算协方差 + 观测噪声协方差)
kfp->Kg = kfp->Now_P / (kfp->Now_P + kfp->config.R);
// 更新最优值方程:k时刻状态变量的最优值 = 状态变量的预测值 + 卡尔曼增益 * (测量值 - 状态变量的预测值)
kfp->out = kfp->out + kfp->Kg * (input - kfp->out); // 因为这一次的预测值就是上一次的输出值
// 更新协方差方程: 本次的系统协方差付给 kfp->LastP 威下一次运算准备。
kfp->LastP = (1 - kfp->Kg) * kfp->Now_P;
return kfp->out;
}
//======================================== KalmanFilter Utils (old) ==============================
// static float p_last = 0;
// static float x_last = 0;
// //过程噪音 买的破温度计有多破。。
// #define P_Q 0.1
// static float kalman_P_Q = 0.002f;
// //测量噪声
// #define M_R 0.01
// static float kalman_M_R = 0.6f;
// /*
// Q:过程噪声,Q增大,动态响应变快,收敛稳定性变坏
// R:测量噪声,R增大,动态响应变慢,收敛稳定性变好
// 其中p的初值可以随便取,但是不能为0(为0的话卡尔曼滤波器就认为已经是最优滤波器了)
// q,r的值需要我们试出来,讲白了就是(买的破温度计有多破,以及你的超人力有多强)
// r参数调整滤波后的曲线与实测曲线的相近程度,r越小越接近。
// q参数调滤波后的曲线平滑程度,q越小越平滑。
// */
// // 2. 以高度为例 定义卡尔曼结构体并初始化参数
// // KFP KFP_height = {0.02, 0, 0, 0, 0.001, 0.543};
// void kalmanFilter_Set_R(float input){
// kalman_M_R = input;
// }
// void kalmanFilter_Set_Q(float input){
// kalman_P_Q = input;
// }
// float KalmanFilter_A(float ResrcData)
// {
// float R = kalman_M_R;
// float Q = kalman_P_Q;
// float x_mid = x_last;
// float x_now;
// float p_mid ;
// float p_now;
// float kg;
// //这里p_last 等于 kalmanFilter_A 的p直接取0
// x_mid=x_last; //x_last=x(k-1|k-1),x_mid=x(k|k-1)
// p_mid=p_last+Q; //p_mid=p(k|k-1),p_last=p(k-1|k-1),Q=噪声
// /*
// * 卡尔曼滤波的五个重要公式
// */
// kg=p_mid/(p_mid+R); //kg为kalman filter,R 为噪声
// x_now=x_mid+kg*(ResrcData-x_mid); //估计出的最优值
// p_now=(1-kg)*p_mid; //最优值对应的covariance
// p_last = p_now; //更新covariance 值
// x_last = x_now; //更新系统状态值
// return x_now;
// }
// float kalmanFilter_B(float inData)
// {
// static float prevData=0;
// //其中p的初值可以随便取,但是不能为0(为0的话卡尔曼滤波器就认为已经是最优滤波器了)
// static float p=0.01, q=P_Q, r=M_R, kGain=0;
// p = p+q;
// kGain = p/(p+r);
// inData = prevData+(kGain*(inData-prevData));
// p = (1-kGain)*p;
// prevData = inData;
// return inData;
// }