feat: Add Vibration Percentile Analysis

[nerdCopter]

Priority: Lower

Why Useful:

• Shows statistical distribution of vibrations (90th, 95th, 99th percentile)
• Better than just peak values
• Identifies consistent vs transient noise

Implementation:
Cumulative distribution plots for gyro noise energy

Impact: Advanced statistical analysis; useful for understanding noise consistency.

[coderabbitai]

📝 CodeRabbit Plan Mode

Generate an implementation plan and prompts that you can use with your favorite coding agent.

✅ Successfully generated implementation plan.

💡 Use @coderabbitai <feedback> to replan

Examples

• Example 1
• Example 2

---

🔗 Similar Issues

Related Issues

• feat: Add Coherence Plot (Input-Output Correlation) #106
• feat: Add Motor Command Symmetry Analysis #107
• feat: Add Phase Margin Analysis (Bode Plot) #104
• feat: Add Dynamic Notch Filter Activity Over Time #109
• feat: Add Setpoint Derivative (Rate of Change) Plot #103

👤 Suggested Assignees

• nerdCopter

---

🧪 Issue enrichment is currently in open beta.

You can configure auto-planning by selecting labels in the issue_enrichment configuration.

To disable automatic issue enrichment, add the following to your .coderabbit.yaml:

issue_enrichment:
  auto_enrich:
    enabled: false

💬 Have feedback or questions? Drop into our discord!

[nerdCopter]

@coderabbitai, we recently found that bode-plots were essentially useless for real flights. confirm or disprove that "Vibration Percentile Analysis" is useful or not for real Betaflight flight logs. deep research and fact-check all findings.

[coderabbitai]

🧩 Analysis chain

🏁 Script executed:

#!/bin/bash

# Search for percentile-related code
echo "=== Searching for percentile references ==="
rg -i "percentile" -C 3

echo -e "\n=== Searching for vibration analysis code ==="
rg -i "vibration.*analysis" -C 3

echo -e "\n=== Search for vibration and percentile together ==="
rg -i "vibration.*percentile|percentile.*vibration" -C 5

echo -e "\n=== List all files with 'percentile' in name or content ==="
fd -i percentile

echo -e "\n=== Check for vibration-related plotting functions ==="
rg -l "vibration" src/plot_functions/

echo -e "\n=== Look for statistical analysis functions ==="
ls -la src/data_analysis/ | grep -i "stat\|percentile\|vibr"

Length of output: 3360

---

🏁 Script executed:

#!/bin/bash

# Search more broadly for vibration-related features
echo "=== Search for any vibration-related code ==="
rg -i "vibr" -C 2

echo -e "\n=== Search for noise analysis ==="
rg -i "noise.*analys|analys.*noise" -C 2

echo -e "\n=== Check for gyro noise or vibration metrics ==="
rg -i "gyro.*noise|noise.*gyro|vibration" -C 2

echo -e "\n=== List all data analysis files ==="
ls -la src/data_analysis/

echo -e "\n=== Check README for feature list ==="
grep -i "vibration\|noise\|percentile" README.md -C 3 || echo "No vibration/noise/percentile mentions in README"

echo -e "\n=== Search issues for vibration analysis context ==="
echo "Checking if there are any references to vibration analysis in recent development"

Length of output: 2251

---

🌐 Web query:

vibration percentile analysis flight logs drone betaflight

💡 Result:

Short guide — how to do vibration‑percentile analysis on Betaflight Blackbox logs and what it tells you

What it is

• Vibration‑percentile analysis summarizes how often vibration magnitude (or PSD power at a given frequency) exceeds levels of interest (e.g., median, 90th, 95th percentiles). It helps you quantify worst‑case vibration exposure and compare changes after hardware/tuning fixes. [1][4]

Quick steps

1. Record Blackbox at a high rate (1–2 kHz recommended for vibration/PSD work). [5]
2. Convert the .blf/.bblog to CSV with blackbox_decode from betaflight/blackbox-tools. [2][3]
3. Choose your metric(s):
   • Time-domain: gyro axis values or gyro magnitude (sqrt(gyro_x^2+gyro_y^2+gyro_z^2)) or RMS over short windows (e.g., 10–50 ms).
   • Frequency‑domain: PSD magnitude at frequencies of interest (use Blackbox Explorer PSD view or compute FFT per window). [4]
4. Compute percentiles (e.g., 50th, 90th, 95th, 99th) across your chosen metric over the flight. Higher percentiles show the tail/worst vibration events.
5. Interpret: high gyro RMS/percentiles → mechanical/prop/arm resonance or soft‑mounting issues; spikes at discrete frequencies (from PSD) → motor/prop harmonics or frame resonance. Use PSD vs RPM to correlate to motor harmonics. [4][1]

Minimal example (after converting to CSV)

• Python/pandas example to compute percentiles of gyro magnitude:

import pandas as pd
df = pd.read_csv('log.csv') # output from blackbox_decode
df['gyro_mag'] = (df.gx2 + df.gy2 + df.gz**2)**0.5
for p in [50,90,95,99]:
print(f'{p}th percentile:', df['gyro_mag'].quantile(p/100))

Tips and tools

• Use Blackbox Explorer for interactive PSD, PSD vs RPM heatmaps and quick visual checks. [1][4]
• Use blackbox_decode (blackbox-tools) to export CSV for custom analysis and batch processing. [2][3]
• Record multiple representative flights (hover, aggressive, full‑throttle) and compare percentile changes after fixes. [5]

