Function bodies 301 total

fn compute_min<T: Float>(spectra: &[Vec<T>]) -> Vec<T> {
    if spectra.is_empty() {
        return vec![];
    }

    let n_bins = spectra[0].len();

    (0..n_bins)
        .map(|i| {
            spectra
                .iter()
                .map(|spec| spec[i])
                .fold(T::infinity(), |acc, val| {
                    if val < acc { val } else { acc }
                })
        })
        .collect()
}

fn median_of_slice<T: Float>(values: &mut [T]) -> T {
    values.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let len = values.len();

    if len % 2 == 0 {
        // Even length: average of middle two
        let mid = len / 2;
        (values[mid - 1] + values[mid]) / T::from(2.0).unwrap()
    } else {
        // Odd length: middle value
        values[len / 2]
    }
}

fn normalize_psd<T: Float>(
    mut psd: Vec<T>,
    window: &[T],
    sample_rate: T,
    nfft: usize,
) -> Vec<T> {
    // Scaling factor: nfft^2 / (fs * sum(window^2))
    // The nfft^2 accounts for FFT magnitude scaling
    let window_power = window_energy(window);
    let nfft_float = T::from(nfft).unwrap();
    let scale = (nfft_float * nfft_float) / (sample_rate * window_power);

    // Apply scaling
    for val in psd.iter_mut() {
        *val = *val * scale;
    }

    psd
}

    fn test_welch_basic() {
        // Test with a simple CW tone
        let sample_rate = 1e6;
        let freq = 1e5; // 100 kHz tone
        let block_size = 4096;

        let mut cw = CW::new(freq, sample_rate, block_size);
        let signal = cw.generate_block();

        let nperseg = 512;
        let (freqs, psd) = welch(
            &signal.iter().cloned().collect::<Vec<_>>(),
            sample_rate,
            nperseg,
            None,
            None,
            WindowType::Hann,
            None,
        );

        // Should have positive frequencies
        assert!(freqs.len() > 0);
        assert!(psd.len() == freqs.len());

        // Find peak
        let (peak_idx, peak_power) = psd
            .iter()
            .enumerate()
            .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
            .unwrap();

        let peak_freq = freqs[peak_idx];

        // Peak should be near the CW frequency
        assert!(
            (peak_freq - freq).abs() < 5e3,

    fn test_welch_cw_tone() {
        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping Welch CW tone plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Welch PSD: CW Tone Test ===");

        let sample_rate = 1e6_f64;
        let freq = 2e5_f64; // 200 kHz
        let block_size = 8192;

        let mut cw = CW::new(freq, sample_rate, block_size);
        let signal = cw.generate_block();

        let nperseg = 1024;
        let (freqs, psd) = welch(
            &signal.iter().cloned().collect::<Vec<_>>(),
            sample_rate,
            nperseg,
            None,
            None,
            WindowType::Hann,
            None,
        );

        println!("Sample rate: {} Hz", sample_rate);
        println!("CW frequency: {} Hz", freq);
        println!("Segment length: {}", nperseg);
        println!("Number of frequency bins: {}", freqs.len());

    fn test_welch_two_tones() {
        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping Welch two tones plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Welch PSD: Two Tones Test ===");

        let sample_rate = 1e6_f64;
        let freq1 = 1.5e5_f64; // 150 kHz
        let freq2 = 3e5_f64;   // 300 kHz
        let block_size = 8192;

        // Generate two CW tones
        let mut cw1 = CW::new(freq1, sample_rate, block_size);
        let mut cw2 = CW::new(freq2, sample_rate, block_size);
        let signal1 = cw1.generate_block::<f64>();
        let signal2 = cw2.generate_block::<f64>();

        // Add them together
        let mut combined: Vec<Complex<f64>> = Vec::with_capacity(block_size);
        for i in 0..block_size {
            combined.push(signal1[i] + signal2[i]);
        }

        let nperseg = 1024;
        let (freqs, psd) = welch(

    fn test_welch_awgn() {
        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping Welch AWGN plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Welch PSD: AWGN Test ===");

        let sample_rate = 1e6_f64;
        let block_size = 16384;
        let noise_power = 1.0;

        let mut awgn = AWGN::new_from_seed(sample_rate, block_size, noise_power, 42);
        let signal = awgn.generate_block();

        let nperseg = 2048;
        let (freqs, psd) = welch(
            &signal.iter().cloned().collect::<Vec<_>>(),
            sample_rate,
            nperseg,
            None,
            None,
            WindowType::Hann,
            None,
        );

        println!("Sample rate: {} Hz", sample_rate);
        println!("Noise power: {}", noise_power);
        println!("Segment length: {}", nperseg);

        // Convert to dB
        use crate::vecto

    fn test_welch_fsk_signal() {
        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping Welch FSK signal plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Welch PSD: FSK Signal Test ===");

        let sample_rate = 1e6_f64;
        let symbol_rate = 1e5_f64;
        let carrier_freq = 2.5e5_f64;
        let deviation = 5e4_f64;
        let block_size = 16384;

        // Generate FSK signal
        let mut fsk: FskCarrier<f64> = FskCarrier::new(
            sample_rate,
            symbol_rate,
            carrier_freq,
            deviation,
            block_size,
            Some(42),
        );
        let signal = fsk.generate_block();

        // Add some noise
        let snr_db = 20.0;
        let signal_power = 1.0; // FSK has unit magnitude
        let snr_linear = 10.0_f64.powf(snr_db / 10.0);
        let noise_power = signal_power / snr_lin

    fn test_welch_psk_signal() {
        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping Welch PSK signal plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Welch PSD: QPSK Signal Test ===");

        let sample_rate = 10e6_f64;
        let symbol_rate = 1e6_f64;
        let rolloff = 0.35;
        let block_size = 16384;
        let filter_taps = 51;

        let mut psk = PskCarrier::new(
            sample_rate,
            symbol_rate,
            ModType::_QPSK,
            rolloff,
            block_size,
            filter_taps,
            Some(42),
        );
        let signal = psk.generate_block();

        let nperseg = 2048;
        let (freqs, psd) = welch(
            &signal.iter().cloned().collect::<Vec<_>>(),
            sample_rate,
            nperseg,
            None,
            None,
            WindowType::Hann,
            None,

    fn test_welch_window_comparison() {
        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping Welch window comparison plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Welch PSD: Window Comparison Test ===");

        let sample_rate = 1e6_f64;
        let freq = 2e5_f64;
        let block_size = 8192;

        let mut cw = CW::new(freq, sample_rate, block_size);
        let signal = cw.generate_block();
        let signal_vec: Vec<Complex<f64>> = signal.iter().cloned().collect();

        let nperseg = 1024;

        // Compute PSD with different windows
        let (freqs_hann, psd_hann) = welch(
            &signal_vec,
            sample_rate,
            nperseg,
            None,
            None,
            WindowType::Hann,
            None,
        );

        let (_, psd_hamming) = welch(
            &signal_vec,
            sample_rate,
   

    fn test_welch_averaging_comparison() {
        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping Welch averaging comparison plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Welch PSD: Averaging Method Comparison Test ===");

        let sample_rate = 1e6_f64;
        let freq = 2e5_f64;
        let block_size = 16384;

        // Generate CW with noise
        let mut cw = CW::new(freq, sample_rate, block_size);
        let signal = cw.generate_block::<f64>();

        let snr_db = 15.0;
        let signal_power = 1.0;
        let snr_linear = 10.0_f64.powf(snr_db / 10.0);
        let noise_power = signal_power / snr_linear;

        let mut awgn = AWGN::new_from_seed(sample_rate, block_size, noise_power, 77);
        let noise = awgn.generate_block();

        let mut noisy_signal: Vec<Complex<f64>> = Vec::with_capacity(block_size);
        for i in 

    fn test_welch_parameter_effects() {
        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping Welch parameter effects plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Welch PSD: Parameter Effects Test ===");

        let sample_rate = 1e6_f64;
        let freq = 2e5_f64;
        let block_size = 16384;

        let mut cw = CW::new(freq, sample_rate, block_size);
        let signal = cw.generate_block();
        let signal_vec: Vec<Complex<f64>> = signal.iter().cloned().collect();

        // Test different segment lengths
        let (freqs_256, psd_256) = welch(
            &signal_vec,
            sample_rate,
            256,
            None,
            None,
            WindowType::Hann,
            None,
        );

        let (freqs_512, psd_512) = welch(
            &signal_vec,
            sample_rate,
            512,
            None,
  

    fn test_welch_overlap_effects() {
        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping Welch overlap effects plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Welch PSD: Overlap Effects Test ===");

        let sample_rate = 1e6_f64;
        let freq = 2e5_f64;
        let block_size = 16384;

        let mut cw = CW::new(freq, sample_rate, block_size);
        let signal = cw.generate_block();
        let signal_vec: Vec<Complex<f64>> = signal.iter().cloned().collect();

        let nperseg = 1024;

        // Test different overlap amounts
        let (freqs_0, psd_0) = welch(
            &signal_vec,
            sample_rate,
            nperseg,
            Some(0), // 0% overlap
            None,
            WindowType::Hann,
            None,
        );

        let (freqs_50, psd_50) = welch(
            &signal_vec,
            sample_rate,

    fn test_parseval_theorem_conservation() {
        // Test that total power is conserved between time and frequency domains
        // This validates the PSD normalization is correct per Parseval's theorem:
        // sum(|x[n]|^2) / N = integral(PSD(f) df) from 0 to fs

        let sample_rate = 1e6_f64;
        let freq = 2e5_f64; // 200 kHz tone
        let block_size = 8192;

        // Generate a simple CW tone
        let mut cw = CW::new(freq, sample_rate, block_size);
        let signal = cw.generate_block::<f64>();

        // Calculate time-domain power (average power)
        let time_power: f64 = signal
            .iter()
            .map(|s| s.norm_sqr())
            .sum::<f64>() / block_size as f64;

        // Calculate frequency-domain power by integrating PSD
        let nperseg = 1024;
        let (freqs, psd) = welch(
            &signal.iter().cloned().collect::<Vec<_>>(),
            sample_rate,
            nperseg,
            None,
            None,
       

    fn test_parseval_theorem_with_noise() {
        // Test Parseval's theorem with a signal containing both tone and AWGN
        // This is more realistic and tests the PSD normalization under typical conditions

        let sample_rate = 1e6_f64;
        let freq = 2e5_f64; // 200 kHz tone
        let block_size = 16384;

        // Generate CW tone
        let mut cw = CW::new(freq, sample_rate, block_size);
        let signal = cw.generate_block::<f64>();

        // Add AWGN with known power
        let snr_db = 10.0;
        let signal_power = 1.0;
        let snr_linear = 10.0_f64.powf(snr_db / 10.0);
        let noise_power = signal_power / snr_linear;

        let mut awgn = AWGN::new_from_seed(sample_rate, block_size, noise_power, 42);
        let noise = awgn.generate_block::<f64>();

        let mut noisy_signal: Vec<Complex<f64>> = Vec::with_capacity(block_size);
        for i in 0..block_size {
            noisy_signal.push(signal[i] + noise[i]);
        }

        // Ca

pub fn generate_window<T: Float>(window_type: WindowType, size: usize) -> Vec<T> {
    match window_type {
        WindowType::Rectangular => rectangular_window(size),
        WindowType::Hann => hann_window(size),
        WindowType::Hamming => hamming_window(size),
        WindowType::Blackman => blackman_window(size),
    }
}

fn hann_window<T: Float>(size: usize) -> Vec<T> {
    if size == 1 {
        return vec![T::one()];
    }

    let pi = T::from(PI).unwrap();
    let two = T::from(2.0).unwrap();
    let half = T::from(0.5).unwrap();
    let n_minus_1 = T::from(size - 1).unwrap();

    (0..size)
        .map(|i| {
            let t = T::from(i).unwrap();
            let arg = (two * pi * t) / n_minus_1;
            half * (T::one() - arg.cos())
        })
        .collect()
}

fn hamming_window<T: Float>(size: usize) -> Vec<T> {
    if size == 1 {
        return vec![T::one()];
    }

    let pi = T::from(PI).unwrap();
    let two = T::from(2.0).unwrap();
    let alpha = T::from(0.54).unwrap();
    let beta = T::from(0.46).unwrap();
    let n_minus_1 = T::from(size - 1).unwrap();

    (0..size)
        .map(|i| {
            let t = T::from(i).unwrap();
            let arg = (two * pi * t) / n_minus_1;
            alpha - beta * arg.cos()
        })
        .collect()
}

fn blackman_window<T: Float>(size: usize) -> Vec<T> {
    if size == 1 {
        return vec![T::one()];
    }

    let pi = T::from(PI).unwrap();
    let two = T::from(2.0).unwrap();
    let four = T::from(4.0).unwrap();
    let a0 = T::from(0.42).unwrap();
    let a1 = T::from(0.5).unwrap();
    let a2 = T::from(0.08).unwrap();
    let n_minus_1 = T::from(size - 1).unwrap();

    (0..size)
        .map(|i| {
            let t = T::from(i).unwrap();
            let arg1 = (two * pi * t) / n_minus_1;
            let arg2 = (four * pi * t) / n_minus_1;
            a0 - a1 * arg1.cos() + a2 * arg2.cos()
        })
        .collect()
}

    fn test_window_properties() {
        // Test basic properties of windows
        let size = 256;

        // Rectangular window
        let rect: Vec<f64> = generate_window(WindowType::Rectangular, size);
        assert_eq!(rect.len(), size);
        assert!((rect[0] - 1.0).abs() < 1e-10);
        assert!((rect[size - 1] - 1.0).abs() < 1e-10);

        // Hann window
        let hann: Vec<f64> = generate_window(WindowType::Hann, size);
        assert_eq!(hann.len(), size);
        assert!(hann[0].abs() < 1e-10, "Hann window should start near 0");
        assert!(hann[size - 1].abs() < 1e-10, "Hann window should end near 0");
        assert!(hann[size / 2] > 0.99, "Hann window peak should be near 1");

        // Hamming window
        let hamming: Vec<f64> = generate_window(WindowType::Hamming, size);
        assert_eq!(hamming.len(), size);
        assert!((hamming[0] - 0.08).abs() < 0.01, "Hamming window should start near 0.08");
        assert!(hamming[size / 2] > 0.99, "Hammin

    fn test_window_energy_sum() {
        let size = 128;

        // Rectangular window
        let rect: Vec<f64> = generate_window(WindowType::Rectangular, size);
        let rect_energy = window_energy(&rect);
        let rect_sum = window_sum(&rect);
        assert!((rect_energy - size as f64).abs() < 1e-10);
        assert!((rect_sum - size as f64).abs() < 1e-10);

        // Hann window
        let hann: Vec<f64> = generate_window(WindowType::Hann, size);
        let hann_energy = window_energy(&hann);
        let hann_sum = window_sum(&hann);
        // Hann window sum should be approximately N/2
        assert!((hann_sum - (size as f64 / 2.0)).abs() < 1.0);
        // Hann window energy should be approximately 3N/8
        assert!((hann_energy - (3.0 * size as f64 / 8.0)).abs() < 1.0);
    }

    fn test_window_shapes() {
        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping window shapes plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Window Shapes Test ===");

        let size = 256;
        let indices: Vec<f64> = (0..size).map(|i| i as f64).collect();

        // Generate all window types
        let rect: Vec<f64> = generate_window(WindowType::Rectangular, size);
        let hann: Vec<f64> = generate_window(WindowType::Hann, size);
        let hamming: Vec<f64> = generate_window(WindowType::Hamming, size);
        let blackman: Vec<f64> = generate_window(WindowType::Blackman, size);

        println!("Generated {} point windows", size);

        // Plot using plotly
        use plotly::{Plot, Scatter};
        use plotly::common::Mode;
        use plotly::layout::{Axis, Layout};

        let trace_rect = Scatter::new(indices.clone(),

    fn test_window_frequency_response() {
        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping window frequency response plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Window Frequency Response Test ===");

        use num_complex::Complex;
        use crate::fft::{fft, fftshift, fftfreqs};
        use crate::vector_ops;

        let size = 1024;

        // Generate windows
        let rect: Vec<f64> = generate_window(WindowType::Rectangular, size);
        let hann: Vec<f64> = generate_window(WindowType::Hann, size);
        let hamming: Vec<f64> = generate_window(WindowType::Hamming, size);
        let blackman: Vec<f64> = generate_window(WindowType::Blackman, size);

        // Compute FFT of each window
        let compute_window_fft = |window: &[f64]| -> Vec<f64> {
            let mut window_complex: Vec<Complex<f64>> = window
               

pub fn bpsk_map<T: Float>() -> HashMap<u8, Complex<T>> {
    let mut map = HashMap::new();
    let zero = T::from(0).unwrap();
    let one = T::from(1).unwrap();

    map.insert(0b0, Complex::new(one, zero));
    map.insert(0b1, Complex::new(-one, zero));

    map
}

pub fn qpsk_gray_map<T: Float>() -> HashMap<u8, Complex<T>> {
    let mut map = HashMap::new();
    let scale = T::from(1.0).unwrap() / T::from(2.0).unwrap().sqrt();

    // Gray coding: 00, 01, 11, 10
    map.insert(0b00, Complex::new(scale, scale));      // Q1
    map.insert(0b01, Complex::new(-scale, scale));     // Q2
    map.insert(0b11, Complex::new(-scale, -scale));    // Q3
    map.insert(0b10, Complex::new(scale, -scale));     // Q4

    map
}

pub fn psk8_gray_map<T: Float>() -> HashMap<u8, Complex<T>> {
    let mut map = HashMap::new();

    // Gray coding for 8-PSK, starting at π/8
    let angles = [
        (0b000, PI / 8.0),
        (0b001, 3.0 * PI / 8.0),
        (0b011, 5.0 * PI / 8.0),
        (0b010, 7.0 * PI / 8.0),
        (0b110, 9.0 * PI / 8.0),
        (0b111, 11.0 * PI / 8.0),
        (0b101, 13.0 * PI / 8.0),
        (0b100, 15.0 * PI / 8.0),
    ];

    for (bits, angle) in angles.iter() {
        let i = T::from(angle.cos()).unwrap();
        let q = T::from(angle.sin()).unwrap();
        map.insert(*bits, Complex::new(i, q));
    }

    map
}

pub fn apsk16_gray_map<T: Float>() -> HashMap<u8, Complex<T>> {
    let mut map = HashMap::new();

    // DVB-S2 radii ratio γ = 2.6 for 16-APSK
    let r1 = T::from(0.5).unwrap();
    let r2 = T::from(1.0).unwrap();

    // Inner ring: 4 symbols at 0°, 90°, 180°, 270°
    let inner_angles = [
        (0b0000, 0.0),
        (0b0001, PI / 2.0),
        (0b0011, PI),
        (0b0010, 3.0 * PI / 2.0),
    ];

    for (bits, angle) in inner_angles.iter() {
        let i = r1 * T::from(angle.cos()).unwrap();
        let q = r1 * T::from(angle.sin()).unwrap();
        map.insert(*bits, Complex::new(i, q));
    }

    // Outer ring: 12 symbols
    let outer_angles = [
        (0b0110, PI / 12.0),
        (0b0111, 3.0 * PI / 12.0),
        (0b0101, 5.0 * PI / 12.0),
        (0b0100, 7.0 * PI / 12.0),
        (0b1100, 9.0 * PI / 12.0),
        (0b1101, 11.0 * PI / 12.0),
        (0b1111, 13.0 * PI / 12.0),
        (0b1110, 15.0 * PI / 12.0),
        (0b1010, 17.0 * PI / 12.0),
        (0b1011, 

pub fn qam16_gray_map<T: Float>() -> HashMap<u8, Complex<T>> {
    let mut map = HashMap::new();

    // Normalized to unit average power
    let scale = T::from(1.0 / 10.0_f64.sqrt()).unwrap();
    let levels = [-3.0, -1.0, 1.0, 3.0];

    for (idx, &i_val) in levels.iter().enumerate() {
        for (jdx, &q_val) in levels.iter().enumerate() {
            // Gray code mapping for 16-QAM
            let i_bits = match idx {
                0 => 0b10,  // -3
                1 => 0b11,  // -1
                2 => 0b01,  //  1
                3 => 0b00,  //  3
                _ => 0,
            };
            let q_bits = match jdx {
                0 => 0b10,  // -3
                1 => 0b11,  // -1
                2 => 0b01,  //  1
                3 => 0b00,  //  3
                _ => 0,
            };

            let bits = (i_bits << 2) | q_bits;
            let i = scale * T::from(i_val).unwrap();
            let q = scale * T::from(q_val).unwrap();
            map.insert(bits, Co

pub fn qam32_gray_map<T: Float>() -> HashMap<u8, Complex<T>> {
    let mut map = HashMap::new();

    // 32-QAM uses a cross constellation pattern
    // Normalized approximately for unit average power
    let scale = T::from(0.3).unwrap();

    // Define the cross-shaped constellation points with Gray coding
    let points = [
        // Center cross (horizontal)
        (0b00000, 0.0, 0.0),
        (0b00001, 2.0, 0.0),
        (0b00011, 4.0, 0.0),
        (0b00010, -2.0, 0.0),
        (0b00110, -4.0, 0.0),

        // Center cross (vertical)
        (0b00111, 0.0, 2.0),
        (0b00101, 0.0, 4.0),
        (0b00100, 0.0, -2.0),
        (0b01100, 0.0, -4.0),

        // Upper right quadrant
        (0b01101, 2.0, 2.0),
        (0b01111, 4.0, 2.0),
        (0b01110, 2.0, 4.0),

        // Upper left quadrant
        (0b01010, -2.0, 2.0),
        (0b01011, -4.0, 2.0),
        (0b01001, -2.0, 4.0),

        // Lower right quadrant
        (0b01000, 2.0, -2.0),
        (0b11000, 4.0, -2.0

pub fn qam64_gray_map<T: Float>() -> HashMap<u8, Complex<T>> {
    let mut map = HashMap::new();

    // Normalized to unit average power
    let scale = T::from(1.0 / 42.0_f64.sqrt()).unwrap();
    let levels = [-7.0, -5.0, -3.0, -1.0, 1.0, 3.0, 5.0, 7.0];

    for (idx, &i_val) in levels.iter().enumerate() {
        for (jdx, &q_val) in levels.iter().enumerate() {
            // Gray code mapping for 64-QAM
            let i_bits = match idx {
                0 => 0b100,  // -7
                1 => 0b101,  // -5
                2 => 0b111,  // -3
                3 => 0b110,  // -1
                4 => 0b010,  //  1
                5 => 0b011,  //  3
                6 => 0b001,  //  5
                7 => 0b000,  //  7
                _ => 0,
            };
            let q_bits = match jdx {
                0 => 0b100,  // -7
                1 => 0b101,  // -5
                2 => 0b111,  // -3
                3 => 0b110,  // -1
                4 => 0b010,  //  1
                5 =

    fn test_bpsk_map() {
        let map = bpsk_map::<f32>();
        assert_eq!(map.len(), 2);

        for val in map.values() {
            let power = val.norm_sqr();
            assert!((power - 1.0).abs() < 1e-10);
        }
    }

    fn test_qpsk_gray_map() {
        let map = qpsk_gray_map::<f64>();
        assert_eq!(map.len(), 4);

        // Verify unit power normalization
        for val in map.values() {
            let power = val.norm_sqr();
            assert!((power - 1.0).abs() < 1e-10);
        }
    }

    fn test_psk8_gray_map() {
        let map = psk8_gray_map::<f64>();
        assert_eq!(map.len(), 8);

        // Verify unit circle
        for val in map.values() {
            let radius = val.norm();
            assert!((radius - 1.0).abs() < 1e-10);
        }
    }

    fn test_qpsk_constellation_plot() {
        use std::env;
        use crate::generate::random_bit_generator::BitGenerator;
        use crate::plot::plot_constellation;

        let plot = env::var("PLOT").unwrap_or_else(|_| "false".to_string());
        if plot.to_lowercase() != "true" {
            println!("Skipping constellation plot (set PLOT=true to enable)");
            return;
        }

        // Get QPSK constellation map
        let map = qpsk_gray_map::<f64>();

        // Generate actual QPSK symbols
        let mut bit_gen = BitGenerator::new_from_seed(42);
        let num_symbols = 1000;
        let mut symbols = Vec::with_capacity(num_symbols);

        for _ in 0..num_symbols {
            let bits = bit_gen.next_2_bits();
            if let Some(symbol) = map.get(&bits) {
                symbols.push(*symbol);
            }
        }

        // Plot the constellation using the utility function
        plot_constellation(&symbols, "QPSK Constellation");
    }

pub fn should_plot() -> bool {
    env::var("PLOT")
        .unwrap_or_else(|_| "false".to_string())
        .to_lowercase()
        == "true"
}

pub fn max<T: Clone + Copy + PartialOrd>(vector: &[T]) -> (usize, T) {
    let mut max_idx = 0;
    let mut max_val = vector[0].clone();

    for (i, &val) in vector.iter().enumerate().skip(1) {
        if val > max_val {
            max_val = val.clone();
            max_idx = i;
        }
    }

    (max_idx, max_val)
}

pub fn to_db<T: Float + 'static>(vector: &[T]) -> Vec<T> {
    let mut v = Vec::new();
    v.reserve(vector.len());

    let ten = T::from(10.0).unwrap();
    if TypeId::of::<T>() == TypeId::of::<f32>() {
        v = vector.iter().map(|val| ten * val.log10()).collect();
    } else if TypeId::of::<T>() == TypeId::of::<f64>() {
        v = vector.iter().map(|val| ten * val.log10()).collect();
    }

    v
}

pub fn add<T>(a: &[T], b: &[T]) -> Vec<T>
where
    T: std::ops::Add<Output = T> + Copy,
{
    assert_eq!(a.len(), b.len());
    a.iter().zip(b.iter()).map(|(x, y)| *x + *y).collect()
}

    fn test_add_vectors() {
        let a = vec![1, 2, 3];
        let b = vec![4, 5, 6];
        let c = add(&a, &b);
        assert_eq!(c, vec![5, 7, 9]);
    }

    fn process_simd_batches<B>(
        &self,
        other: &VectorSimd<T>,
        batch_op: impl Fn(B, B) -> B,
        scalar_op: impl Fn(T, T) -> T,
    ) -> VectorSimd<T>
    where
        B: SimdBatch<Scalar = T> + Pod,
    {
        assert_eq!(
            self.len(),
            other.len(),
            "Vector sizes must match for element-wise operations"
        );

        // Allocate without initialization for better performance
        let mut result = Vec::with_capacity(self.len());
        unsafe {
            result.set_len(self.len());
        }

        let lanes = B::LANES;
        let simd_end = (self.len() / lanes) * lanes;

        // Process SIMD batches using zero-copy cast_slice
        let a_simd: &[B] = cast_slice(&self.data[..simd_end]);
        let b_simd: &[B] = cast_slice(&other.data[..simd_end]);
        let result_simd: &mut [B] = cast_slice_mut(&mut result[..simd_end]);

        // Use iterator-based approach for better optimization
        result_si

    fn process_simd_batches_scalar<B>(
        &self,
        scalar: T,
        batch_op: impl Fn(B, B) -> B,
        scalar_op: impl Fn(T, T) -> T,
    ) -> VectorSimd<T>
    where
        B: SimdBatch<Scalar = T> + Pod,
    {
        // Allocate without initialization for better performance
        let mut result = Vec::with_capacity(self.len());
        unsafe {
            result.set_len(self.len());
        }

        let lanes = B::LANES;
        let simd_end = (self.len() / lanes) * lanes;
        let scalar_batch = B::splat(scalar);

        // Process SIMD batches using zero-copy cast_slice
        let a_simd: &[B] = cast_slice(&self.data[..simd_end]);
        let result_simd: &mut [B] = cast_slice_mut(&mut result[..simd_end]);

        // Use iterator-based approach for better optimization
        result_simd.iter_mut()
            .zip(a_simd.iter())
            .for_each(|(r, a)| *r = batch_op(*a, scalar_batch));

        // Process remainder
        for i in simd_end..sel

    fn process_simd_batches_scalar_inplace<B>(
        &mut self,
        scalar: T,
        batch_op: impl Fn(B, B) -> B,
        scalar_op: impl Fn(T, T) -> T,
    )
    where
        B: SimdBatch<Scalar = T> + Pod,
    {
        let lanes = B::LANES;
        let simd_end = (self.len() / lanes) * lanes;
        let scalar_batch = B::splat(scalar);

        // Process SIMD batches using zero-copy cast_slice
        let data_simd: &mut [B] = cast_slice_mut(&mut self.data[..simd_end]);

        // Use iterator-based approach for better optimization
        data_simd.iter_mut()
            .for_each(|d| *d = batch_op(*d, scalar_batch));

        // Process remainder
        for i in simd_end..self.len() {
            self.data[i] = scalar_op(self.data[i], scalar);
        }
    }

    pub fn linspace(start: f32, stop: f32, length: usize) -> Self {
        if length == 0 {
            return Self::new();
        }
        if length == 1 {
            return Self::from_vec(vec![start]);
        }

        let step = (stop - start) / (length as f32 - 1.0);
        let data: Vec<f32> = (0..length)
            .map(|i| start + (i as f32) * step)
            .collect();
        Self { data }
    }

    fn add(self, rhs: Self) -> Self::Output {
        use simd_config::F32Batch;
        self.process_simd_batches::<F32Batch>(
            &rhs,
            |a, b| SimdBatch::add(a, b),
            |a, b| a + b,
        )
    }

    fn add(self, rhs: Self) -> Self::Output {
        use simd_config::F32Batch;
        self.process_simd_batches::<F32Batch>(
            rhs,
            |a, b| SimdBatch::add(a, b),
            |a, b| a + b,
        )
    }

    fn sub(self, rhs: Self) -> Self::Output {
        use simd_config::F32Batch;
        self.process_simd_batches::<F32Batch>(
            &rhs,
            |a, b| SimdBatch::sub(a, b),
            |a, b| a - b,
        )
    }

    fn sub(self, rhs: Self) -> Self::Output {
        use simd_config::F32Batch;
        self.process_simd_batches::<F32Batch>(
            rhs,
            |a, b| SimdBatch::sub(a, b),
            |a, b| a - b,
        )
    }

    fn mul(self, rhs: Self) -> Self::Output {
        use simd_config::F32Batch;
        self.process_simd_batches::<F32Batch>(
            &rhs,
            |a, b| SimdBatch::mul(a, b),
            |a, b| a * b,
        )
    }

    fn mul(self, rhs: Self) -> Self::Output {
        use simd_config::F32Batch;
        self.process_simd_batches::<F32Batch>(
            rhs,
            |a, b| SimdBatch::mul(a, b),
            |a, b| a * b,
        )
    }

    fn div(self, rhs: Self) -> Self::Output {
        use simd_config::F32Batch;
        self.process_simd_batches::<F32Batch>(
            &rhs,
            |a, b| SimdBatch::div(a, b),
            |a, b| a / b,
        )
    }