2 * This file is part of Cleanflight and Betaflight.
4 * Cleanflight and Betaflight are free software. You can redistribute
5 * this software and/or modify this software under the terms of the
6 * GNU General Public License as published by the Free Software
7 * Foundation, either version 3 of the License, or (at your option)
10 * Cleanflight and Betaflight are distributed in the hope that they
11 * will be useful, but WITHOUT ANY WARRANTY; without even the implied
12 * warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 * See the GNU General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this software.
18 * If not, see <http://www.gnu.org/licenses/>.
21 /* original work by Rav
23 * 2018_07 updated by ctzsnooze to post filter, wider Q, different peak detection
24 * coding assistance and advice from DieHertz, Rav, eTracer
25 * test pilots icr4sh, UAV Tech, Flint723
27 * 2021_02 updated by KarateBrot: switched FFT with SDFT, multiple notches per axis
28 * test pilots: Sugar K, bizmar
36 #ifdef USE_GYRO_DATA_ANALYSE
37 #include "build/debug.h"
39 #include "common/filter.h"
40 #include "common/maths.h"
41 #include "common/sdft.h"
42 #include "common/time.h"
43 #include "common/utils.h"
45 #include "drivers/accgyro/accgyro.h"
46 #include "drivers/time.h"
48 #include "sensors/gyro.h"
52 #include "gyroanalyse.h"
54 // SDFT_SAMPLE_SIZE defaults to 72 (common/sdft.h).
55 // We get 36 frequency bins from 72 consecutive data values, called SDFT_BIN_COUNT (common/sdft.h)
56 // Bin 0 is DC and can't be used.
57 // Only bins 1 to 35 are usable.
59 // A gyro sample is collected every PID loop.
60 // maxSampleCount recent gyro values are accumulated and averaged
61 // to ensure that 72 samples are collected at the right rate for the required SDFT bandwidth.
63 // For an 8k PID loop, at default 600hz max, 6 sequential gyro data points are averaged, SDFT runs 1333Hz.
64 // Upper limit of SDFT is half that frequency, eg 666Hz by default.
65 // At 8k, if user sets a max of 300Hz, int(8000/600) = 13, sdftSampleRateHz = 615Hz, range 307Hz.
66 // Note that lower max requires more samples to be averaged, increasing precision but taking longer to get enough samples.
67 // For Bosch at 3200Hz gyro, max of 600, int(3200/1200) = 2, sdftSampleRateHz = 1600, range to 800hz.
68 // For Bosch on XClass, better to set a max of 300, int(3200/600) = 5, sdftSampleRateHz = 640, range to 320Hz.
70 // When sampleCount reaches maxSampleCount, the averaged gyro value is put into the corresponding SDFT.
71 // At 8k, with 600Hz max, maxSampleCount = 6, this happens every 6 * 0.125us, or every 0.75ms.
72 // Hence to completely replace all 72 samples of the SDFT input buffer with clean new data takes 54ms.
74 // The SDFT code is split into steps. It takes 4 PID loops to calculate the SDFT, track peaks and update the filters for one axis.
75 // Since there are three axes, it takes 12 PID loops to completely update all axes.
76 // At 8k, any one axis gets updated at 8000 / 12 or 666hz or every 1.5ms
77 // In this time, 2 points in the SDFT buffer will have changed.
78 // At 4k, it takes twice as long to update an axis, i.e. each axis updates only every 3ms.
79 // Four points in the buffer will have changed in that time, and each point will be the average of three samples.
80 // Hence output jitter at 4k is about four times worse than at 8k. At 2k output jitter is quite bad.
82 // Each SDFT output bin has width sdftSampleRateHz/72, ie 18.5Hz per bin at 1333Hz.
83 // Usable bandwidth is half this, ie 666Hz if sdftSampleRateHz is 1333Hz, i.e. bin 1 is 18.5Hz, bin 2 is 37.0Hz etc.
85 #define DYN_NOTCH_SMOOTH_HZ 4
86 #define DYN_NOTCH_CALC_TICKS (XYZ_AXIS_COUNT * STEP_COUNT) // 3 axes and 4 steps per axis
87 #define DYN_NOTCH_OSD_MIN_THROTTLE 20
93 STEP_CALC_FREQUENCIES
,
99 typedef struct peak_s
{
106 static sdft_t FAST_DATA_ZERO_INIT sdft
[XYZ_AXIS_COUNT
];
107 static peak_t FAST_DATA_ZERO_INIT peaks
[DYN_NOTCH_COUNT_MAX
];
108 static float FAST_DATA_ZERO_INIT sdftData
[SDFT_BIN_COUNT
];
109 static uint16_t FAST_DATA_ZERO_INIT sdftSampleRateHz
;
110 static float FAST_DATA_ZERO_INIT sdftResolutionHz
;
111 static uint8_t FAST_DATA_ZERO_INIT sdftStartBin
;
112 static uint8_t FAST_DATA_ZERO_INIT sdftEndBin
;
113 static float FAST_DATA_ZERO_INIT sdftMeanSq
;
114 static uint16_t FAST_DATA_ZERO_INIT dynNotchBandwidthHz
;
115 static uint16_t FAST_DATA_ZERO_INIT dynNotchMinHz
;
116 static uint16_t FAST_DATA_ZERO_INIT dynNotchMaxHz
;
117 static uint16_t FAST_DATA_ZERO_INIT dynNotchMaxFFT
;
118 static float FAST_DATA_ZERO_INIT smoothFactor
;
119 static uint8_t FAST_DATA_ZERO_INIT numSamples
;
121 void gyroDataAnalyseInit(gyroAnalyseState_t
*state
, uint32_t targetLooptimeUs
)
123 // initialise even if FEATURE_DYNAMIC_FILTER not set, since it may be set later
124 dynNotchBandwidthHz
= gyroConfig()->dyn_notch_bandwidth_hz
;
125 dynNotchMinHz
= gyroConfig()->dyn_notch_min_hz
;
126 dynNotchMaxHz
= MAX(2 * dynNotchMinHz
, gyroConfig()->dyn_notch_max_hz
);
128 // gyroDataAnalyse() is running at targetLoopRateHz (which is PID loop rate aka. 1e6f/gyro.targetLooptimeUs)
129 const int32_t targetLoopRateHz
= lrintf((1.0f
/ targetLooptimeUs
) * 1e6f
);
130 numSamples
= MAX(1, targetLoopRateHz
/ (2 * dynNotchMaxHz
)); // 600hz, 8k looptime, 6.00
132 sdftSampleRateHz
= targetLoopRateHz
/ numSamples
;
133 // eg 8k, user max 600hz, int(8000/1200) = 6 (6.666), sdftSampleRateHz = 1333hz, range 666Hz
134 // eg 4k, user max 600hz, int(4000/1200) = 3 (3.333), sdftSampleRateHz = 1333hz, range 666Hz
135 // eg 2k, user max 600hz, int(2000/1200) = 1 (1.666) sdftSampleRateHz = 2000hz, range 1000Hz
136 // eg 2k, user max 400hz, int(2000/800) = 2 (2.5) sdftSampleRateHz = 1000hz, range 500Hz
137 // eg 1k, user max 600hz, int(1000/1200) = 1 (max(1,0.8333)) sdftSampleRateHz = 1000hz, range 500Hz
138 // the upper limit of DN is always going to be the Nyquist frequency (= sampleRate / 2)
140 sdftResolutionHz
= (float)sdftSampleRateHz
/ SDFT_SAMPLE_SIZE
; // 13.3hz per bin at 8k
141 sdftStartBin
= MAX(2, lrintf(dynNotchMinHz
/ sdftResolutionHz
+ 0.5f
)); // can't use bin 0 because it is DC.
142 sdftEndBin
= MIN(SDFT_BIN_COUNT
- 1, lrintf(dynNotchMaxHz
/ sdftResolutionHz
+ 0.5f
)); // can't use more than SDFT_BIN_COUNT bins.
143 smoothFactor
= pt1FilterGain(DYN_NOTCH_SMOOTH_HZ
, DYN_NOTCH_CALC_TICKS
/ (float)targetLoopRateHz
); // minimum PT1 k value
145 for (uint8_t axis
= 0; axis
< XYZ_AXIS_COUNT
; axis
++) {
146 sdftInit(&sdft
[axis
], sdftStartBin
, sdftEndBin
, numSamples
);
149 state
->maxSampleCount
= numSamples
;
150 state
->maxSampleCountRcp
= 1.0f
/ state
->maxSampleCount
;
152 for (uint8_t axis
= 0; axis
< XYZ_AXIS_COUNT
; axis
++) {
153 for (uint8_t p
= 0; p
< gyro
.notchFilterDynCount
; p
++) {
154 // any init value is fine, but evenly spreading centerFreqs across frequency range makes notch filters stick to peaks quicker
155 state
->centerFreq
[axis
][p
] = (p
+ 0.5f
) * (dynNotchMaxHz
- dynNotchMinHz
) / (float)gyro
.notchFilterDynCount
+ dynNotchMinHz
;
160 // Collect gyro data, to be downsampled and analysed in gyroDataAnalyse() function
161 void gyroDataAnalysePush(gyroAnalyseState_t
*state
, const uint8_t axis
, const float sample
)
163 state
->oversampledGyroAccumulator
[axis
] += sample
;
166 static void gyroDataAnalyseUpdate(gyroAnalyseState_t
*state
);
168 // Downsample and analyse gyro data
169 FAST_CODE
void gyroDataAnalyse(gyroAnalyseState_t
*state
)
171 // samples should have been pushed by `gyroDataAnalysePush`
172 // if gyro sampling is > 1kHz, accumulate and average multiple gyro samples
173 if (state
->sampleCount
== state
->maxSampleCount
) {
174 state
->sampleCount
= 0;
176 // calculate mean value of accumulated samples
177 for (uint8_t axis
= 0; axis
< XYZ_AXIS_COUNT
; axis
++) {
178 const float sample
= state
->oversampledGyroAccumulator
[axis
] * state
->maxSampleCountRcp
;
179 state
->downsampledGyroData
[axis
] = sample
;
181 DEBUG_SET(DEBUG_FFT
, 2, lrintf(sample
));
183 state
->oversampledGyroAccumulator
[axis
] = 0;
186 // We need DYN_NOTCH_CALC_TICKS ticks to update all axes with newly sampled value
187 // recalculation of filters takes 4 calls per axis => each filter gets updated every DYN_NOTCH_CALC_TICKS calls
188 // at 8kHz PID loop rate this means 8kHz / 4 / 3 = 666Hz => update every 1.5ms
189 // at 4kHz PID loop rate this means 4kHz / 4 / 3 = 333Hz => update every 3ms
190 state
->updateTicks
= DYN_NOTCH_CALC_TICKS
;
194 // SDFT processing in batches to synchronize with incoming downsampled data
195 for (uint8_t axis
= 0; axis
< XYZ_AXIS_COUNT
; axis
++) {
196 sdftPushBatch(&sdft
[axis
], &state
->downsampledGyroData
[axis
], &state
->sampleCount
);
198 state
->sampleCount
++;
200 // Find frequency peaks and update filters
201 if (state
->updateTicks
> 0) {
202 gyroDataAnalyseUpdate(state
);
203 --state
->updateTicks
;
207 // Find frequency peaks and update filters
208 static FAST_CODE_NOINLINE
void gyroDataAnalyseUpdate(gyroAnalyseState_t
*state
)
210 uint32_t startTime
= 0;
211 if (debugMode
== (DEBUG_FFT_TIME
)) {
212 startTime
= micros();
215 DEBUG_SET(DEBUG_FFT_TIME
, 0, state
->updateStep
);
217 switch (state
->updateStep
) {
219 case STEP_WINDOW
: // 6us @ F722
221 sdftWinSq(&sdft
[state
->updateAxis
], sdftData
);
223 // Calculate mean square over frequency range (= average power of vibrations)
225 for (uint8_t bin
= (sdftStartBin
+ 1); bin
< sdftEndBin
; bin
++) { // don't use startBin or endBin because they are not windowed properly
226 sdftMeanSq
+= sdftData
[bin
]; // sdftData is already squared (see sdftWinSq)
228 sdftMeanSq
/= sdftEndBin
- sdftStartBin
- 1;
230 DEBUG_SET(DEBUG_FFT_TIME
, 1, micros() - startTime
);
234 case STEP_DETECT_PEAKS
: // 6us @ F722
236 // Get memory ready for new peak data on current axis
237 for (uint8_t p
= 0; p
< gyro
.notchFilterDynCount
; p
++) {
239 peaks
[p
].value
= 0.0f
;
242 // Search for N biggest peaks in frequency spectrum
243 for (uint8_t bin
= (sdftStartBin
+ 1); bin
< sdftEndBin
; bin
++) {
244 // Check if bin is peak
245 if ((sdftData
[bin
] > sdftData
[bin
- 1]) && (sdftData
[bin
] > sdftData
[bin
+ 1])) {
246 // Check if peak is big enough to be one of N biggest peaks.
247 // If so, insert peak and sort peaks in descending height order
248 for (uint8_t p
= 0; p
< gyro
.notchFilterDynCount
; p
++) {
249 if (sdftData
[bin
] > peaks
[p
].value
) {
250 for (uint8_t k
= gyro
.notchFilterDynCount
- 1; k
> p
; k
--) {
251 peaks
[k
] = peaks
[k
- 1];
254 peaks
[p
].value
= sdftData
[bin
];
258 bin
++; // If bin is peak, next bin can't be peak => jump it
262 // Sort N biggest peaks in ascending bin order (example: 3, 8, 25, 0, 0, ..., 0)
263 for (uint8_t p
= gyro
.notchFilterDynCount
- 1; p
> 0; p
--) {
264 for (uint8_t k
= 0; k
< p
; k
++) {
265 // Swap peaks but ignore swapping void peaks (bin = 0). This leaves
266 // void peaks at the end of peaks array without moving them
267 if (peaks
[k
].bin
> peaks
[k
+ 1].bin
&& peaks
[k
+ 1].bin
!= 0) {
268 peak_t temp
= peaks
[k
];
269 peaks
[k
] = peaks
[k
+ 1];
275 DEBUG_SET(DEBUG_FFT_TIME
, 1, micros() - startTime
);
279 case STEP_CALC_FREQUENCIES
: // 4us @ F722
281 for (uint8_t p
= 0; p
< gyro
.notchFilterDynCount
; p
++) {
283 // Only update state->centerFreq if there is a peak (ignore void peaks) and if peak is above noise floor
284 if (peaks
[p
].bin
!= 0 && peaks
[p
].value
> sdftMeanSq
) {
286 // accumulate sdftSum and sdftWeightedSum from peak bin, and shoulder bins either side of peak
287 float squaredData
= peaks
[p
].value
; // peak data already squared (see sdftWinSq)
288 float sdftSum
= squaredData
;
289 float sdftWeightedSum
= squaredData
* peaks
[p
].bin
;
291 // accumulate upper shoulder unless it would be sdftEndBin
292 uint8_t shoulderBin
= peaks
[p
].bin
+ 1;
293 if (shoulderBin
< sdftEndBin
) {
294 squaredData
= sdftData
[shoulderBin
]; // sdftData already squared (see sdftWinSq)
295 sdftSum
+= squaredData
;
296 sdftWeightedSum
+= squaredData
* shoulderBin
;
299 // accumulate lower shoulder unless lower shoulder would be bin 0 (DC)
300 if (peaks
[p
].bin
> 1) {
301 shoulderBin
= peaks
[p
].bin
- 1;
302 squaredData
= sdftData
[shoulderBin
]; // sdftData already squared (see sdftWinSq)
303 sdftSum
+= squaredData
;
304 sdftWeightedSum
+= squaredData
* shoulderBin
;
307 // get centerFreq in Hz from weighted bins
308 float centerFreq
= dynNotchMaxHz
;
309 float sdftMeanBin
= 0;
312 sdftMeanBin
= (sdftWeightedSum
/ sdftSum
);
313 centerFreq
= sdftMeanBin
* sdftResolutionHz
;
314 centerFreq
= constrainf(centerFreq
, dynNotchMinHz
, dynNotchMaxHz
);
315 // In theory, the index points to the centre frequency of the bin.
316 // at 1333hz, bin widths are 13.3Hz, so bin 2 (26.7Hz) has the range 20Hz to 33.3Hz
317 // Rav feels that maybe centerFreq = (sdftMeanBin + 0.5) * sdftResolutionHz is better
318 // empirical checking shows that not adding 0.5 works better
320 // PT1 style dynamic smoothing moves rapidly towards big peaks and slowly away, up to 8x faster
321 // DYN_NOTCH_SMOOTH_HZ = 4 & dynamicFactor = 1 .. 8 => PT1 -3dB cutoff frequency = 4Hz .. 41Hz
322 const float dynamicFactor
= constrainf(peaks
[p
].value
/ sdftMeanSq
, 1.0f
, 8.0f
);
323 state
->centerFreq
[state
->updateAxis
][p
] += smoothFactor
* dynamicFactor
* (centerFreq
- state
->centerFreq
[state
->updateAxis
][p
]);
328 if(calculateThrottlePercentAbs() > DYN_NOTCH_OSD_MIN_THROTTLE
) {
329 for (uint8_t p
= 0; p
< gyro
.notchFilterDynCount
; p
++) {
330 dynNotchMaxFFT
= MAX(dynNotchMaxFFT
, state
->centerFreq
[state
->updateAxis
][p
]);
334 if (state
->updateAxis
== gyro
.gyroDebugAxis
) {
335 for (uint8_t p
= 0; p
< gyro
.notchFilterDynCount
&& p
< 3; p
++) {
336 DEBUG_SET(DEBUG_FFT_FREQ
, p
, lrintf(state
->centerFreq
[state
->updateAxis
][p
]));
338 DEBUG_SET(DEBUG_DYN_LPF
, 1, lrintf(state
->centerFreq
[state
->updateAxis
][0]));
341 DEBUG_SET(DEBUG_FFT_TIME
, 1, micros() - startTime
);
345 case STEP_UPDATE_FILTERS
: // 7us @ F722
347 for (uint8_t p
= 0; p
< gyro
.notchFilterDynCount
; p
++) {
348 // Only update notch filter coefficients if the corresponding peak got its center frequency updated in the previous step
349 if (peaks
[p
].bin
!= 0 && peaks
[p
].value
> sdftMeanSq
) {
350 // Choose notch Q in such a way that notch bandwidth stays constant (improves prop wash handling)
351 float dynamicQ
= state
->centerFreq
[state
->updateAxis
][p
] / (float)dynNotchBandwidthHz
;
352 dynamicQ
= constrainf(dynamicQ
, 2.0f
, 10.0f
);
353 biquadFilterUpdate(&gyro
.notchFilterDyn
[state
->updateAxis
][p
], state
->centerFreq
[state
->updateAxis
][p
], gyro
.targetLooptime
, dynamicQ
, FILTER_NOTCH
);
357 DEBUG_SET(DEBUG_FFT_TIME
, 1, micros() - startTime
);
359 state
->updateAxis
= (state
->updateAxis
+ 1) % XYZ_AXIS_COUNT
;
363 state
->updateStep
= (state
->updateStep
+ 1) % STEP_COUNT
;
367 uint16_t getMaxFFT(void) {
368 return dynNotchMaxFFT
;
371 void resetMaxFFT(void) {
375 #endif // USE_GYRO_DATA_ANALYSE