References
[1] Betaflight Blackbox log docs and Blackbox Explorer.
[2] betaflight/blackbox-tools (blackbox_decode).
[3] Betaflight Blackbox internals / logging details.
[4] Blackbox Explorer PSD guide (PSD & PSD vs RPM).
[5] Practical Blackbox settings and recommendations (Oscar Liang).

---

🌐 Web query:

gyro vibration percentile statistics RMS analysis validity flight data

💡 Result:

Short answer: both RMS and percentile statistics are valid and complementary for evaluating gyro vibration in flight data — RMS (or PSD-integrated RMS) quantifies broad-band energy/noise (good for sensor noise and fatigue/energy metrics), while percentiles (e.g., 90th/95th/99th of instantaneous magnitude or short-window RMS) expose non‑Gaussian peaks, duty‑cycle effects and shocks that RMS averages away. Use PSD/spectral analysis plus time-domain percentiles to get a complete, valid picture, and always check sensor bandwidth, noise density, sampling, filtering, temperature drift and mounting. [1][2][3]

Practical checklist (concise):

• Use PSD to compute band‑limited RMS that matches sensor bandwidth / mission band. This ties to manufacturer noise density. [1][2]
• Compute short-window RMS (e.g., 1–10 ms or relevant structural window) and report percentiles (50/90/95/99) of that short‑window RMS to capture transient/shock content. Percentiles of instantaneous magnitude (or envelope) are useful for duty cycle. [2][3]
• Inspect time series and spectrograms for narrowband tones (rotors, harmonics); RMS hides tonal peaks — percentiles and crest factor reveal them. [1][3]
• Validate results against sensor specs: noise density, scale‑factor drift vs temperature, bandwidth, cross‑axis sensitivity, and aliasing from sampling. If sensor noise floor ~ signal RMS, interpretation is limited. [1]
• Use complementary metrics where needed: PSD, kurtosis/crest factor, peak counts, and cumulative RMS vs frequency (cumulative PSD) for damage/response prediction. [2][3]
• When reporting, state window length, filter bandwidth, sampling rate, percentile definition (e.g., 95th of 100 ms RMS), and any calibration/temperature compensation. This ensures reproducibility and validity. [1][2]

Key references:

• Flight Test Engineering Reference Handbook — motion & vibration analysis (PSD↔RMS conversions, windowing and conversion guidance). [2]
• Applied Sciences (2024) analysis of MEMS accelerometers/gyros in flight recorders — compares RMS, PSD, temperature effects and sensor limits. [1]
• Analog Devices technical article — estimating gyro error from vibration and practical considerations (g‑sensitivity, duty cycle, noise). [3]

If you want, tell me the sensor model, sampling rate, and a short excerpt of your flight time series (or sample PSD/RMS numbers) and I’ll give concrete percentile windows, filters and a validity assessment.

---

🌐 Web query:

Betaflight gyro noise floor vibration analysis spectrum vs time domain

💡 Result:

Short summary — what each domain shows and how to use them in Betaflight:

• Time domain (gyro trace, motor/PIDsum traces): shows when impulses/oscillations occur, how big they are in actual gyro units, and whether filtering or PID changes reduce visible shakes. Use Blackbox to view raw gyro, D-term and motor traces; inspect waveform shape and envelope. Useful for detecting transient events, resonances, and overall control behaviour. [2][5]

• Frequency domain / spectrum (FFT / PSD / spectrogram): shows energy vs frequency — isolates motor/prop/ESC harmonics, frame resonances, and where noise-floor sits. Use Blackbox Explorer’s Power Spectral Density or spectrogram views to find narrow peaks (RPM harmonics) and broadband noise. PSD also lets you plot PSD vs RPM or throttle to confirm RPM-coupled peaks. Interpret peaks (Hz) and their growth into D-term/PIDsum to decide notches or LPFs. [1][2][3]

Practical workflow (concise steps)

1. Log at full debug (raw gyro traces) and high sample rate (PID loop / Blackbox ≥ gyro rate, e.g. 1 kHz logging for 4k PID). Enable debug_mode(notch) to capture pre/post-notch traces. Fly a steady hover and a throttle sweep. [2]
2. Inspect time-domain gyro trace: note amplitude, axis, and when shakes appear (transient vs continuous). If shakes look like narrow-band sinusoidal, expect RPM-related source; broadband fuzz suggests electrical noise or flutter. [2][5]
3. Open PSD / spectrogram for the same traces: find dominant peaks and their frequencies, check PSD vs RPM/throttle to see which peaks track motor RPM (motor/prop). Peaks that don’t track RPM likely frame/soft-mount resonance or electronic noise. Measure the noise floor level and relative peak amplitude — important for filter design. [1][2]
4. Check D-term and PIDsum spectra vs gyro: D-term amplifies high‑freq noise — any remaining gyro peaks will be larger in D. Use this to choose D filters (PT1 + notch). Fix gyro filters first, then D filters. [3]
5. Apply fixes iteratively: add RPM filters/notches for motor harmonics, adjust dyn_lpf/dyn_notch or fixed notches, raise gyro LPF only as needed (beware added delay), soften mounts or add mass for structural resonance, and add power cap/cleaning for electrical noise. Re-log and re-check PSD and time traces. [3][4][5]

Key interpretation tips (quick)

• RPM-correlated narrow lines → motor/prop harmonics → use RPM filter or notches. [1][2]
• Narrow peak not tracking RPM but at low frequency → frame or soft-mount resonance → change mounting stiffness/mass or move notch to that frequency. [4]
• Broadband increase or spikes at high freq → ESC/BLHeli/drive electrical noise or aliasing; check ESC PWM frequency, add capacitors, try different ESC firmware/refresh rates. [5]
• Always compare raw gyro vs post-filter vs D-term vs motor/PIDsum to see how noise propagates through control loop. [2][3]

