diff --git a/OVERVIEW.md b/OVERVIEW.md
index 6c5d412..eac549a 100644
--- a/OVERVIEW.md
+++ b/OVERVIEW.md
@@ -182,6 +182,11 @@ When `--step` flag is not used, all plots below are generated:
 - **`*_Throttle_Freq_Heatmap_comparative.png`** — System noise characteristics across throttle levels and frequencies
 - **`*_PID_Activity_stacked.png`** — P, I, D term activity over time for each axis (Roll, Pitch, Yaw). Displays all three PID components on the same time-domain plot with unified Y-axis scaling for visual comparison. Each term shows min/avg/max statistics in the legend. Useful for visualizing PID contribution balance during flight and identifying control issues (persistent P-term offset, I-term wind direction, D-term phase lag).
 
+#### Generated Reports
+
+- **`*_report.md`** — Structured markdown flight report written alongside PNGs on every run. Content is assembled from typed result structs returned by each analysis pass — no CSV re-reading. Sections: Metadata (firmware revision, craft name, PIDs, sample rate, gyroUnfilt source warning), Filter Configuration (per-axis LPF1/LPF2/IMUF/Pseudo-Kalman table, Dynamic Notch, RPM filter), PID Tuning P:D ratios, Step Response Analysis (Roll/Pitch: peak value, assessment, setpoint authority, P:D recommendations), Gyro Analysis (per-axis filtering delay with confidence, spectrum peaks), D-Term Analysis (per-axis filtering delay with N/A disambiguation, spectrum peaks), Motor Oscillation table, and relative links to all generated PNGs. Optimal P Estimation and Bode Analysis sections are included when those features produce results.
+
+
 #### P:D Ratio Recommendations
 
 The system provides intelligent P:D tuning recommendations based on step-response peak analysis:
diff --git a/README.md b/README.md
index 20d539e..31248f0 100644
--- a/README.md
+++ b/README.md
@@ -104,6 +104,10 @@ Arguments can be in any order. Wildcards (e.g., *.csv) are shell-expanded and wo
 - `*_Gyro_PSD_Spectrogram_comparative.png` — Gyro spectrogram (PSD vs. time)
 - `*_Throttle_Freq_Heatmap_comparative.png` — Throttle/frequency heatmap analysis
 
+#### Markdown Report (always generated)
+
+- `*_report.md` — Structured flight report written alongside PNGs on every run. Sections: Metadata (firmware, PIDs, sample rate, gyroUnfilt source), Filter Configuration (LPF1/LPF2/IMUF/Pseudo-Kalman table, Dynamic Notch, RPM filter), PID Tuning, Step Response Analysis (Roll/Pitch with P:D assessment and setpoint authority), Gyro Analysis (filtering delay, confidence, spectrum peaks per axis), D-Term Analysis (filtering delay with N/A reason, spectrum peaks), Motor Oscillation, and links to all generated PNGs. Optimal P Estimation and Bode Analysis sections appear when those features are active.
+
 #### Console Output:
 - Current P:D ratio and peak analysis with response assessment
 - Conservative and Moderate tuning recommendations (with D/D-Min/D-Max values)
diff --git a/src/data_analysis/d_term_delay.rs b/src/data_analysis/d_term_delay.rs
index 51216af..5faf591 100644
--- a/src/data_analysis/d_term_delay.rs
+++ b/src/data_analysis/d_term_delay.rs
@@ -12,6 +12,13 @@ use crate::data_analysis::derivative::calculate_derivative;
 use crate::data_analysis::filter_delay;
 use crate::data_input::log_data::LogRowData;
 
+/// Per-axis outcome from D-term delay calculation, including why delay is unavailable.
+pub struct DTermAxisDelay {
+    pub result: Option<filter_delay::DelayResult>,
+    /// Human-readable reason when result is None (None when result is Some).
+    pub na_reason: Option<&'static str>,
+}
+
 /// Calculate filtering delay comparison between unfiltered and filtered D-terms.
 ///
 /// This function computes the delay between the unfiltered D-term (calculated as derivative of gyroUnfilt)
@@ -44,23 +51,36 @@ use crate::data_input::log_data::LogRowData;
 pub fn calculate_d_term_filtering_delay_comparison(
     log_data: &[LogRowData],
     sample_rate: f64,
-) -> Vec<Option<filter_delay::DelayResult>> {
-    // Validate input parameters
+) -> Vec<DTermAxisDelay> {
+    let make_empty = |reason: &'static str| -> Vec<DTermAxisDelay> {
+        (0..AXIS_NAMES.len())
+            .map(|_| DTermAxisDelay {
+                result: None,
+                na_reason: Some(reason),
+            })
+            .collect()
+    };
+
     if !sample_rate.is_finite() || sample_rate <= 0.0 {
         eprintln!(
             "Error: Invalid sample rate ({}) for D-term delay analysis",
             sample_rate
         );
-        return vec![None; AXIS_NAMES.len()];
+        return make_empty("Invalid sample rate");
     }
 
     if log_data.is_empty() {
         eprintln!("Error: Empty log data provided for D-term delay analysis");
-        return vec![None; AXIS_NAMES.len()];
+        return make_empty("No log data");
     }
 
-    // Initialize with None for all axes to preserve axis alignment
-    let mut results: Vec<Option<filter_delay::DelayResult>> = vec![None; AXIS_NAMES.len()];
+    // Initialize with no-data reason; overwritten when data is found
+    let mut results: Vec<DTermAxisDelay> = (0..AXIS_NAMES.len())
+        .map(|_| DTermAxisDelay {
+            result: None,
+            na_reason: Some("No D-term data"),
+        })
+        .collect();
 
     // First, check data availability for diagnosis
     let mut gyro_unfilt_available = [false; 3];
@@ -112,8 +132,12 @@ pub fn calculate_d_term_filtering_delay_comparison(
 
     // Use AXIS_NAMES.iter().enumerate() for consistency with other parts of the codebase
     for (axis_idx, axis_name) in AXIS_NAMES.iter().enumerate() {
-        // Skip if either data type is unavailable
-        if !gyro_unfilt_available[axis_idx] || !d_term_available[axis_idx] {
+        if !gyro_unfilt_available[axis_idx] {
+            results[axis_idx].na_reason = Some("No gyroUnfilt data");
+            continue;
+        }
+        if !d_term_available[axis_idx] {
+            // na_reason already "No D-term data" from init
             continue;
         }
 
@@ -141,6 +165,7 @@ pub fn calculate_d_term_filtering_delay_comparison(
                 gyro_unfilt_data.len(),
                 d_term_filtered_data.len()
             );
+            results[axis_idx].na_reason = Some("Insufficient samples");
             continue;
         }
 
@@ -157,6 +182,7 @@ pub fn calculate_d_term_filtering_delay_comparison(
                 "  {}: D-term appears disabled (max abs value: {:.2e}, likely D gain = 0)",
                 axis_name, d_term_max_abs
             );
+            results[axis_idx].na_reason = Some("D gain disabled");
             continue;
         }
 
@@ -175,6 +201,7 @@ pub fn calculate_d_term_filtering_delay_comparison(
                 "  {}: D-term has no variation (std dev: {:.2e}, likely D gain = 0)",
                 axis_name, d_term_std_dev
             );
+            results[axis_idx].na_reason = Some("D gain disabled");
             continue;
         }
 
@@ -203,12 +230,16 @@ pub fn calculate_d_term_filtering_delay_comparison(
                     result.delay_ms,
                     result.confidence * 100.0
                 );
-                results[axis_idx] = Some(result);
+                results[axis_idx] = DTermAxisDelay {
+                    result: Some(result),
+                    na_reason: None,
+                };
             } else {
                 println!(
                     "  {}: D-term delay calculation failed - correlation below D-term threshold",
                     axis_name
                 );
+                results[axis_idx].na_reason = Some("Low signal correlation");
             }
         } else if let Some(result) = calculate_d_term_filtering_delay_enhanced_xcorr(
             &Array1::from_vec(d_term_filtered_data),
@@ -221,12 +252,16 @@ pub fn calculate_d_term_filtering_delay_comparison(
                 result.delay_ms,
                 result.confidence * 100.0
             );
-            results[axis_idx] = Some(result);
+            results[axis_idx] = DTermAxisDelay {
+                result: Some(result),
+                na_reason: None,
+            };
         } else {
             println!(
                 "  {}: D-term delay calculation failed - correlation below D-term threshold",
                 axis_name
             );
+            results[axis_idx].na_reason = Some("Low signal correlation");
         }
     }
 