References

• Blackbox Explorer: Power spectral density charts (Betaflight docs). [1]
• Blackbox logging and analyzer usage (Betaflight docs — debug traces, logging rate, analyzer). [2]
• Gyro and D-term filtering guidance (Betaflight tuning notes). [3]
• Soft-mounting and vibration reduction (Betaflight docs). [4]
• Practical Blackbox analysis tutorial (Oscar Liang). [5]

If you want, I can: (A) walk through your Blackbox log and mark the dominant peaks and recommended filter/notch frequencies, or (B) give a one-page checklist tailored to 3", 5" or cinewhoop builds. Which do you want?

---

@nerdCopter, after deep research, I can confirm that Vibration Percentile Analysis IS useful and valid for real Betaflight flight logs—unlike Bode plots.

Key Finding: It Doesn't Exist in Your Codebase Yet

I searched the entire codebase and found no "Vibration Percentile Analysis" feature. The only percentile usage is for Y-axis scaling in plots (5th/95th percentiles to exclude outliers for better visualization).

Why It's Valid (Unlike Bode Plots)

Fundamental difference:

Feature	Requires Controlled Input?	Works with Normal Flights?
Bode Plots	✅ Yes (chirp/PRBS)	❌ No (low coherence, unreliable)
Vibration Percentile Analysis	❌ No	✅ Yes (any operational data)

What Makes Percentile Analysis Valid

Based on research from flight test engineering and sensor analysis literature:

1. No special input required - Works with ANY flight log (hover, acro, freestyle)
2. Captures what PSD misses - Transient shocks, worst-case peaks, duty-cycle effects that RMS/spectral analysis averages away
3. Complementary to spectrum analysis - Your existing PSD/spectrum plots show frequency content; percentiles show magnitude distribution and tail behavior
4. Industry standard - Compute percentiles (50th, 90th, 95th, 99th) of:
   • Gyro magnitude: sqrt(gyro_x² + gyro_y² + gyro_z²)
   • Short-window RMS (e.g., 10-50ms windows)
   • PSD power at specific frequencies

Practical Use Cases

• Before/after hardware changes: Compare 95th percentile gyro magnitude to quantify vibration reduction
• Worst-case analysis: 99th percentile shows tail events (hard maneuvers, resonances)
• Flight mode comparison: Compare cruise vs acro vibration exposure

Would You Like This Implemented?

This would be a straightforward addition that complements your existing spectrum analysis features. It could provide:

• Per-axis gyro magnitude percentiles
• Overall magnitude percentiles
• Time-series plots showing percentile bands
• Comparison metrics across multiple logs

Should I create an implementation plan for this feature?

---

🧠 Learnings used

Learnt from: nerdCopter
Repo: nerdCopter/BlackBox_CSV_Render PR: 0
File: :0-0
Timestamp: 2026-01-08T14:48:17.022Z
Learning: For the BlackBox_CSV_Render project, Bode plot transfer function estimation (PR `#118`) is mathematically correct but fundamentally incompatible with normal operational flight logs due to low coherence (requires controlled system-identification inputs like chirp/PRBS on secured craft). Normal flights produce unreliable results due to nonlinearities, closed-loop conditions, and nonstationary maneuvers. The feature should be opt-in only (via --bode flag) with documentation clarifying it's for formal system-identification experiments, not routine flight log analysis. Existing PSD/spectrum plots are the appropriate tool for operational flight logs.

---

Learnt from: nerdCopter
Repo: nerdCopter/BlackBox_CSV_Render PR: 0
File: :0-0
Timestamp: 2025-09-09T19:10:32.659Z
Learning: [monologue]
[task review_after_revert_1/1 analyze_current_state_1/1]
The user has reverted to commit 919100d and force pushed, expressing frustration with the ongoing issues with the measured filter response analysis. They're requesting a full review to understand the current state and provide proper guidance.

[observations]
- User has been dealing with persistent issues with measured cutoff frequencies showing unrealistic values (first 15Hz, then stuck at 50Hz)
- Multiple attempts at fixes have been made but haven't resolved the core issue
- User has now reverted to an earlier commit (919100d) and wants a comprehensive review
- From the learning, I know the user wants the measured filter response analysis to be purely data-driven without artificial frequency constraints
- They want to measure actual cutoff frequencies from unfiltered vs filtered spectral data, even if results show unexpectedly low cutoffs

[analysis]
The user is frustrated because we've been applying incremental fixes that haven't resolved the fundamental issue. A full review means I need to:
1. Examine the current state of the code after the revert
2. Understand what commit 919100d contains
3. Provide a comprehensive analysis of the measured filter response implementation
4. Identify the root causes of the cutoff frequency detection issues
5. Provide a clear, working solution

The fact they reverted suggests the incremental approach wasn't working and we need to start fresh with a proper understanding of the codebase.

[actions]
- Use run_scripts to examine the current code state after revert
- Analyze the complete filter response measurement implementation
- Identify fundamental issues with the approach
- Provide comprehensive recommendations for fixing the system
</monologue>

I understand your frustration! Let me conduct a comprehensive review of the current codebase state after the revert to properly understand and fix the measured filter response analysis.

---

🏁 Script executed:
```
#!/bin/bash
# Check current commit state and what 919100d contains
echo "=== Current commit information ==="
git log --oneline -5

echo -e "\n=== Examine the current measured filter response implementation ==="
cat src/data_analysis/filter_response.rs | head -50

echo -e "\n=== Check what constants are currently defined ==="
grep -A 30 "Measured filter response analysis constants" src/constants.rs || echo "No measured filter constants found"

echo -e "\n=== Look for measure_filter_response function ==="
rg -n -A 20 "pub fn measure_filter_response" src/data_analysis/filter_response.rs || echo "measure_filter_response function not found"

echo -e "\n=== Check if --measure flag exists in main.rs ==="
rg -n -C 5 "measure" src/main.rs || echo "No measure flag found"
```

Length of output: 8172

---

🏁 Script executed:
```
#!/bin/bash
# Let's examine the complete measure_filter_response implementation
echo "=== Complete measure_filter_response function ==="
sed -n '/pub fn measure_filter_response/,/^}/p' src/data_analysis/filter_response.rs

echo -e "\n=== Look at how it's called from the plotting function ==="
rg -n -A 10 -B 10 "measure_filter_response" src/plot_functions/plot_gyro_spectrums.rs

echo -e "\n=== Check the cutoff detection function ==="
rg -n -A 15 -B 5 "find_cutoff_frequency" src/data_analysis/filter_response.rs

echo -e "\n=== Check MeasuredFilterResponse struct ==="
rg -n -A 10 -B 5 "struct MeasuredFilterResponse" src/data_analysis/filter_response.rs
```

Length of output: 5586

---

🏁 Script executed:
```
#!/bin/bash
# Let's see the complete find_cutoff_frequency function
echo "=== Complete find_cutoff_frequency function ==="
sed -n '/fn find_cutoff_frequency/,/^}/p' src/data_analysis/filter_response.rs

echo -e "\n=== Let's also see how the transfer function data is filtered ==="
rg -n -A 5 -B 5 "FILTER_ANALYSIS_MIN_FREQ_HZ.*FILTER_ANALYSIS_MAX_FREQ_HZ" src/data_analysis/filter_response.rs

echo -e "\n=== Check what the current constants are ==="
grep -n "FILTER_ANALYSIS_MIN_FREQ_HZ\|FILTER_ANALYSIS_MAX_FREQ_HZ" src/constants.rs

echo -e "\n=== Let's see what the actual frequency range covers for typical data ==="
echo "Current constants:"
echo "MIN: 10Hz, MAX: 500Hz"
echo "This means analysis only looks at 10-500Hz range"
```

Length of output: 1279

---

<!-- [monologue]
[task full_review_1/1 analyze_core_issues_1/1]
Now I can see the fundamental problem clearly after this full review. The current implementation has several critical flaws:

[observations]
- The cutoff detection algorithm is completely wrong - it's using a "closest match" approach instead of looking for an actual -3dB crossing
- The `find_cutoff_frequency` function finds the frequency point with the magnitude ratio CLOSEST to 0.707, not where it actually crosses 0.707
- This means if the transfer function never actually crosses 0.707 properly, it just picks whatever frequency happens to be numerically closest
- The transfer function calculation looks reasonable (no artificial caps), but the detection logic is fundamentally flawed
- The frequency filtering (10-500Hz) is reasonable
- The user wants a purely data-driven approach that measures actual cutoff frequencies from the data

[analysis]
The reason the user has been getting strange results (15Hz, then 50Hz) is because:
1. The algorithm isn't looking for a proper -3dB crossing point
2. Instead, it's just finding whatever frequency point has a magnitude ratio closest to 0.707
3. If the data is noisy or the transfer function doesn't behave as expected, this "closest match" approach will pick essentially random points
4. This is not how filter cutoff detection should work at all

For proper cutoff detection, we need to:
1. Look for where the transfer function actually crosses from above 0.707 to below 0.707 (or vice versa)
2. Use interpolation to find the exact crossing point
3. Handle cases where there might be multiple crossings or no clear crossing

The current "closest match" algorithm is fundamentally broken for this purpose.

[actions]
- Explain the core issue with the current cutoff detection algorithm
- Provide a complete, correct implementation that actually looks for -3dB crossings
- Give clear recommendations for a proper data-driven approach
- Acknowledge the user's frustration and provide a definitive solution
</monologue>

nerdCopter I completely understand your frustration! After this full review, I've identified the **fundamental flaw** in the current implementation. The issue isn't with windowing, frequency ranges, or caps - it's with the **cutoff detection algorithm itself**.

## The Core Problem

The current `find_cutoff_frequency` function is completely wrong for filter analysis:

```rust
// CURRENT (BROKEN) ALGORITHM:
for &(freq, ratio) in transfer_function {
    let diff = (ratio - target_ratio).abs();
    if diff < best_diff {  // ← This just finds the CLOSEST value to 0.707
        best_diff = diff;   //   NOT an actual crossing!
        best_freq = freq;
    }
}
```

**This algorithm doesn't detect cutoff frequencies - it just finds whatever frequency happens to have a magnitude ratio numerically closest to 0.707!** No wonder you're getting random results like 15Hz or 50Hz.

## The Correct Solution

Here's a complete, proper implementation for measured filter response analysis:

---