@@ -235,7 +270,7 @@ pub fn calculate_d_term_filtering_delay_comparison(
         .iter()
         .enumerate()
         .filter_map(|(idx, result)| {
-            if result.is_some() {
+            if result.result.is_some() {
                 Some(AXIS_NAMES[idx])
             } else {
                 None
diff --git a/src/data_analysis/filter_delay.rs b/src/data_analysis/filter_delay.rs
index 658e363..ff6778d 100644
--- a/src/data_analysis/filter_delay.rs
+++ b/src/data_analysis/filter_delay.rs
@@ -188,6 +188,7 @@ pub fn calculate_average_filtering_delay_comparison(
     sample_rate: f64,
 ) -> DelayAnalysisResult {
     let mut all_results: Vec<DelayResult> = Vec::new();
+    let mut axis_delays: Vec<Option<(f32, f32)>> = vec![None; AXIS_NAMES.len()];
 
     // First, diagnose data availability
     println!("=== Gyro Data Availability Diagnostic ===");
@@ -274,6 +275,13 @@ pub fn calculate_average_filtering_delay_comparison(
                 }
             }
 
+            // Capture per-axis delay for callers that need it (e.g. report)
+            if let Some(r) = axis_results
+                .iter()
+                .find(|r| r.method == "Enhanced Cross-Correlation")
+            {
+                axis_delays[axis] = Some((r.delay_ms, r.confidence));
+            }
             all_results.extend(axis_results);
         }
     }
@@ -298,17 +306,20 @@ pub fn calculate_average_filtering_delay_comparison(
             DelayAnalysisResult {
                 average_delay: Some(avg_delay),
                 results: method_summaries,
+                axis_delays,
             }
         } else {
             DelayAnalysisResult {
                 average_delay: None,
                 results: method_summaries,
+                axis_delays,
             }
         }
     } else {
         DelayAnalysisResult {
             average_delay: None,
             results: Vec::new(),
+            axis_delays,
         }
     }
 }
@@ -325,6 +336,8 @@ pub struct DelayResult {
 pub struct DelayAnalysisResult {
     pub average_delay: Option<f32>,
     pub results: Vec<DelayResult>,
+    /// Per-axis (delay_ms, confidence 0–1); index matches AXIS_NAMES order
+    pub axis_delays: Vec<Option<(f32, f32)>>,
 }
 
 /// Calculate filtering delay using enhanced cross-correlation method only
diff --git a/src/lib.rs b/src/lib.rs
index 6d46635..8cfabd7 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -10,4 +10,5 @@ pub mod font_config;
 pub mod pid_context;
 pub mod plot_framework;
 pub mod plot_functions;
+pub mod report;
 pub mod types;
diff --git a/src/main.rs b/src/main.rs
index 6339372..94db8f2 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -9,6 +9,7 @@ mod font_config;
 mod pid_context;
 mod plot_framework;
 mod plot_functions;
+mod report;
 mod types;
 
 use std::collections::{BTreeMap, HashSet};
@@ -19,6 +20,7 @@ use std::path::{Path, PathBuf};
 
 use ndarray::Array1;
 
+use crate::axis_names::AXIS_COUNT;
 use crate::data_analysis::torque_inertia_profiler::{extract_punch_ratios, AircraftProfile};
 use crate::types::StepResponseResults;
 
@@ -182,6 +184,7 @@ use crate::pid_context::PidContext;
 // Data analysis imports
 use crate::data_analysis::calc_step_response;
 use crate::data_analysis::calc_step_response::{compute_setpoint_authority, SetpointAuthority};
+use crate::data_analysis::filter_response;
 
 /// Expand input paths to a list of CSV files.
 /// If a path is a file, validate CSV extension before adding.
@@ -647,6 +650,12 @@ fn process_file(
     println!("Note: Optimal P:D ratio varies per aircraft. Check step response for overshoot/undershoot.");
     println!();
 
+    let pd_ratios_for_report: [Option<f64>; AXIS_COUNT] = [
+        pid_metadata.roll.calculate_pd_ratio(),
+        pid_metadata.pitch.calculate_pd_ratio(),
+        pid_metadata.yaw.calculate_pd_ratio(),
+    ];
+
     let mut has_nonzero_f_term_data = [false; 3];
     for axis in 0..crate::axis_names::AXIS_NAMES.len() {
         if f_term_header_found[axis]
@@ -755,6 +764,14 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
     let mut recommended_d_min_aggressive: [Option<u32>; 3] = [None, None, None];
     let mut recommended_d_max_aggressive: [Option<u32>; 3] = [None, None, None];
 
+    // Setpoint authority per axis (captured for report)
+    let mut setpoint_authority_names: [Option<&'static str>; AXIS_COUNT] =
+        std::array::from_fn(|_| None);
+    let mut setpoint_authority_means: [Option<f32>; AXIS_COUNT] = std::array::from_fn(|_| None);
+
+    // Step response warnings per axis (captured for report)
+    let mut step_warnings: [Vec<String>; AXIS_COUNT] = std::array::from_fn(|_| Vec::new());
+
     if let Some(sr) = sample_rate {
         for axis_index in 0..crate::axis_names::AXIS_NAMES.len() {
             let axis_name = crate::axis_names::AXIS_NAMES[axis_index];
@@ -812,7 +829,7 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
         println!("      Always test in a safe environment. Conservative = safer first step.");
         println!("      Moderate = for experienced pilots (test carefully to avoid hot motors).");
         println!();
-        for axis_index in 0..2 {
+        for axis_index in 0..crate::axis_names::ROLL_PITCH_AXIS_COUNT {
             // Only Roll (0) and Pitch (1)
             let axis_name = crate::axis_names::AXIS_NAMES[axis_index];
 
@@ -977,6 +994,8 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
                                         valid_window_max_setpoints.as_slice().unwrap_or(&[]),
                                     )
                                     .unwrap_or((SetpointAuthority::Low, 0.0));
+                                    setpoint_authority_names[axis_index] = Some(authority.name());
+                                    setpoint_authority_means[axis_index] = Some(authority_mean);
                                     println!(
                                         "  Setpoint Authority: {} (mean={:.0}dps \u{22a2}\u{2265}{}dps)",
                                         authority.name(),
@@ -996,6 +1015,9 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
                                             println!(
                                                 "      - Check for bent props, loose hardware, or damaged motors"
                                             );
+                                            step_warnings[axis_index].push(format!(
+                                                "Severe overshoot (Peak={peak_value:.2}): P may be too high, or check for bent props/loose hardware/damaged motors"
+                                            ));
                                         }
 
                                         // Check for unreasonable P:D ratios
@@ -1006,11 +1028,17 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
                                             println!(
                                                 "      Consider increasing P instead of only adding D"
                                             );
+                                            step_warnings[axis_index].push(format!(
+                                                "Recommended P:D ratio ({rec_ratio:.2}) is very low — consider increasing P instead of only adding D"
+                                            ));
                                         } else if rec_ratio
                                             > crate::constants::MAX_REASONABLE_PD_RATIO
                                         {
                                             println!("  ⚠️  WARNING: Recommended P:D ratio ({rec_ratio:.2}) is very high");
                                             println!("      Consider decreasing P or checking for overdamped response");
+                                            step_warnings[axis_index].push(format!(
+                                                "Recommended P:D ratio ({rec_ratio:.2}) is very high — consider decreasing P or checking for overdamped response"
+                                            ));
                                         }
 
                                         // Show conservative recommendation
@@ -1152,6 +1180,66 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
         println!();
     }
 
+    // Collect step response analysis into typed report structs
+    let mut step_reports: Vec<report::StepAxisReport> = Vec::new();
+    {
+        let dmax_enabled = pid_metadata.is_dmax_enabled();
+        for axis_index in 0..crate::axis_names::ROLL_PITCH_AXIS_COUNT {
+            if let (Some(peak_value), Some(current_pd_ratio), Some(assessment)) = (
+                peak_values[axis_index],
+                current_pd_ratios[axis_index],
+                assessments[axis_index],
+            ) {
+                let conservative =
+                    recommended_pd_conservative[axis_index].map(|pd| report::DTermRec {
+                        pd_ratio: pd,
+                        d: recommended_d_conservative[axis_index],
+                        d_min: recommended_d_min_conservative[axis_index],
+                        d_max: recommended_d_max_conservative[axis_index],
+                    });
+                let moderate = recommended_pd_aggressive[axis_index].map(|pd| report::DTermRec {
+                    pd_ratio: pd,
+                    d: recommended_d_aggressive[axis_index],
+                    d_min: recommended_d_min_aggressive[axis_index],
+                    d_max: recommended_d_max_aggressive[axis_index],
+                });
+                let aggressive = if assessment == "Significant overshoot" {
+                    let aggressive_pd =
+                        current_pd_ratio * crate::constants::PD_RATIO_AGGRESSIVE_MULTIPLIER;
+                    let (rec_d, rec_d_min, rec_d_max) = if axis_index == 0 {
+                        pid_metadata
+                            .roll
+                            .calculate_goal_d_with_range(aggressive_pd, dmax_enabled)
+                    } else {
+                        pid_metadata
+                            .pitch
+                            .calculate_goal_d_with_range(aggressive_pd, dmax_enabled)
+                    };
+                    Some(report::DTermRec {
+                        pd_ratio: aggressive_pd,
+                        d: rec_d,
+                        d_min: rec_d_min,
+                        d_max: rec_d_max,
+                    })
+                } else {
+                    None
+                };
+                step_reports.push(report::StepAxisReport {
+                    axis_name: crate::axis_names::AXIS_NAMES[axis_index],
+                    peak_value,
+                    assessment,
+                    current_pd_ratio,
+                    conservative,
+                    moderate,
+                    aggressive,
+                    setpoint_authority_name: setpoint_authority_names[axis_index],
+                    setpoint_authority_mean: setpoint_authority_means[axis_index],
+                    warnings: std::mem::take(&mut step_warnings[axis_index]),
+                });
+            }
+        }
+    }
+
     // Optimal P Estimation Analysis (if enabled)
     // Store results for both console output and PNG overlay
     let mut optimal_p_analyses: [Option<
@@ -1305,6 +1393,8 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
         }
     }
 
+    let optimal_p_for_report = optimal_p_analyses.clone();
+
     // Create RAII guard BEFORE changing directory if needed
     let _cwd_guard = if let Some(output_dir) = output_dir {
         // Create guard to save current directory BEFORE changing it
@@ -1417,8 +1507,8 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
         )?;
     }
 
-    if plot_config.gyro_spectrums {
-        plot_gyro_spectrums(
+    let gyro_analysis = if plot_config.gyro_spectrums {
+        Some(plot_gyro_spectrums(
             &all_log_data,
             &root_name_string,
             sample_rate,
@@ -1426,8 +1516,10 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
             analysis_opts.show_butterworth,
             using_debug_fallback,
             debug_mode_label,
-        )?;
-    }
+        )?)
+    } else {
+        None
+    };
 
     if plot_config.d_term_psd {
         plot_d_term_psd(
@@ -1441,7 +1533,7 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
         )?;
     }
 
-    if plot_config.d_term_spectrums {
+    let dterm_results = if plot_config.d_term_spectrums {
         plot_d_term_spectrums(
             &all_log_data,
             &root_name_string,
@@ -1450,12 +1542,16 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
             analysis_opts.show_butterworth,
             using_debug_fallback,
             debug_mode_label,
-        )?;
-    }
+        )?
+    } else {
+        vec![]
+    };
 
-    if plot_config.motor_spectrums {
-        plot_motor_spectrums(&all_log_data, &root_name_string, sample_rate)?;
-    }
+    let motor_results = if plot_config.motor_spectrums {
+        plot_motor_spectrums(&all_log_data, &root_name_string, sample_rate)?
+    } else {
+        vec![]
+    };
 
     if plot_config.psd {
         plot_psd(
@@ -1467,7 +1563,7 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
         )?;
     }
 
-    if plot_config.bode {
+    let bode_results = if plot_config.bode {
         eprintln!();
         eprintln!("⚠️  WARNING: Bode plots are designed for controlled test flights with system-identification inputs.");
         eprintln!(
@@ -1479,8 +1575,10 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
             &root_name_string,
             sample_rate,
             analysis_opts.debug_mode,
-        )?;
-    }
+        )?
+    } else {
+        vec![]
+    };
 
     if plot_config.psd_db_heatmap {
         plot_psd_db_heatmap(
@@ -1510,7 +1608,107 @@ INFO ({input_file_str}): Skipping Step Response input data filtering: {reason}."
         plot_pid_activity(&all_log_data, &root_name_string, Some(&header_metadata))?;
     }
 
+    // --- Filter configuration (from header metadata, independent of CSV data) ---
+    let filter_config = Some(filter_response::parse_filter_config(&header_metadata));
+    let dynamic_notch = filter_response::extract_dynamic_notch_range(Some(&header_metadata));
+    let rpm_filter = filter_response::extract_rpm_filter_config(Some(&header_metadata));
+
+    // --- Collect generated PNG filenames ---
+    let mut png_links: Vec<String> = Vec::new();
+
+    if plot_config.step_response {
+        // Step response filename includes duration and optional dps suffix — scan for it.
+        let prefix = format!("{root_name_string}_Step_Response_stacked_plot_");
+        if let Ok(entries) = std::fs::read_dir(".") {
+            let mut matches: Vec<String> = entries
+                .flatten()
+                .map(|e| e.file_name().to_string_lossy().into_owned())
+                .filter(|n| n.starts_with(&prefix) && n.ends_with(".png"))
+                .collect();
+            matches.sort();
+            png_links.extend(matches);
+        }
+    }
+    if plot_config.pidsum_error_setpoint {
+        png_links.push(format!(
+            "{root_name_string}_PIDsum_PIDerror_Setpoint_stacked.png"
+        ));
+    }
+    if plot_config.setpoint_vs_gyro {
+        png_links.push(format!("{root_name_string}_SetpointVsGyro_stacked.png"));
+    }
+    if plot_config.setpoint_derivative {
+        png_links.push(format!("{root_name_string}_SetpointDerivative_stacked.png"));
+    }
+    if plot_config.gyro_vs_unfilt {
+        png_links.push(format!("{root_name_string}_GyroVsUnfilt_stacked.png"));
+    }
+    if plot_config.gyro_spectrums {
+        png_links.push(format!("{root_name_string}_Gyro_Spectrums_comparative.png"));
+    }
+    if plot_config.d_term_psd {
+        png_links.push(format!("{root_name_string}_D_Term_PSD_comparative.png"));
+    }
+    if plot_config.d_term_spectrums {
+        png_links.push(format!(
+            "{root_name_string}_D_Term_Spectrums_comparative.png"
+        ));
+    }
+    if plot_config.motor_spectrums {
+        png_links.push(format!("{root_name_string}_Motor_Spectrums_stacked.png"));
+    }
+    if plot_config.psd {
+        png_links.push(format!("{root_name_string}_Gyro_PSD_comparative.png"));
+    }
+    if plot_config.psd_db_heatmap {
+        png_links.push(format!(
+            "{root_name_string}_Gyro_PSD_Spectrogram_comparative.png"
+        ));
+    }
+    if plot_config.throttle_freq_heatmap {
+        png_links.push(format!(
+            "{root_name_string}_Throttle_Freq_Heatmap_comparative.png"
+        ));
+    }
+    if plot_config.d_term_heatmap {
+        png_links.push(format!("{root_name_string}_D_Term_Heatmap_comparative.png"));
+    }
+    if plot_config.bode {
+        png_links.push(format!("{root_name_string}_Bode_Analysis.png"));
+    }
+    if plot_config.pid_activity {
+        png_links.push(format!("{root_name_string}_PID_Activity_stacked.png"));
+    }
+
+    // --- Markdown Report ---
+    // Must run after all plots so png_links is complete.
+    let report_filename = format!("{root_name_string}_report.md");
+    let report_path = std::path::Path::new(&report_filename);
+    println!("\n--- Generating Report: {report_filename} ---");
+    let flight_report = report::FlightReport {
+        root_name: root_name_string.clone(),
+        sample_rate,
+        header_metadata,
+        pd_ratios: pd_ratios_for_report,
+        step_reports,
+        optimal_p: optimal_p_for_report,
+        gyro_analysis,
+        dterm_results,
+        bode_results,
+        motor_results,
+        png_links,
+        filter_config,
+        dynamic_notch,
+        rpm_filter,
+        debug_fallback: using_debug_fallback,
+        debug_mode_name: debug_mode_label,
+    };
+    report::generate_markdown_report(&flight_report, report_path)
+        .map_err(|e| format!("Report generation failed: {e}"))?;
+    println!("  [OK] Report written.");
+
     // CWD restoration happens automatically when _cwd_guard goes out of scope
+
     println!("--- Finished processing file: {input_file_str} ---");
     Ok(())
 }
diff --git a/src/plot_functions/plot_bode.rs b/src/plot_functions/plot_bode.rs
index 0cc0321..ba46954 100644
--- a/src/plot_functions/plot_bode.rs
+++ b/src/plot_functions/plot_bode.rs
@@ -28,6 +28,12 @@ use crate::font_config::{
 /// Minimum coherence threshold for filtering Bode plot data
 const MIN_COHERENCE_FOR_PLOT: f64 = 0.1;
 
+/// Per-axis result from Bode analysis containing stability margins
+pub struct BodeAxisResult {
+    pub axis_name: String,
+    pub margins: StabilityMargins,
+}
+
 /// Plot Bode analysis for all three axes (Roll, Pitch, Yaw)
 ///
 /// Generates a single 3×3 grid plot with magnitude, phase, and coherence for all axes
@@ -36,12 +42,12 @@ pub fn plot_bode_analysis(
     root_name: &str,
     sample_rate: Option<f64>,
     debug_mode: bool,
-) -> Result<(), Box<dyn Error>> {
+) -> Result<Vec<BodeAxisResult>, Box<dyn Error>> {
     let sr_value = if let Some(sr) = sample_rate {
         sr
     } else {
         println!("\nINFO: Skipping Bode Plot: Sample rate could not be determined.");
-        return Ok(());
+        return Ok(vec![]);
     };
 
     // Estimate transfer functions for all three axes
@@ -154,7 +160,7 @@ pub fn plot_bode_analysis(
     // Early exit if no valid axes
     if tf_results.is_empty() {
         println!("\nINFO: Skipping Bode Plot: No valid transfer function data for any axis.");
-        return Ok(());
+        return Ok(vec![]);
     }
 
     // Generate single combined plot filename
@@ -180,7 +186,14 @@ pub fn plot_bode_analysis(
         }
     }
 
-    Ok(())
+    Ok(tf_results
+        .iter()
+        .enumerate()
+        .map(|(i, tf)| BodeAxisResult {
+            axis_name: tf.axis_name.clone(),
+            margins: margins_results[i].clone(),
+        })
+        .collect())
 }
 
 /// Create a grid Bode plot (1 to 3 axes × 3 plot types)
diff --git a/src/plot_functions/plot_d_term_psd.rs b/src/plot_functions/plot_d_term_psd.rs
index 8d4b68c..034f071 100644
--- a/src/plot_functions/plot_d_term_psd.rs
+++ b/src/plot_functions/plot_d_term_psd.rs
@@ -59,7 +59,7 @@ pub fn plot_d_term_psd(
         d_term_delay::calculate_d_term_filtering_delay_comparison(log_data, sr_value);
 
     // Check if any delay calculations succeeded - if not, don't show delay in legends
-    let any_delay_calculated = delay_by_axis.iter().any(|result| result.is_some());
+    let any_delay_calculated = delay_by_axis.iter().any(|result| result.result.is_some());
 
     let mut global_max_y_unfilt = f64::NEG_INFINITY;
     let mut global_max_y_filt = f64::NEG_INFINITY;
@@ -329,7 +329,7 @@ pub fn plot_d_term_psd(
 
         // Get delay string for this axis for legend display
         let delay_str = if any_delay_calculated {
-            if let Some(result) = delay_by_axis.get(axis_idx).and_then(|r| r.as_ref()) {
+            if let Some(result) = delay_by_axis.get(axis_idx).and_then(|d| d.result.as_ref()) {
                 format!(
                     "Delay: {:.1}ms(c:{:.0}%)",
                     result.delay_ms,
diff --git a/src/plot_functions/plot_d_term_spectrums.rs b/src/plot_functions/plot_d_term_spectrums.rs
index 0390937..05d6cc2 100644
--- a/src/plot_functions/plot_d_term_spectrums.rs
+++ b/src/plot_functions/plot_d_term_spectrums.rs
@@ -21,6 +21,16 @@ use crate::plot_framework::{
 use crate::plot_functions::peak_detection::find_and_sort_peaks_with_threshold;
 use plotters::style::RGBColor;
 
+/// Per-axis D-term spectrum analysis: peaks (sorted by amplitude) and filtering delay
+pub struct DTermAxisResult {
+    pub axis_name: &'static str,
+    pub peaks: Vec<(f64, f64)>, // (freq_hz, amplitude); [0] = primary, rest = subordinates
+    pub delay_ms: Option<f32>,
+    pub delay_confidence: Option<f32>, // 0.0–1.0
+    /// Why delay is N/A (None when delay_ms is Some)
+    pub na_reason: Option<&'static str>,
+}
+
 /// Generates a stacked plot with two columns per axis, showing Unfiltered D-term and Filtered D-term spectrums (linear amplitude).
 /// Now includes filter response curve overlays based on header metadata.
 /// Unfiltered D-term is calculated as the derivative of gyroUnfilt.
@@ -33,12 +43,12 @@ pub fn plot_d_term_spectrums(
     show_butterworth: bool,
     using_debug_fallback: bool,
     debug_mode_name: Option<&str>,
-) -> Result<(), Box<dyn Error>> {
+) -> Result<Vec<DTermAxisResult>, Box<dyn Error>> {
     // Clone debug mode name to move into closures
     let debug_mode_name_owned = debug_mode_name.map(|s| s.to_string());
     // Input validation
     if log_data.is_empty() {
-        return Ok(()); // No data to process
+        return Ok(vec![]); // No data to process
     }
 
     if root_name.is_empty() {
@@ -52,11 +62,11 @@ pub fn plot_d_term_spectrums(
             sr
         } else {
             println!("\nINFO: Skipping D-Term Spectrum Plot: Invalid sample rate provided.");
-            return Ok(());
+            return Ok(vec![]);
         }
     } else {
         println!("\nINFO: Skipping D-Term Spectrum Plot: Sample rate could not be determined.");
-        return Ok(());
+        return Ok(vec![]);
     };
 
     // Calculate filtering delay using enhanced cross-correlation on D-terms
@@ -80,6 +90,8 @@ pub fn plot_d_term_spectrums(
 
     // Store axis spectrum data
     let mut axis_spectrums: Vec<AxisSpectrum> = Vec::new();
+    // Per-axis analysis results for the report
+    let mut dterm_results: Vec<DTermAxisResult> = Vec::new();
 
     // Iterate safely over the minimum of AXIS_NAMES.len() and the fixed array size
     let axis_count = AXIS_NAMES.len().min(3);
@@ -263,8 +275,20 @@ pub fn plot_d_term_spectrums(
         let _filt_primary_peak = filt_primary_peak;
         let filt_peaks = Vec::new();
 
+        // Capture per-axis data for report before peaks are moved into plot configs
+        let axis_delay = delay_by_axis.get(axis_idx);
+        let delay_result = axis_delay.and_then(|d| d.result.as_ref());
+        let na_reason = axis_delay.and_then(|d| d.na_reason);
+        dterm_results.push(DTermAxisResult {
+            axis_name,
+            peaks: unfilt_peaks.clone(),
+            delay_ms: delay_result.map(|r| r.delay_ms),
+            delay_confidence: delay_result.map(|r| r.confidence),
+            na_reason,
+        });
+
         // Get delay string for this axis for legend display
-        let delay_str = if let Some(result) = delay_by_axis.get(axis_idx).and_then(|r| r.as_ref()) {
+        let delay_str = if let Some(result) = delay_result {
             format!(
                 "Delay: {:.1}ms(c:{:.0}%)",
                 result.delay_ms,
@@ -541,7 +565,7 @@ pub fn plot_d_term_spectrums(
 
     if overall_max_y_amplitude <= 0.0 {
         println!("  No valid D-term spectrum data found. Skipping D-term spectrum plot.");
-        return Ok(());
+        return Ok(dterm_results);
     }
 
     draw_dual_spectrum_plot(
@@ -558,7 +582,7 @@ pub fn plot_d_term_spectrums(
     )?;
 
     println!("  D-term spectrum plot saved as '{}'", output_file);
-    Ok(())
+    Ok(dterm_results)
 }
 
 // Removed duplicate function: calculate_d_term_filtering_delay_comparison
diff --git a/src/plot_functions/plot_gyro_spectrums.rs b/src/plot_functions/plot_gyro_spectrums.rs
index 6c50ee6..4968ebe 100644
--- a/src/plot_functions/plot_gyro_spectrums.rs
+++ b/src/plot_functions/plot_gyro_spectrums.rs
@@ -22,6 +22,19 @@ use crate::plot_functions::peak_detection::find_and_sort_peaks_with_threshold;
 use crate::types::AllFFTData;
 use plotters::style::RGBColor;
 
+/// Per-axis spectrum peaks and filtering delay (unfiltered gyro)
+pub struct GyroSpectrumAxisResult {
+    pub axis_name: &'static str,
+    pub peaks: Vec<(f64, f64)>, // (freq_hz, amplitude); [0] = primary, rest = subordinates
+    pub delay_ms: Option<f32>,
+    pub delay_confidence: Option<f32>, // 0.0–1.0
+}
+
+/// Gyro spectrum analysis: per-axis peaks and filtering delay
+pub struct GyroAnalysisResult {
+    pub axes: Vec<GyroSpectrumAxisResult>,
+}
+
 /// Generates a stacked plot with two columns per axis, showing Unfiltered and Filtered Gyro spectrums.
 /// Now includes filter response curve overlays based on header metadata.
 pub fn plot_gyro_spectrums(
@@ -32,7 +45,7 @@ pub fn plot_gyro_spectrums(
     show_butterworth: bool,
     using_debug_fallback: bool,
     debug_mode_name: Option<&str>,
-) -> Result<(), Box<dyn Error>> {
+) -> Result<GyroAnalysisResult, Box<dyn Error>> {
     // Clone debug mode name to move into closures
     let debug_mode_name_owned = debug_mode_name.map(|s| s.to_string());
     let output_file = format!("{root_name}_Gyro_Spectrums_comparative.png");
@@ -42,12 +55,13 @@ pub fn plot_gyro_spectrums(
         sr
     } else {
         println!("\nINFO: Skipping Gyro Spectrum Plot: Sample rate could not be determined.");
-        return Ok(());
+        return Ok(GyroAnalysisResult { axes: vec![] });
     };
 
     // Calculate filtering delay using enhanced cross-correlation
     let delay_analysis =
         filter_delay::calculate_average_filtering_delay_comparison(log_data, sr_value);
+    let axis_delays = delay_analysis.axis_delays.clone();
     let delay_comparison_results = if !delay_analysis.results.is_empty() {
         Some(delay_analysis.results)
     } else {
@@ -207,6 +221,25 @@ pub fn plot_gyro_spectrums(
         overall_max_y_amplitude = SPECTRUM_Y_AXIS_FLOOR;
     }
 
+    // Extract per-axis peaks and delay before the closure takes ownership of all_fft_raw_data
+    let gyro_axes: Vec<GyroSpectrumAxisResult> = (0..axis_count)
+        .map(|axis_idx| {
+            let peaks = all_fft_raw_data[axis_idx]
+                .as_ref()
+                .map_or(vec![], |(_, unfilt_peaks, _, _)| unfilt_peaks.clone());
+            let (delay_ms, delay_confidence) = axis_delays
+                .get(axis_idx)
+                .and_then(|d| *d)
+                .map_or((None, None), |(d, c)| (Some(d), Some(c)));
+            GyroSpectrumAxisResult {
+                axis_name: AXIS_NAMES[axis_idx],
+                peaks,
+                delay_ms,
+                delay_confidence,
+            }
+        })
+        .collect();
+
     draw_dual_spectrum_plot(&output_file, root_name, plot_type_name, move |axis_index| {
         if let Some((unfilt_series_data, unfilt_peaks, filt_series_data, filt_peaks)) =
             all_fft_raw_data[axis_index].as_ref().cloned()
@@ -578,7 +611,9 @@ pub fn plot_gyro_spectrums(
                 filtered: None,
             })
         }
-    })
+    })?;
+
+    Ok(GyroAnalysisResult { axes: gyro_axes })
 }
 
 // src/plot_functions/plot_gyro_spectrums.rs
diff --git a/src/plot_functions/plot_gyro_vs_unfilt.rs b/src/plot_functions/plot_gyro_vs_unfilt.rs
index 63f7c10..8d90ff3 100644
--- a/src/plot_functions/plot_gyro_vs_unfilt.rs
+++ b/src/plot_functions/plot_gyro_vs_unfilt.rs
@@ -35,6 +35,7 @@ pub fn plot_gyro_vs_unfilt(
         DelayAnalysisResult {
             average_delay: None,
             results: Vec::new(),
+            axis_delays: Vec::new(),
         }
     };
 
diff --git a/src/plot_functions/plot_motor_spectrums.rs b/src/plot_functions/plot_motor_spectrums.rs
index 875b6d0..f61043b 100644
--- a/src/plot_functions/plot_motor_spectrums.rs
+++ b/src/plot_functions/plot_motor_spectrums.rs
@@ -30,20 +30,29 @@ const MOTOR_COLORS: [RGBColor; 8] = [
 /// Type alias for motor spectrum data: (frequencies, amplitudes, max_amplitude)
 type MotorSpectrumData = (Vec<f64>, Vec<f64>, f64);
 
+/// Per-motor result from oscillation analysis in the MOTOR_OSCILLATION_FREQ range
+pub struct MotorOscillationResult {
+    pub motor_idx: usize,
+    pub max_amplitude: Option<f64>,
+    pub oscillation_detected: bool,
+    pub peak_in_range: Option<f64>,
+    pub avg_in_range: Option<f64>,
+}
+
 /// Generates stacked motor spectrum plots showing frequency content of each motor output.
 /// Useful for identifying motor oscillations, ESC noise, and saturation issues.
 pub fn plot_motor_spectrums(
     log_data: &[LogRowData],
     root_name: &str,
     sample_rate: Option<f64>,
-) -> Result<(), Box<dyn Error>> {
+) -> Result<Vec<MotorOscillationResult>, Box<dyn Error>> {
     let output_file = format!("{root_name}_Motor_Spectrums_stacked.png");
 
     let sr_value = if let Some(sr) = sample_rate {
         sr
     } else {
         println!("\nINFO: Skipping Motor Spectrum Plot: Sample rate could not be determined.");
-        return Ok(());
+        return Ok(vec![]);
     };
 
     // Determine motor count from first row
@@ -51,7 +60,7 @@ pub fn plot_motor_spectrums(
 
     if motor_count == 0 {
         println!("\nINFO: Skipping Motor Spectrum Plot: No motor data available.");
-        return Ok(());
+        return Ok(vec![]);
     }
 
     println!(
@@ -134,9 +143,15 @@ pub fn plot_motor_spectrums(
     }
 
     // Check for oscillation issues (peaks > 3× average in 50-200 Hz range)
+    let mut motor_osc_results: Vec<MotorOscillationResult> = Vec::new();
     for (motor_idx, spectrum_data) in motor_spectrums.iter().enumerate() {
+        let max_amplitude = spectrum_data.as_ref().map(|(_, _, max)| *max);
+        let mut oscillation_detected = false;
+        let mut peak_in_range: Option<f64> = None;
+        let mut avg_in_range: Option<f64> = None;
+
         if let Some((frequencies, amplitudes, _)) = spectrum_data {
-            let freq_range_50_200: Vec<(f64, f64)> = frequencies
+            let freq_range: Vec<(f64, f64)> = frequencies
                 .iter()
                 .zip(amplitudes.iter())
                 .filter(|(f, _)| {
@@ -145,24 +160,32 @@ pub fn plot_motor_spectrums(
                 .map(|(f, a)| (*f, *a))
                 .collect();
 
-            if !freq_range_50_200.is_empty() {
-                let avg_amplitude: f64 = freq_range_50_200.iter().map(|(_, a)| a).sum::<f64>()
-                    / freq_range_50_200.len() as f64;
-                let max_in_range = freq_range_50_200
-                    .iter()
-                    .map(|(_, a)| *a)
-                    .fold(0.0f64, f64::max);
+            if !freq_range.is_empty() {
+                let avg: f64 =
+                    freq_range.iter().map(|(_, a)| a).sum::<f64>() / freq_range.len() as f64;
+                let max = freq_range.iter().map(|(_, a)| *a).fold(0.0f64, f64::max);
+                peak_in_range = Some(max);
+                avg_in_range = Some(avg);
 
-                if max_in_range > MOTOR_OSCILLATION_THRESHOLD_MULTIPLIER * avg_amplitude
-                    && max_in_range > MOTOR_OSCILLATION_ABSOLUTE_THRESHOLD
+                if max > MOTOR_OSCILLATION_THRESHOLD_MULTIPLIER * avg
+                    && max > MOTOR_OSCILLATION_ABSOLUTE_THRESHOLD
                 {
+                    oscillation_detected = true;
                     println!(
                         "  ⚠ Motor {}: Potential oscillation detected in {:.0}-{:.0} Hz range (peak {:.1} >> avg {:.1})",
-                        motor_idx, MOTOR_OSCILLATION_FREQ_MIN_HZ, MOTOR_OSCILLATION_FREQ_MAX_HZ, max_in_range, avg_amplitude
+                        motor_idx, MOTOR_OSCILLATION_FREQ_MIN_HZ, MOTOR_OSCILLATION_FREQ_MAX_HZ, max, avg
                     );
                 }
             }
         }
+
+        motor_osc_results.push(MotorOscillationResult {
+            motor_idx,
+            max_amplitude,
+            oscillation_detected,
+            peak_in_range,
+            avg_in_range,
+        });
     }
 
     // Generate stacked plots with dynamic row count for motors
@@ -274,5 +297,5 @@ pub fn plot_motor_spectrums(
         println!("  Skipping '{}': No motor data to plot.", output_file);
     }
 
-    Ok(())
+    Ok(motor_osc_results)
 }
diff --git a/src/plot_functions/plot_setpoint_vs_gyro.rs b/src/plot_functions/plot_setpoint_vs_gyro.rs
index d242f82..9f66537 100644
--- a/src/plot_functions/plot_setpoint_vs_gyro.rs
+++ b/src/plot_functions/plot_setpoint_vs_gyro.rs
@@ -31,6 +31,7 @@ pub fn plot_setpoint_vs_gyro(
         DelayAnalysisResult {
             average_delay: None,
             results: Vec::new(),
+            axis_delays: vec![None; AXIS_NAMES.len()],
         }
     };
 
diff --git a/src/report.rs b/src/report.rs
new file mode 100644
index 0000000..de529b8
--- /dev/null
+++ b/src/report.rs
@@ -0,0 +1,574 @@
+// src/report.rs
+// Structured markdown report capturing computed flight analysis outputs.
+// Derives all content from typed structs returned by analysis functions —
+// no re-reading of CSV data, no duplication of println logic.
+
+use std::error::Error;
+use std::fmt::Write as FmtWrite;
+use std::fs;
+use std::path::Path;
+
+use crate::axis_names::{AXIS_COUNT, AXIS_NAMES};
+use crate::constants::{MOTOR_OSCILLATION_FREQ_MAX_HZ, MOTOR_OSCILLATION_FREQ_MIN_HZ};
+use crate::data_analysis::filter_response::{
+    AllFilterConfigs, DynamicNotchConfig, RpmFilterConfig,
+};
+use crate::data_analysis::optimal_p_estimation::{OptimalPAnalysis, PRecommendation};
+use crate::data_analysis::transfer_function_estimation::Confidence;
+use crate::plot_functions::plot_bode::BodeAxisResult;
+use crate::plot_functions::plot_d_term_spectrums::DTermAxisResult;
+use crate::plot_functions::plot_gyro_spectrums::GyroAnalysisResult;
+use crate::plot_functions::plot_motor_spectrums::MotorOscillationResult;
+
+/// D-term recommendation for one tier (conservative, moderate, or aggressive)
+pub struct DTermRec {
+    pub pd_ratio: f64,
+    pub d: Option<u32>,
+    pub d_min: Option<u32>,
+    pub d_max: Option<u32>,
+}
+
+/// Step response analysis result for one axis
+pub struct StepAxisReport {
+    pub axis_name: &'static str,
+    pub peak_value: f64,
+    pub assessment: &'static str,
+    pub current_pd_ratio: f64,
+    pub conservative: Option<DTermRec>,
+    pub moderate: Option<DTermRec>,
+    pub aggressive: Option<DTermRec>,
+    pub setpoint_authority_name: Option<&'static str>,
+    pub setpoint_authority_mean: Option<f32>,
+    pub warnings: Vec<String>,
+}
+
+/// Aggregated analysis results collected from one flight log
+pub struct FlightReport {
+    pub root_name: String,
+    pub sample_rate: Option<f64>,
+    pub header_metadata: Vec<(String, String)>,
+    pub pd_ratios: [Option<f64>; AXIS_COUNT],
+    pub step_reports: Vec<StepAxisReport>,
+    pub optimal_p: [Option<OptimalPAnalysis>; AXIS_COUNT],
+    pub gyro_analysis: Option<GyroAnalysisResult>,
+    pub dterm_results: Vec<DTermAxisResult>,
+    pub bode_results: Vec<BodeAxisResult>,
+    pub motor_results: Vec<MotorOscillationResult>,
+    pub png_links: Vec<String>,
+    pub filter_config: Option<AllFilterConfigs>,
+    pub dynamic_notch: Option<DynamicNotchConfig>,
+    pub rpm_filter: Option<RpmFilterConfig>,
+    /// True when gyroUnfilt came from debug channels instead of dedicated gyroUnfilt columns.
+    pub debug_fallback: bool,
+    pub debug_mode_name: Option<&'static str>,
+}
+
+/// Generate a structured markdown report and write it to `output_path`.
+pub fn generate_markdown_report(
+    report: &FlightReport,
+    output_path: &Path,
+) -> Result<(), Box<dyn Error>> {
+    let mut md = String::new();
+
+    writeln!(md, "# BlackBox Flight Report: {}", report.root_name)?;
+    writeln!(md)?;
+
+    // --- Metadata ---
+    writeln!(md, "## Metadata")?;
+    writeln!(md)?;
+    match report.sample_rate {
+        Some(sr) => writeln!(md, "- **Sample Rate:** {:.1} Hz", sr)?,
+        None => writeln!(md, "- **Sample Rate:** Unknown")?,
+    }
+    let interesting_keys = [
+        "Firmware revision",
+        "Craft name",
+        "looptime",
+        "gyro_lpf_hz",
+        "dterm_lpf_hz",
+        "pid_process_denom",
+        "rollPID",
+        "pitchPID",
+        "yawPID",
+    ];
+    for (k, v) in &report.header_metadata {
+        if interesting_keys.iter().any(|ik| k.eq_ignore_ascii_case(ik)) {
+            writeln!(md, "- **{}:** {}", k, v)?;
+        }
+    }
+    if report.debug_fallback {
+        let mode_str = report
+            .debug_mode_name
+            .map_or(String::new(), |m| format!(" ({})", m));
+        writeln!(
+            md,
+            "- **⚠ gyroUnfilt source:** debug[0-2] fallback{} — unfiltered gyro spectrum and filtering delay derived from debug channels, not dedicated gyroUnfilt columns",
+            mode_str
+        )?;
+    }
+    writeln!(md)?;
+
+    // --- Filter Configuration ---
+    if let Some(ref fc) = report.filter_config {
+        let axis_labels = ["Roll", "Pitch", "Yaw"];
+        // Only emit section if at least one axis has any filter configured
+        let has_any = (0..3).any(|i| {
+            fc.gyro[i].lpf1.is_some()
+                || fc.gyro[i].lpf2.is_some()
+                || fc.gyro[i].dynamic_lpf1.is_some()
+                || fc.gyro[i].imuf.is_some()
+        });
+        if has_any {
+            writeln!(md, "## Filter Configuration")?;
+            writeln!(md)?;
+            writeln!(md, "| Axis | LPF1 | LPF2 | IMUF / Pseudo-Kalman |")?;
+            writeln!(md, "|------|------|------|----------------------|")?;
+            for (i, label) in axis_labels.iter().enumerate() {
+                let lpf1 =
+                    fmt_filter_stage(fc.gyro[i].lpf1.as_ref(), fc.gyro[i].dynamic_lpf1.as_ref());
+                let lpf2 = fc.gyro[i]
+                    .lpf2
+                    .as_ref()
+                    .map_or("—".into(), fmt_static_filter);
+                let imuf = fc.gyro[i].imuf.as_ref().map_or("—".into(), fmt_imuf);
+                writeln!(md, "| {} | {} | {} | {} |", label, lpf1, lpf2, imuf)?;
+            }
+            writeln!(md)?;
+
+            if let Some(ref dn) = report.dynamic_notch {
+                writeln!(
+                    md,
+                    "- **Dynamic Notch:** {:.0}–{:.0} Hz, Q={:.0}, {} notch{}",
+                    dn.min_hz,
+                    dn.max_hz,
+                    dn.q_factor,
+                    dn.notch_count,
+                    if dn.notch_count == 1 { "" } else { "es" }
+                )?;
+            }
+            if let Some(ref rpm) = report.rpm_filter {
+                writeln!(
+                    md,
+                    "- **RPM Filter:** {} harmonic{}, Q={:.0}, min {:.0} Hz",
+                    rpm.harmonics,
+                    if rpm.harmonics == 1 { "" } else { "s" },
+                    rpm.q_factor,
+                    rpm.min_hz
+                )?;
+            }
+            if report.dynamic_notch.is_some() || report.rpm_filter.is_some() {
+                writeln!(md)?;
+            }
+        }
+    }
+
+    // --- PID Tuning ---
+    writeln!(md, "## PID Tuning")?;
+    writeln!(md)?;
+    writeln!(md, "| Axis | P:D Ratio |")?;
+    writeln!(md, "|------|-----------|")?;
+    let axis_labels = ["Roll", "Pitch", "Yaw (informational)"];
+    for (i, label) in axis_labels.iter().enumerate() {
+        match report.pd_ratios[i] {
+            Some(r) => writeln!(md, "| {} | {:.2} |", label, r)?,
+            None => writeln!(md, "| {} | N/A |", label)?,
+        }
+    }
+    writeln!(md)?;
+
+    // --- Step Response Analysis ---
+    if !report.step_reports.is_empty() {
+        writeln!(md, "## Step Response Analysis")?;
+        writeln!(md)?;
+        for axis_report in &report.step_reports {
+            writeln!(md, "### {}", axis_report.axis_name)?;
+            writeln!(md)?;
+            writeln!(md, "- **Peak Value:** {:.3}", axis_report.peak_value)?;
+            writeln!(md, "- **Assessment:** {}", axis_report.assessment)?;
+            writeln!(md, "- **Current P:D:** {:.2}", axis_report.current_pd_ratio)?;
+            if let (Some(name), Some(mean)) = (
+                axis_report.setpoint_authority_name,
+                axis_report.setpoint_authority_mean,
+            ) {
+                writeln!(
+                    md,
+                    "- **Setpoint Authority:** {} (mean {:.0} dps)",
+                    name, mean
+                )?;
+            }
+            if let Some(rec) = &axis_report.conservative {
+                writeln!(
+                    md,
+                    "- **Recommendation (conservative):** {}",
+                    fmt_dterm_rec(rec)
+                )?;
+            }
+            if let Some(rec) = &axis_report.moderate {
+                writeln!(
+                    md,
+                    "- **Recommendation (moderate):** {}",
+                    fmt_dterm_rec(rec)
+                )?;
+            }
+            if let Some(rec) = &axis_report.aggressive {
+                writeln!(
+                    md,
+                    "- **Recommendation (aggressive):** {}",
+                    fmt_dterm_rec(rec)
+                )?;
+            }
+            if axis_report.conservative.is_none()
+                && axis_report.moderate.is_none()
+                && axis_report.aggressive.is_none()
+            {
+                writeln!(
+                    md,
+                    "- **Recommendation:** No obvious tuning adjustments needed"
+                )?;
+            }
+            for w in &axis_report.warnings {
+                writeln!(md, "- **⚠ Warning:** {}", w)?;
+            }
+            writeln!(md)?;
+        }
+    }
+
+    // --- Gyro Analysis (per-axis filtering delay + spectrum peaks) ---
+    if let Some(gyro) = &report.gyro_analysis {
+        let has_gyro_data = gyro
+            .axes
+            .iter()
+            .any(|a| !a.peaks.is_empty() || a.delay_ms.is_some());
+        if has_gyro_data {
+            writeln!(md, "## Gyro Analysis")?;
+            writeln!(md)?;
+            writeln!(
+                md,
+                "| Axis | Delay (ms) | Confidence | Peak | Freq (Hz) | Amplitude |"
+            )?;
+            writeln!(
+                md,
+                "|------|-----------|------------|------|-----------|-----------|"
+            )?;
+            for axis in &gyro.axes {
+                let delay = axis.delay_ms.map_or("N/A".into(), |v| format!("{:.2}", v));
+                let conf = axis
+                    .delay_confidence
+                    .map_or("N/A".into(), |v| format!("{:.0}%", v * 100.0));
+                if axis.peaks.is_empty() {
+                    writeln!(
+                        md,
+                        "| {} | {} | {} | N/A | N/A | N/A |",
+                        axis.axis_name, delay, conf
+                    )?;
+                } else {
+                    for (idx, (freq, amp)) in axis.peaks.iter().enumerate() {
+                        let label = if idx == 0 {
+                            "Primary".to_string()
+                        } else {
+                            format!("Sub {}", idx)
+                        };
+                        let (d_col, c_col) = if idx == 0 {
+                            (delay.clone(), conf.clone())
+                        } else {
+                            ("".into(), "".into())
+                        };
+                        writeln!(
+                            md,
+                            "| {} | {} | {} | {} | {:.1} | {:.2} |",
+                            axis.axis_name, d_col, c_col, label, freq, amp
+                        )?;
+                    }
+                }
+            }
+            writeln!(md)?;
+        }
+    }
+
+    // --- D-Term Analysis (per-axis delay + spectrum peaks) ---
+    let has_dterm_data = report
+        .dterm_results
+        .iter()
+        .any(|r| !r.peaks.is_empty() || r.delay_ms.is_some());
+    if has_dterm_data {
+        writeln!(md, "## D-Term Analysis")?;
+        writeln!(md)?;
+        writeln!(
+            md,
+            "| Axis | Delay (ms) | Confidence | Peak | Freq (Hz) | Amplitude |"
+        )?;
+        writeln!(
+            md,
+            "|------|-----------|------------|------|-----------|-----------|"
+        )?;
+        for r in &report.dterm_results {
+            let delay = r.delay_ms.map_or_else(
+                || {
+                    r.na_reason
+                        .map_or("N/A".into(), |reason| format!("N/A ({})", reason))
+                },
+                |v| format!("{:.1}", v),
+            );
+            let conf = r
+                .delay_confidence
+                .map_or("N/A".into(), |v| format!("{:.0}%", v * 100.0));
+            if r.peaks.is_empty() {
+                writeln!(
+                    md,
+                    "| {} | {} | {} | N/A | N/A | N/A |",
+                    r.axis_name, delay, conf
+                )?;
+            } else {
+                for (idx, (freq, amp)) in r.peaks.iter().enumerate() {
+                    let label = if idx == 0 {
+                        "Primary".to_string()
+                    } else {
+                        format!("Sub {}", idx)
+                    };
+                    let (d_col, c_col) = if idx == 0 {
+                        (delay.clone(), conf.clone())
+                    } else {
+                        ("".into(), "".into())
+                    };
+                    writeln!(
+                        md,
+                        "| {} | {} | {} | {} | {:.1} | {:.2} |",
+                        r.axis_name, d_col, c_col, label, freq, amp
+                    )?;
+                }
+            }
+        }
+        writeln!(md)?;
+    }
+
+    // --- Optimal P Estimation ---
+    let has_optimal_p = report.optimal_p.iter().any(|o| o.is_some());
+    if has_optimal_p {
+        writeln!(md, "## Optimal P Estimation")?;
+        writeln!(md)?;
+        for (axis_index, opt) in report.optimal_p.iter().enumerate() {
+            if let Some(analysis) = opt {
+                writeln!(md, "### {}", AXIS_NAMES[axis_index])?;
+                writeln!(md)?;
+                writeln!(md, "- **Current P:** {}", analysis.current_p)?;
+                if let Some(d) = analysis.current_d {
+                    writeln!(md, "- **Current D:** {}", d)?;
+                }
+                writeln!(
+                    md,
+                    "- **Td measured:** {:.1} ms ({} samples)",
+                    analysis.td_stats.mean_ms, analysis.td_stats.num_samples
+                )?;
+                writeln!(
+                    md,
+                    "- **Td target:** {:.1} ms ± {:.1} ms",
+                    analysis.td_target_ms, analysis.td_tolerance_ms
+                )?;
+                writeln!(
+                    md,
+                    "- **Td deviation:** {:.1}% ({})",
+                    analysis.td_deviation_percent,
+                    analysis.td_deviation.name()
+                )?;
+                writeln!(md, "- **Noise level:** {}", analysis.noise_level.name())?;
+                let rec_text = match &analysis.recommendation {
+                    PRecommendation::Optimal { reasoning } => {
+                        format!("No change — {}", reasoning)
+                    }
+                    PRecommendation::Increase {
+                        conservative_p,
+                        reasoning,
+                    } => format!("Increase P to {} — {}", conservative_p, reasoning),
+                    PRecommendation::Decrease {
+                        recommended_p,
+                        reasoning,
+                    } => format!("Decrease P to {} — {}", recommended_p, reasoning),
+                    PRecommendation::Investigate { issue } => format!("Investigate — {}", issue),
+                };
+                writeln!(md, "- **Recommendation:** {}", rec_text)?;
+                if analysis.source_files > 1 {
+                    writeln!(
+                        md,
+                        "- **Source:** {} files, {} throttle-punch events",
+                        analysis.source_files, analysis.source_events
+                    )?;
+                } else {
+                    writeln!(
+                        md,
+                        "- **Source:** {} throttle-punch events",
+                        analysis.source_events
+                    )?;
+                }
+                writeln!(md)?;
+            }
+        }
+    }
+
+    // --- Bode Analysis ---
+    if !report.bode_results.is_empty() {
+        writeln!(md, "## Bode Analysis")?;
+        writeln!(md)?;
+        writeln!(
+            md,
+            "> ⚠ Bode analysis is designed for controlled system-identification test flights."
+        )?;
+        writeln!(md)?;
+        writeln!(
+            md,
+            "| Axis | Phase Margin | Gain Margin | Gain Crossover | Bandwidth | Confidence |"
+        )?;
+        writeln!(
+            md,
+            "|------|-------------|-------------|----------------|-----------|------------|"
+        )?;
+        for r in &report.bode_results {
+            let m = &r.margins;
+            let pm = m
+                .phase_margin_deg
+                .map_or("N/A".into(), |v| format!("{:.1}°", v));
+            let gm = m
+                .gain_margin_db
+                .map_or("N/A".into(), |v| format!("{:.1} dB", v));
+            let gc = m
+                .gain_crossover_hz
+                .map_or("N/A".into(), |v| format!("{:.2} Hz", v));
+            let bw = m
+                .bandwidth_hz
+                .map_or("N/A".into(), |v| format!("{:.2} Hz", v));
+            let conf = match m.confidence {
+                Confidence::High => "High",
+                Confidence::Medium => "Medium",
+                Confidence::Low => "Low",
+            };
+            writeln!(
+                md,
+                "| {} | {} | {} | {} | {} | {} |",
+                r.axis_name, pm, gm, gc, bw, conf
+            )?;
+        }
+        writeln!(md)?;
+    }
+
+    // --- Motor Oscillation ---
+    if !report.motor_results.is_empty() {
+        writeln!(md, "## Motor Oscillation")?;
+        writeln!(md)?;
+        writeln!(
+            md,
+            "Analysis range: {:.0}–{:.0} Hz",
+            MOTOR_OSCILLATION_FREQ_MIN_HZ, MOTOR_OSCILLATION_FREQ_MAX_HZ
+        )?;
+        writeln!(md)?;
+        writeln!(
+            md,
+            "| Motor | Max Amplitude | Oscillation | Peak in Range | Avg in Range |"
+        )?;
+        writeln!(
+            md,
+            "|-------|--------------|-------------|---------------|-------------|"
+        )?;
+        for r in &report.motor_results {
+            let max_amp = r
+                .max_amplitude
+                .map_or("N/A".into(), |v| format!("{:.2}", v));
+            let osc = if r.oscillation_detected {
+                "⚠ Detected"
+            } else {
+                "None"
+            };
+            let peak = r
+                .peak_in_range
+                .map_or("N/A".into(), |v| format!("{:.2}", v));
+            let avg = r.avg_in_range.map_or("N/A".into(), |v| format!("{:.2}", v));
+            writeln!(
+                md,
+                "| {} | {} | {} | {} | {} |",
+                r.motor_idx, max_amp, osc, peak, avg
+            )?;
+        }
+        writeln!(md)?;
+    }
+
+    // --- Generated Plots ---
+    if !report.png_links.is_empty() {
+        writeln!(md, "## Generated Plots")?;
+        writeln!(md)?;
+        for name in &report.png_links {
+            writeln!(md, "- [{}]({})", name, name)?;
+        }
+        writeln!(md)?;
+    }
+
+    fs::write(output_path, md)?;
+    Ok(())
+}
+
+fn fmt_static_filter(f: &crate::data_analysis::filter_response::FilterConfig) -> String {
+    match f.q_factor {
+        Some(q) => format!(
+            "{} @ {:.0} Hz (Q={:.2})",
+            f.filter_type.name(),
+            f.cutoff_hz,
+            q
+        ),
+        None => format!("{} @ {:.0} Hz", f.filter_type.name(), f.cutoff_hz),
+    }
+}
+
+fn fmt_filter_stage(
+    lpf1: Option<&crate::data_analysis::filter_response::FilterConfig>,
+    dyn_lpf: Option<&crate::data_analysis::filter_response::DynamicFilterConfig>,
+) -> String {
+    if let Some(f) = lpf1 {
+        fmt_static_filter(f)
+    } else if let Some(d) = dyn_lpf {
+        format!(
+            "Dyn {} {:.0}–{:.0} Hz",
+            d.filter_type.name(),
+            d.min_cutoff_hz,
+            d.max_cutoff_hz
+        )
+    } else {
+        "—".into()
+    }
+}
+
+fn fmt_imuf(f: &crate::data_analysis::filter_response::ImufFilterConfig) -> String {
+    if f.lowpass_cutoff_hz > 0.0 {
+        let stages = match f.ptn_order {
+            1 => "PT1",
+            2 => "2×PT1",
+            3 => "3×PT1",
+            4 => "4×PT1",
+            _ => "PT?",
+        };
+        let rev = f.revision.map_or("IMUF".into(), |r| format!("IMUF v{}", r));
+        let w_part = f
+            .pseudo_kalman_w
+            .map_or(String::new(), |w| format!(", w={:.0}", w));
+        format!(
+            "{}: {} Combined={:.0} Hz per-stage={:.0} Hz (Q={:.1}{})",
+            rev, stages, f.lowpass_cutoff_hz, f.effective_cutoff_hz, f.q_factor, w_part
+        )
+    } else {
+        let w_part = f
+            .pseudo_kalman_w
+            .map_or(String::new(), |w| format!(", w={:.0}", w));
+        format!("Pseudo-Kalman Q={:.1}{}", f.q_factor, w_part)
+    }
+}
+
+fn fmt_dterm_rec(rec: &DTermRec) -> String {
+    if rec.d_min.is_some() || rec.d_max.is_some() {
+        let d_min_s = rec.d_min.map_or("N/A".to_string(), |v| v.to_string());
+        let d_max_s = rec.d_max.map_or("N/A".to_string(), |v| v.to_string());
+        format!(
+            "P:D={:.2} (D-Min≈{}, D-Max≈{})",
+            rec.pd_ratio, d_min_s, d_max_s
+        )
+    } else if let Some(d) = rec.d {
+        format!("P:D={:.2} (D≈{})", rec.pd_ratio, d)
+    } else {
+        format!("P:D={:.2}", rec.pd_ratio)
+    }
+}