🏁 Script executed:
```
#!/bin/bash
# Create a corrected implementation
cat > fixed_filter_response.rs << 'EOF'
// Replace the find_cutoff_frequency function with this corrected version:

fn find_cutoff_frequency(
    transfer_function: &[(f64, f64)],
    target_ratio: f64,
) -> Result<f64, Box<dyn std::error::Error>> {
    if transfer_function.len() < 2 {
        return Err("Need at least 2 points for cutoff detection".into());
    }

    // Sort transfer function by frequency to ensure proper ordering
    let mut tf_sorted: Vec<(f64, f64)> = transfer_function.to_vec();
    tf_sorted.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));

    // Look for where the transfer function crosses the target ratio (0.707)
    // We want the frequency where the response goes from above to below the threshold
    for i in 1..tf_sorted.len() {
        let (f1, r1) = tf_sorted[i - 1];
        let (f2, r2) = tf_sorted[i];

        // Check for valid data
        if !f1.is_finite() || !f2.is_finite() || !r1.is_finite() || !r2.is_finite() {
            continue;
        }

        // Look for downward crossing: previous point above threshold, current below
        if r1 >= target_ratio && r2 <= target_ratio {
            // Linear interpolation to find exact crossing point
            if (r2 - r1).abs() < 1e-12 {
                // Avoid division by nearly zero
                return Ok(f1);
            }
            
            let t = (target_ratio - r1) / (r2 - r1);
            let cutoff_freq = f1 + t * (f2 - f1);
            
            if cutoff_freq > 0.0 && cutoff_freq.is_finite() {
                return Ok(cutoff_freq);
            }
        }
    }

    // If no clear downward crossing found, try to find the best approximation
    // Look for the point closest to target_ratio, but only if it's reasonable
    let mut best_freq = None;
    let mut best_diff = f64::INFINITY;
    
    for &(freq, ratio) in &tf_sorted {
        if freq > 10.0 && freq < 1000.0 && ratio > 0.1 && ratio < 2.0 {
            let diff = (ratio - target_ratio).abs();
            if diff < best_diff {
                best_diff = diff;
                best_freq = Some(freq);
            }
        }
    }
    
    if let Some(freq) = best_freq {
        if best_diff < 0.2 { // Only accept if reasonably close to target
            return Ok(freq);
        }
    }

    Err("Could not find cutoff frequency crossing".into())
}

// Also improve the measure_filter_response function:

pub fn measure_filter_response(
    unfiltered_spectrum: &[(f64, f64)],
    filtered_spectrum: &[(f64, f64)],
    _sample_rate: f64,
) -> Result<MeasuredFilterResponse, Box<dyn std::error::Error>> {
    // Validate input data
    if filtered_spectrum.len() != unfiltered_spectrum.len() {
        return Err("Filtered and unfiltered spectra must have same length".into());
    }
    if filtered_spectrum.len() < MIN_SPECTRUM_POINTS_FOR_ANALYSIS {
        return Err(format!(
            "Need at least {} frequency points for reliable analysis",
            MIN_SPECTRUM_POINTS_FOR_ANALYSIS
        ).into());
    }

    // Build transfer function with better validation
    let mut transfer_function = Vec::new();
    let mut skipped_count = 0;

    for i in 0..filtered_spectrum.len() {
        let (freq, filt_mag) = filtered_spectrum[i];
        let (_, unfilt_mag) = unfiltered_spectrum[i];

        // Apply frequency bounds and validate data
        if freq >= FILTER_ANALYSIS_MIN_FREQ_HZ
            && freq <= FILTER_ANALYSIS_MAX_FREQ_HZ
            && freq.is_finite()
            && filt_mag.is_finite()
            && unfilt_mag.is_finite()
            && unfilt_mag > 1e-12  // Avoid division by very small numbers
            && filt_mag >= 0.0
        {
            let ratio = filt_mag / unfilt_mag;
            
            // Allow ratios > 1.0 but reject extreme outliers
            if ratio.is_finite() && ratio > 0.0 && ratio < 10.0 {
                transfer_function.push((freq, ratio));
            } else {
                skipped_count += 1;
            }
        } else {
            skipped_count += 1;
        }
    }

    if transfer_function.len() < 10 {
        return Err(format!(
            "Insufficient valid data points: {} (skipped {})",
            transfer_function.len(), skipped_count
        ).into());
    }

    // Find cutoff frequency using proper crossing detection
    let cutoff_hz = find_cutoff_frequency(&transfer_function, CUTOFF_MAGNITUDE_RATIO)?;

    // Estimate filter order (keep existing implementation)
    let filter_order = estimate_filter_order_from_transfer(&transfer_function, cutoff_hz)
        .unwrap_or(DEFAULT_FILTER_ORDER);

    // Calculate confidence (keep existing implementation)  
    let confidence = calculate_measurement_confidence(&transfer_function);

    Ok(MeasuredFilterResponse {
        cutoff_hz,
        filter_order,
        confidence,
        transfer_function,
    })
}
EOF

echo "Created corrected implementation - the key changes:"
echo "1. Proper crossing detection instead of 'closest match'"
echo "2. Linear interpolation for accurate cutoff frequency"
echo "3. Better data validation and outlier rejection"
echo "4. Fallback logic for edge cases"
echo ""
echo "This should give you realistic, data-driven cutoff measurements!"
```

Length of output: 309

---



</details>

nerdCopter I completely understand your frustration! After this full review, I've identified the **fundamental flaw** in the current implementation that explains all your issues.

## The Core Problem

The current `find_cutoff_frequency` function is completely wrong for filter analysis:

```rust
// CURRENT (BROKEN) ALGORITHM:
for &(freq, ratio) in transfer_function {
    let diff = (ratio - target_ratio).abs();  // ← Just finds CLOSEST to 0.707
    if diff < best_diff {                     //   NOT an actual crossing!
        best_freq = freq;
    }
}
```

**This doesn't detect cutoff frequencies - it just finds whatever frequency happens to be numerically closest to 0.707!** That's why you've been getting random results.

## Complete Solution

Replace the current `find_cutoff_frequency` and `measure_filter_response` functions in `src/data_analysis/filter_response.rs` with this corrected implementation:

```rust
fn find_cutoff_frequency(
    transfer_function: &[(f64, f64)],
    target_ratio: f64,
) -> Result<f64, Box<dyn std::error::Error>> {
    if transfer_function.len() < 2 {
        return Err("Need at least 2 points for cutoff detection".into());
    }

    // Sort by frequency to ensure proper ordering
    let mut tf_sorted: Vec<(f64, f64)> = transfer_function.to_vec();
    tf_sorted.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));

    // Look for actual -3dB crossing (above to below threshold)
    for i in 1..tf_sorted.len() {
        let (f1, r1) = tf_sorted[i - 1];
        let (f2, r2) = tf_sorted[i];

        if !f1.is_finite() || !f2.is_finite() || !r1.is_finite() || !r2.is_finite() {
            continue;
        }

        // Look for downward crossing: above → below threshold
        if r1 >= target_ratio && r2 <= target_ratio {
            if (r2 - r1).abs() < 1e-12 {
                return Ok(f1);
            }
            
            // Linear interpolation for exact crossing
            let t = (target_ratio - r1) / (r2 - r1);
            let cutoff_freq = f1 + t * (f2 - f1);
            
            if cutoff_freq > 0.0 && cutoff_freq.is_finite() {
                return Ok(cutoff_freq);
            }
        }
    }

    // Fallback: find reasonable approximation
    let mut best_freq = None;
    let mut best_diff = f64::INFINITY;
    
    for &(freq, ratio) in &tf_sorted {
        if freq > 10.0 && freq < 1000.0 && ratio > 0.1 && ratio < 2.0 {
            let diff = (ratio - target_ratio).abs();
            if diff < best_diff {
                best_diff = diff;
                best_freq = Some(freq);
            }
        }
    }
    
    if let Some(freq) = best_freq {
        if best_diff < 0.2 { // Only accept if reasonably close
            return Ok(freq);
        }
    }

    Err("Could not find cutoff frequency crossing".into())
}

pub fn measure_filter_response(
    unfiltered_spectrum: &[(f64, f64)],
    filtered_spectrum: &[(f64, f64)],
    _sample_rate: f64,
) -> Result<MeasuredFilterResponse, Box<dyn std::error::Error>> {
    if filtered_spectrum.len() != unfiltered_spectrum.len() {
        return Err("Filtered and unfiltered spectra must have same length".into());
    }
    if filtered_spectrum.len() < MIN_SPECTRUM_POINTS_FOR_ANALYSIS {
        return Err(format!(
            "Need at least {} frequency points for reliable analysis",
            MIN_SPECTRUM_POINTS_FOR_ANALYSIS
        ).into());
    }

    // Build transfer function with proper validation
    let mut transfer_function = Vec::new();

    for i in 0..filtered_spectrum.len() {
        let (freq, filt_mag) = filtered_spectrum[i];
        let (_, unfilt_mag) = unfiltered_spectrum[i];

        if freq >= FILTER_ANALYSIS_MIN_FREQ_HZ
            && freq <= FILTER_ANALYSIS_MAX_FREQ_HZ
            && freq.is_finite()
            && filt_mag.is_finite() 
            && unfilt_mag.is_finite()
            && unfilt_mag > 1e-12  // Avoid division by tiny numbers
            && filt_mag >= 0.0
        {
            let ratio = filt_mag / unfilt_mag;
            
            // Allow ratios > 1.0 but reject extreme outliers
            if ratio.is_finite() && ratio > 0.0 && ratio < 10.0 {
                transfer_function.push((freq, ratio));
            }
        }
    }

    if transfer_function.len() < 10 {
        return Err("Insufficient valid data points".into());
    }

    // Use proper crossing detection
    let cutoff_hz = find_cutoff_frequency(&transfer_function, CUTOFF_MAGNITUDE_RATIO)?;

    let filter_order = estimate_filter_order_from_transfer(&transfer_function, cutoff_hz)
        .unwrap_or(DEFAULT_FILTER_ORDER);

    let confidence = calculate_measurement_confidence(&transfer_function);

    Ok(MeasuredFilterResponse {
        cutoff_hz,
        filter_order,
        confidence,
        transfer_function,
    })
}
```

## What This Fixes

1. **Proper crossing detection** - Looks for where transfer function actually crosses -3dB threshold
2. **Linear interpolation** - Gives accurate cutoff frequency between data points
3. **Data-driven approach** - Measures what's actually in your flight data
4. **Robust validation** - Handles edge cases and outliers properly

This should finally give you realistic, measured cutoff frequencies that reflect the actual filtering behavior in your flight logs, as you originally intended!

<!-- 
The user wants the measured filter response analysis to be purely data-driven without artificial frequency constraints. They want to measure actual cutoff frequencies from unfiltered vs filtered spectral data, even if the results show unexpectedly low cutoffs like 15Hz. The measurement system should report what it finds in the real data, not what we expect it should find.

---

Learnt from: nerdCopter
Repo: nerdCopter/BlackBox_CSV_Render PR: 71
File: src/plot_functions/plot_gyro_vs_unfilt.rs:19-27
Timestamp: 2025-10-19T18:24:42.528Z
Learning: In the BlackBox_CSV_Render project, when adding debug mode annotations to gyro and D-term plots, nerdCopter prefers to include the debug mode tag (e.g., "[Debug=MODE]") only in the legend labels, not in the plot titles.

---

Learnt from: nerdCopter
Repo: nerdCopter/BlackBox_CSV_Render PR: 71
File: src/main.rs:888-909
Timestamp: 2025-10-19T18:25:15.413Z
Learning: In the BlackBox_CSV_Render project, gyroUnfilt headers are always present for all three axes or completely absent (with all three axes falling back to debug data). Per-axis fallback detection is not needed—the all-or-nothing approach in main.rs correctly reflects real-world log file behavior.

---

Learnt from: nerdCopter
Repo: nerdCopter/BlackBox_CSV_Render PR: 60
File: src/data_analysis/filter_response.rs:966-1021
Timestamp: 2025-10-05T20:37:53.448Z
Learning: For the BlackBox_CSV_Render measured filter response analysis, the attenuation-based approach `abs(unfilt_mag - filt_mag)` works better than the transfer function ratio approach `filt_mag / unfilt_mag` for real flight data. The ratio approach produces spurious cutoff detections (e.g., 47Hz) due to numerical instability and noisy crossings, while the attenuation approach reliably finds where significant filtering begins.

---

Learnt from: nerdCopter
Repo: nerdCopter/BlackBox_CSV_Render PR: 55
File: LICENSE_COMMERCIAL:1-7
Timestamp: 2025-09-10T13:56:40.607Z
Learning: nerdCopter prefers to use minimal commercial license templates that are consistent across his projects, similar to what was used in the bbl_parser repository PR `#2`.

---

Learnt from: nerdCopter
Repo: nerdCopter/BlackBox_CSV_Render PR: 18
File: src/data_analysis/calc_step_response.rs:0-0
Timestamp: 2025-06-03T18:01:54.978Z
Learning: The user nerdCopter prefers to defer unit testing implementation when suggested for the BlackBox_CSV_Render project.

[coderabbitai]

Coding Plan

Approach

Add a new Vibration Percentile Analysis plot that displays Empirical Cumulative Distribution Functions (ECDFs) for gyro vibration magnitude, showing statistical distribution of vibrations at key percentiles (50th, 90th, 95th, 99th) and comparing filtered vs unfiltered gyro data to identify consistent vs transient noise patterns.

Observations

The project is a Rust-based Betaflight flight log analyzer with 16 existing plot types in src/plot_functions/. Each plot follows a consistent pattern: a public function accepting log_data: &[LogRowData], root_name: &str, and sample_rate: Option<f64>, returning Result<(), Box<dyn Error>>. The framework uses the plotters crate with helper functions in plot_framework.rs for common visualizations.

Gyro data is available as both filtered (gyro) and unfiltered (gyro_unfilt) 3-axis arrays in LogRowData. Existing plots use percentile-based scaling inline (e.g., in plot_setpoint_vs_gyro.rs), and there are established color conventions: LIGHTBLUE for filtered gyro, AMBER for unfiltered.

Assumptions

Assumption 1: Vibration magnitude metric calculation

Options Considered:

• Per-axis absolute magnitude |gyro[axis]| - preserves directional information, simpler calculation
• Combined Euclidean magnitude sqrt(x²+y²+z²) - single scalar per sample, loses directional info
• RMS over time windows - reduces noise but adds complexity

Chosen Option: 1

Rationale: Per-axis absolute magnitude in a stacked 3-axis plot format. This aligns with existing codebase patterns (most plots use Roll/Pitch/Yaw stacked views), preserves directional vibration information useful for identifying axis-specific issues, and is consistent with how users diagnose Betaflight tuning problems.

Assumption 2: ECDF visualization approach

Options Considered:

• True step plot using duplicated points in LineSeries
• Smooth line plot connecting sorted cumulative probability points
• Custom PathElement-based step rendering

Chosen Option: 2

Rationale: Smooth line connecting ECDF points using LineSeries. With thousands of data points, a smooth line connecting sorted cumulative probabilities will appear step-like naturally and is simpler to implement using existing patterns.

💡 User Tips

Regenerate the plan with different choices with @coderabbitai <feedback>.

Implementation Steps

Phase 1: Create the Vibration Percentile Plot Function

This phase implements the core ECDF computation and visualization logic as a new plot module following established codebase patterns.

Task 1: Create the plot module file

Create src/plot_functions/plot_vibration_percentile.rs implementing the main plot function with the standard signature pattern.

• Extract gyro (filtered) and gyro_unfilt (unfiltered) data for all three axes from log_data
• For each axis, compute absolute values of gyro readings, filter out NaN/invalid values, and sort the resulting vectors
• Compute ECDF by assigning cumulative probability (i / n) to each sorted value
• Calculate key percentile values (50th, 90th, 95th, 99th) from the sorted data for annotation

Task 2: Implement ECDF visualization with filtered/unfiltered comparison

Build the visualization using plotters within the same module, rendering per-axis ECDFs in a stacked 3-axis format.

• Use draw_stacked_plot pattern or create a custom 3-row chart layout similar to existing stacked plots
• For each axis chart, draw two LineSeries: filtered gyro ECDF (LIGHTBLUE) and unfiltered gyro ECDF (AMBER)
• X-axis represents gyro magnitude (deg/s), Y-axis represents cumulative probability (0.0 to 1.0)
• Add legend distinguishing "Filtered" and "Unfiltered" series
• Use percentile-based X-axis scaling (99th percentile with headroom) to avoid outliers dominating the view

Task 3: Add percentile annotations to the plot

Annotate key percentiles directly on each axis chart to highlight important statistical thresholds.

• Draw horizontal reference lines at Y = 0.5, 0.9, 0.95, 0.99 with subtle styling (dashed or light gray)
• Add text labels indicating the corresponding X-value (gyro magnitude) at each percentile intersection
• For filtered data, annotate the 95th and 99th percentile values prominently as these are most diagnostic
• Use existing font/styling constants from constants.rs for consistency

Task 4: Generate output file with consistent naming

Construct the output filename following the established naming convention.

• Output filename pattern: {root_name}_Gyro_Vibration_Percentile.png
• Use standard plot dimensions from constants (2560×1400)
• Include informative title: "Gyro Vibration Percentile Distribution (Roll/Pitch/Yaw)"

🤖 Prompt for AI agents

Implement the core ECDF computation and visualization logic as a new plot
module. Create a vibration percentile plot that computes and displays ECDFs for
gyro data across all three axes.

Create `src/plot_functions/plot_vibration_percentile.rs` with the standard
signature (log_data, root_name, sample_rate).

Extract filtered and unfiltered gyro data for all axes from log_data
For each axis, compute absolute values, filter invalid values, and sort
Compute ECDF by assigning cumulative probability (i/n) to sorted values
Calculate 50th, 90th, 95th, 99th percentile values

Implement visualization using plotters with 3-row stacked layout
(Roll/Pitch/Yaw)
For each axis: draw two LineSeries (filtered in LIGHTBLUE, unfiltered in AMBER)
X-axis: gyro magnitude (deg/s), Y-axis: cumulative probability (0.0-1.0)
Use percentile-based X-axis scaling (99th percentile with headroom)
Add legend for filtered/unfiltered series

Add percentile annotations:
Horizontal reference lines at Y = 0.5, 0.9, 0.95, 0.99 with subtle styling
Text labels showing X-values at percentile intersections
Prominently annotate 95th and 99th percentile values for filtered data

Generate output with filename pattern:
`{root_name}_Gyro_Vibration_Percentile.png`
Use standard plot dimensions (2560×1400)
Set title: "Gyro Vibration Percentile Distribution (Roll/Pitch/Yaw)"

Phase 2: Register and Integrate the Plot

This phase wires the new plot into the application's execution flow following the established registration pattern.

Task 1: Add module declaration

Register the new module in the plot_functions module index.

• Add pub mod plot_vibration_percentile; to src/plot_functions/mod.rs

Task 2: Add PlotConfig field and defaults

Extend the PlotConfig struct to include a toggle for the new plot type.

• Add pub vibration_percentile: bool field to the PlotConfig struct in src/main.rs
• Set default value to true in the Default implementation
• Set to false in the none() implementation

Task 3: Import and invoke the plot function

Add the conditional invocation in main.rs following the established pattern.

• Add import statement: use crate::plot_functions::plot_vibration_percentile::plot_vibration_percentile;
• Add conditional invocation block after the D-term heatmap section (~line 1088): if plot_config.vibration_percentile { plot_vibration_percentile(&all_log_data, &root_name_string, sample_rate)?; }
• Position groups this with other vibration/frequency analysis plots before the final PID Activity summary plot

🤖 Prompt for AI agents

Wire the new vibration percentile plot into the application's execution flow.
Register the plot module and add configuration toggle.

Add module declaration to `src/plot_functions/mod.rs`:
Add `pub mod plot_vibration_percentile;`

Extend PlotConfig in `src/main.rs`:
Add `pub vibration_percentile: bool` field to PlotConfig struct
Set to `true` in Default implementation
Set to `false` in none() implementation

Add plot invocation in main.rs:
Import: `use
crate::plot_functions::plot_vibration_percentile::plot_vibration_percentile;`
Add conditional block after D-term heatmap section: `if
plot_config.vibration_percentile { plot_vibration_percentile(&all_log_data,
&root_name_string, sample_rate)?; }`

---

🚀 Next Steps

🤖 All AI agent prompts combined

Task: 1

Implement the core ECDF computation and visualization logic as a new plot
module. Create a vibration percentile plot that computes and displays ECDFs for
gyro data across all three axes.

Create `src/plot_functions/plot_vibration_percentile.rs` with the standard
signature (log_data, root_name, sample_rate).

Extract filtered and unfiltered gyro data for all axes from log_data
For each axis, compute absolute values, filter invalid values, and sort
Compute ECDF by assigning cumulative probability (i/n) to sorted values
Calculate 50th, 90th, 95th, 99th percentile values

Implement visualization using plotters with 3-row stacked layout
(Roll/Pitch/Yaw)
For each axis: draw two LineSeries (filtered in LIGHTBLUE, unfiltered in AMBER)
X-axis: gyro magnitude (deg/s), Y-axis: cumulative probability (0.0-1.0)
Use percentile-based X-axis scaling (99th percentile with headroom)
Add legend for filtered/unfiltered series

Add percentile annotations:
Horizontal reference lines at Y = 0.5, 0.9, 0.95, 0.99 with subtle styling
Text labels showing X-values at percentile intersections
Prominently annotate 95th and 99th percentile values for filtered data

Generate output with filename pattern:
`{root_name}_Gyro_Vibration_Percentile.png`
Use standard plot dimensions (2560×1400)
Set title: "Gyro Vibration Percentile Distribution (Roll/Pitch/Yaw)"
===============================================================================

Task: 2

Wire the new vibration percentile plot into the application's execution flow.
Register the plot module and add configuration toggle.

Add module declaration to `src/plot_functions/mod.rs`:
Add `pub mod plot_vibration_percentile;`

Extend PlotConfig in `src/main.rs`:
Add `pub vibration_percentile: bool` field to PlotConfig struct
Set to `true` in Default implementation
Set to `false` in none() implementation

Add plot invocation in main.rs:
Import: `use
crate::plot_functions::plot_vibration_percentile::plot_vibration_percentile;`
Add conditional block after D-term heatmap section: `if
plot_config.vibration_percentile { plot_vibration_percentile(&all_log_data,
&root_name_string, sample_rate)?; }`

💡 Iterate on the plan with: @coderabbitai <feedback>

Example Feedback
- `@coderabbitai` You can skip phase 3. Add a simple unit test case for phase 2.
- `@coderabbitai` For assumption 1 go ahead with option 3 and replan.

---

💬 Have feedback or questions? Drop into our discord!