Function bodies 301 total

    fn test_rrc_filter_frequency_response() {
        use std::env;
        use crate::fft::{fft, fftshift, fftfreqs};
        use crate::vector_ops;

        let sample_rate_hz = 1e6_f64;
        let symbol_rate_hz = 800e3_f64;
        let filter_taps = 51;
        let rolloff_factor = 0.05_f64;

        // Build the RRC filter
        let rrc = RRCFilter::new(
            filter_taps,
            sample_rate_hz,
            symbol_rate_hz,
            rolloff_factor,
        );
        let filter = rrc.build_filter::<f64>();

        println!("\n=== RRC Filter Analysis ===");
        println!("Filter taps: {}", filter.len());
        println!("Sample rate: {} Hz", sample_rate_hz);
        println!("Symbol rate: {} Hz", symbol_rate_hz);
        println!("Samples per symbol (SPS): {}", sample_rate_hz / symbol_rate_hz);
        println!("Rolloff: {}", rolloff_factor);

        // Print all taps to see the shape
        println!("\nAll filter taps:");
        for i in 0..filter.len() {
 

    pub fn new_from_seed(seed: u64) -> Self {
        Self {
            rng: StdRng::seed_from_u64(seed),
            buffer: 0,
            bits_remaining: 0,
        }
    }

    pub fn new_from_entropy() -> Self {
        Self {
            rng: StdRng::from_entropy(),
            buffer: 0,
            bits_remaining: 0,
        }
    }

    pub fn next_bits(&mut self, num_bits: u8) -> u64 {
        assert!(num_bits > 0 && num_bits <= 64, "num_bits must be 1-64");

        // Refill buffer if we don't have enough bits
        if self.bits_remaining < num_bits {
            self.buffer = self.rng.next_u64();
            self.bits_remaining = 64;
        }

        // Extract the requested bits from the low end
        let mask = if num_bits == 64 {
            u64::MAX
        } else {
            (1u64 << num_bits) - 1
        };
        let result = self.buffer & mask;

        // Shift buffer and update remaining count
        self.buffer >>= num_bits;
        self.bits_remaining -= num_bits;

        return result;
    }

    pub fn get_mod_type_from_enum<T: Float>(mod_enum: ModType) -> Modulation<T> {
        match mod_enum {
            ModType::_BPSK => Modulation::BPSK(_BPSK::new(bpsk_map())),
            ModType::_QPSK => Modulation::QPSK(_QPSK::new(qpsk_gray_map())),
            ModType::_8PSK => Modulation::PSK8(_8PSK::new(psk8_gray_map())),
            ModType::_16APSK => Modulation::APSK16(_16APSK::new(apsk16_gray_map())),
            ModType::_16QAM => Modulation::QAM16(_16QAM::new(qam16_gray_map())),
            ModType::_32QAM => Modulation::QAM32(_32QAM::new(qam32_gray_map())),
            ModType::_64QAM => Modulation::QAM64(_64QAM::new(qam64_gray_map())),
            ModType::_CW => unreachable!("CW is handled separately by Carrier, not through modulation system"),
        }
    }

        pub fn new(bit_symbol_map: HashMap<u8, Complex<T>>) -> Self {
            _BPSK {
                mod_properties: ModProperties {
                    bits_per_symbol: 1,
                    num_symbols: 2,
                    bit_symbol_map
                }
            }
        }

        fn demodulate(&self, symbol: &Complex<T>) -> Option<u8> {
            self.mod_properties.bit_symbol_map
                .iter()
                .min_by(|(_, s1), (_, s2)| {
                    let d1 = (*s1 - *symbol).norm_sqr();
                    let d2 = (*s2 - *symbol).norm_sqr();
                    d1.partial_cmp(&d2).unwrap()
                })
                .map(|(bits, _)| *bits)
        }

        pub fn new(bit_symbol_map: HashMap<u8, Complex<T>>) -> Self {
            _QPSK {
                mod_properties: ModProperties {
                    bits_per_symbol: 2,
                    num_symbols: 4,
                    bit_symbol_map,
                }
            }
        }

        fn demodulate(&self, symbol: &Complex<T>) -> Option<u8> {
            self.mod_properties.bit_symbol_map
                .iter()
                .min_by(|(_, s1), (_, s2)| {
                    let d1 = (*s1 - *symbol).norm_sqr();
                    let d2 = (*s2 - *symbol).norm_sqr();
                    d1.partial_cmp(&d2).unwrap()
                })
                .map(|(bits, _)| *bits)
        }

        pub fn new(bit_symbol_map: HashMap<u8, Complex<T>>) -> Self {
            _8PSK {
                mod_properties: ModProperties {
                    bits_per_symbol: 3,
                    num_symbols: 8,
                    bit_symbol_map,
                }
            }
        }

        fn demodulate(&self, symbol: &Complex<T>) -> Option<u8> {
            self.mod_properties.bit_symbol_map
                .iter()
                .min_by(|(_, s1), (_, s2)| {
                    let d1 = (*s1 - *symbol).norm_sqr();
                    let d2 = (*s2 - *symbol).norm_sqr();
                    d1.partial_cmp(&d2).unwrap()
                })
                .map(|(bits, _)| *bits)
        }

        pub fn new(bit_symbol_map: HashMap<u8, Complex<T>>) -> Self {
            _16APSK {
                mod_properties: ModProperties {
                    bits_per_symbol: 4,
                    num_symbols: 16,
                    bit_symbol_map,
                }
            }
        }

        fn demodulate(&self, symbol: &Complex<T>) -> Option<u8> {
            self.mod_properties.bit_symbol_map
                .iter()
                .min_by(|(_, s1), (_, s2)| {
                    let d1 = (*s1 - *symbol).norm_sqr();
                    let d2 = (*s2 - *symbol).norm_sqr();
                    d1.partial_cmp(&d2).unwrap()
                })
                .map(|(bits, _)| *bits)
        }

        pub fn new(bit_symbol_map: HashMap<u8, Complex<T>>) -> Self {
            _16QAM {
                mod_properties: ModProperties {
                    bits_per_symbol: 4,
                    num_symbols: 16,
                    bit_symbol_map,
                }
            }
        }

        fn demodulate(&self, symbol: &Complex<T>) -> Option<u8> {
            self.mod_properties.bit_symbol_map
                .iter()
                .min_by(|(_, s1), (_, s2)| {
                    let d1 = (*s1 - *symbol).norm_sqr();
                    let d2 = (*s2 - *symbol).norm_sqr();
                    d1.partial_cmp(&d2).unwrap()
                })
                .map(|(bits, _)| *bits)
        }

        pub fn new(bit_symbol_map: HashMap<u8, Complex<T>>) -> Self {
            _32QAM {
                mod_properties: ModProperties {
                    bits_per_symbol: 5,
                    num_symbols: 32,
                    bit_symbol_map,
                }
            }
        }

        fn demodulate(&self, symbol: &Complex<T>) -> Option<u8> {
            self.mod_properties.bit_symbol_map
                .iter()
                .min_by(|(_, s1), (_, s2)| {
                    let d1 = (*s1 - *symbol).norm_sqr();
                    let d2 = (*s2 - *symbol).norm_sqr();
                    d1.partial_cmp(&d2).unwrap()
                })
                .map(|(bits, _)| *bits)
        }

        pub fn new(bit_symbol_map: HashMap<u8, Complex<T>>) -> Self {
            _64QAM {
                mod_properties: ModProperties {
                    bits_per_symbol: 6,
                    num_symbols: 64,
                    bit_symbol_map,
                }
            }
        }

        fn demodulate(&self, symbol: &Complex<T>) -> Option<u8> {
            self.mod_properties.bit_symbol_map
                .iter()
                .min_by(|(_, s1), (_, s2)| {
                    let d1 = (*s1 - *symbol).norm_sqr();
                    let d2 = (*s2 - *symbol).norm_sqr();
                    d1.partial_cmp(&d2).unwrap()
                })
                .map(|(bits, _)| *bits)
        }

        fn test_bpsk_modulate() {
            let bpsk = _BPSK::<f64>::new(crate::symbol_maps::bpsk_map());

            let symbol = bpsk.modulate(0b0).unwrap();
            assert!(symbol.re > 0.0);
            assert_eq!(symbol.im, 0.0);

            let symbol = bpsk.modulate(0b1).unwrap();
            assert!(symbol.re < 0.0);
            assert_eq!(symbol.im, 0.0);
        }

        fn test_bpsk_demodulate() {
            let bpsk = _BPSK::<f64>::new(crate::symbol_maps::bpsk_map());

            let symbol = bpsk.modulate(0b0).unwrap();
            let bits = bpsk.demodulate(&symbol).unwrap();
            assert_eq!(bits, 0b0);

            let symbol = bpsk.modulate(0b1).unwrap();
            let bits = bpsk.demodulate(&symbol).unwrap();
            assert_eq!(bits, 0b1);
        }

        fn test_qpsk_modulate() {
            let qpsk = _QPSK::<f64>::new(crate::symbol_maps::qpsk_gray_map());

            let symbol = qpsk.modulate(0b00).unwrap();
            assert!(symbol.re > 0.0 && symbol.im > 0.0);

            let symbol = qpsk.modulate(0b11).unwrap();
            assert!(symbol.re < 0.0 && symbol.im < 0.0);
        }

        fn test_qpsk_demodulate() {
            let qpsk = _QPSK::<f64>::new(crate::symbol_maps::qpsk_gray_map());

            let symbol = qpsk.modulate(0b01).unwrap();
            let bits = qpsk.demodulate(&symbol).unwrap();
            assert_eq!(bits, 0b01);
        }

        fn test_psk8_modulate() {
            let psk8 = _8PSK::<f64>::new(crate::symbol_maps::psk8_gray_map());

            let symbol = psk8.modulate(0b000).unwrap();
            assert!(symbol.norm() > 0.99 && symbol.norm() < 1.01);
        }

        fn test_qam16_modulate_demodulate() {
            let qam16 = _16QAM::<f64>::new(crate::symbol_maps::qam16_gray_map());

            for bits in 0..16u8 {
                let symbol = qam16.modulate(bits).unwrap();
                let decoded = qam16.demodulate(&symbol).unwrap();
                assert_eq!(bits, decoded);
            }
        }

        fn test_get_mod_type_from_enum() {
            let modulation = get_mod_type_from_enum::<f64>(ModType::_QPSK);

            match modulation {
                Modulation::QPSK(qpsk) => {
                    let symbol = qpsk.modulate(0b00).unwrap();
                    assert!(symbol.norm() > 0.0);
                },
                _ => panic!("Expected QPSK variant"),
            }
        }

pub fn plot_spectrum<T: Float>(freqs: &Vec<T>, spectrum: &Vec<T>, title: &str) {
    // Convert to f64 for plotting
    let freqs_f64: Vec<f64> = freqs.iter().map(|f| f.to_f64().unwrap()).collect();
    let spectrum_f64: Vec<f64> = spectrum.iter().map(|s| s.to_f64().unwrap()).collect();

    let mut plot = Plot::new();
    let trace = Scatter::new(freqs_f64, spectrum_f64);
    plot.add_trace(trace);

    let layout = Layout::new()
        .title(title)
        .x_axis(Axis::new().title("Frequency (Hz)"))
        .y_axis(Axis::new().title("Magnitude (dB)"))
        .auto_size(true);
    plot.set_layout(layout);

    plot.show();
}

pub fn plot_constellation<T: Float + std::fmt::Display>(symbols: &[Complex<T>], title: &str) {
    // Extract I and Q values
    let mut i_vals = Vec::new();
    let mut q_vals = Vec::new();

    for symbol in symbols {
        // Convert to f64 for plotting
        i_vals.push(symbol.re.to_f64().unwrap());
        q_vals.push(symbol.im.to_f64().unwrap());
    }

    // Create scatter plot for constellation points
    let trace = Scatter::new(i_vals.clone(), q_vals.clone())
        .mode(Mode::Markers)
        .name("Constellation");

    let mut plot = Plot::new();
    plot.add_trace(trace);

    // Add axis reference lines at 0
    let x_axis_line = Scatter::new(vec![-1.5, 1.5], vec![0.0, 0.0])
        .mode(Mode::Lines)
        .name("I-axis")
        .show_legend(false);
    let y_axis_line = Scatter::new(vec![0.0, 0.0], vec![-1.5, 1.5])
        .mode(Mode::Lines)
        .name("Q-axis")
        .show_legend(false);

    plot.add_trace(x_axis_line);
    plot.add_trace(y_axis_line);

fn parse_modulation(mod_str: &str) -> PyResult<ModType> {
    match mod_str.to_uppercase().as_str() {
        "BPSK" => Ok(ModType::_BPSK),
        "QPSK" => Ok(ModType::_QPSK),
        "8PSK" => Ok(ModType::_8PSK),
        "16APSK" => Ok(ModType::_16APSK),
        "16QAM" => Ok(ModType::_16QAM),
        "32QAM" => Ok(ModType::_32QAM),
        "64QAM" => Ok(ModType::_64QAM),
        "CW" => Ok(ModType::_CW),
        _ => Err(PyErr::new::<pyo3::exceptions::PyValueError, _>(
            format!(
                "Unknown modulation type: '{}'. Supported types (case-insensitive):\n  \
                 BPSK (Binary PSK), QPSK (Quadrature PSK), 8PSK, 16APSK, \n  \
                 16QAM, 32QAM, 64QAM, CW (Continuous Wave)",
                mod_str
            ),
        )),
    }
}

fn parse_impairment(config: &Bound<PyAny>) -> PyResult<Impairment> {
    // Try as string first
    if let Ok(s) = config.extract::<String>() {
        return match s.to_lowercase().as_str() {
            "digitizer_droop_ad9361" => Ok(Impairment::DigitizerDroopAD9361),
            "digitizer_droop_traditional" => Ok(Impairment::DigitizerDroopTraditional),
            "cosine_taper_digitizer" => Ok(Impairment::CosineTaperDigitizer),
            _ => Err(PyErr::new::<pyo3::exceptions::PyValueError, _>(
                format!(
                    "Unknown impairment: '{}'. Supported strings: digitizer_droop_ad9361, digitizer_droop_traditional, cosine_taper_digitizer. Or use a dict with 'type' key.",
                    s
                ),
            )),
        };
    }

    // Try as dictionary
    if let Ok(dict) = config.downcast::<PyDict>() {
        let impairment_type: String = dict
            .get_item("type")?
            .ok_or_else(|| PyErr::new::<pyo3::exceptions::PyValueE

    fn new(
        modulation: &str,
        bandwidth: f64,
        center_freq: f64,
        snr_db: f64,
        rolloff: f64,
        sample_rate_hz: f64,
        seed: Option<u64>,
    ) -> PyResult<Self> {
        let mod_type = parse_modulation(modulation)?;

        // Validate parameters (same validation as Rust API)
        if bandwidth <= 0.0 || bandwidth > 1.0 {
            return Err(PyErr::new::<pyo3::exceptions::PyValueError, _>(
                "bandwidth must be in range (0.0, 1.0]",
            ));
        }
        if center_freq < -0.5 || center_freq > 0.5 {
            return Err(PyErr::new::<pyo3::exceptions::PyValueError, _>(
                "center_freq must be in range [-0.5, 0.5]",
            ));
        }
        if rolloff < 0.0 || rolloff > 1.0 {
            return Err(PyErr::new::<pyo3::exceptions::PyValueError, _>(
                "rolloff must be in range [0.0, 1.0]",
            ));
        }
        if sample_rate_hz <= 0.0 {
            return Err(P

    fn generate(&self, py: Python, num_samples: usize) -> Py<PyArray1<Complex<f64>>> {
        let samples = self.inner.generate::<f64>(num_samples);

        // Convert ComplexVec to Vec<Complex<f64>>
        let vec: Vec<Complex<f64>> = (0..samples.len())
            .map(|i| samples[i])
            .collect();

        // Transfer ownership to Python
        let array = vec.into_pyarray_bound(py);
        array.into()
    }

    fn generate_clean(&self, py: Python, num_samples: usize) -> Py<PyArray1<Complex<f64>>> {
        let samples = self.inner.generate_clean::<f64>(num_samples);

        // Convert ComplexVec to Vec<Complex<f64>>
        let vec: Vec<Complex<f64>> = (0..samples.len())
            .map(|i| samples[i])
            .collect();

        // Transfer ownership to Python
        let array = vec.into_pyarray_bound(py);
        array.into()
    }

    fn new(carriers: Vec<PyRef<PythonCarrier>>) -> PyResult<Self> {
        let rust_carriers: Vec<Carrier> = carriers
            .into_iter()
            .map(|c| c.inner.clone())
            .collect();

        if rust_carriers.is_empty() {
            return Err(PyErr::new::<pyo3::exceptions::PyValueError, _>(
                "Channel must have at least one carrier",
            ));
        }

        Ok(PythonChannel {
            inner: Channel::new(rust_carriers),
        })
    }

    fn set_noise_floor_linear(&mut self, noise_floor_linear: f64) -> PyResult<()> {
        if noise_floor_linear <= 0.0 {
            return Err(PyErr::new::<pyo3::exceptions::PyValueError, _>(
                "noise_floor_linear must be positive",
            ));
        }
        self.inner.set_noise_floor_linear(noise_floor_linear);
        Ok(())
    }

    fn generate(&self, py: Python, num_samples: usize) -> PyResult<Py<PyArray1<Complex<f64>>>> {
        let samples = self.inner.generate::<f64>(num_samples);

        // Convert ComplexVec to Vec<Complex<f64>>
        let vec: Vec<Complex<f64>> = (0..samples.len())
            .map(|i| samples[i])
            .collect();

        // Transfer ownership to Python
        let array = vec.into_pyarray_bound(py);
        Ok(array.into())
    }

    fn generate_clean(&self, py: Python, num_samples: usize) -> Py<PyArray1<Complex<f64>>> {
        let samples = self.inner.generate_clean::<f64>(num_samples);

        // Convert ComplexVec to Vec<Complex<f64>>
        let vec: Vec<Complex<f64>> = (0..samples.len())
            .map(|i| samples[i])
            .collect();

        // Transfer ownership to Python
        let array = vec.into_pyarray_bound(py);
        array.into()
    }

fn parse_window_type(window_str: &str) -> PyResult<WindowType> {
    match window_str.to_lowercase().as_str() {
        "hann" => Ok(WindowType::Hann),
        "hamming" => Ok(WindowType::Hamming),
        "blackman" => Ok(WindowType::Blackman),
        "rectangular" | "rect" => Ok(WindowType::Rectangular),
        _ => Err(PyErr::new::<pyo3::exceptions::PyValueError, _>(
            format!(
                "Unknown window type: '{}'. Supported types: hann, hamming, blackman, rectangular",
                window_str
            ),
        )),
    }
}

fn welch(
    py: Python,
    signal: PyReadonlyArray1<Complex<f64>>,
    sample_rate: f64,
    nperseg: usize,
    noverlap: Option<usize>,
    nfft: Option<usize>,
    window: &str,
) -> PyResult<(Py<PyArray1<f64>>, Py<PyArray1<f64>>)> {
    // Validate parameters
    if nperseg == 0 {
        return Err(PyErr::new::<pyo3::exceptions::PyValueError, _>(
            "nperseg must be greater than 0",
        ));
    }
    if sample_rate <= 0.0 {
        return Err(PyErr::new::<pyo3::exceptions::PyValueError, _>(
            "sample_rate must be positive",
        ));
    }

    let window_type = parse_window_type(window)?;

    // Extract signal from numpy array
    let signal_vec: Vec<Complex<f64>> = signal.as_array().to_vec();

    // Call the Rust implementation
    let (freqs, psd) = welch_impl(
        &signal_vec,
        sample_rate,
        nperseg,
        noverlap,
        nfft,
        window_type,
        None, // Use default averaging method (Mean)
    );

    // Convert re

fn _signal_kit(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
    m.add_class::<PythonCarrier>()?;
    m.add_class::<PythonChannel>()?;
    m.add_function(wrap_pyfunction!(welch, m)?)?;
    m.add("__version__", "0.1.0")?;

    Ok(())
}

pub fn welch<T>(
    signal: &[Complex<T>],
    sample_rate: T,
    nperseg: usize,
    noverlap: Option<usize>,
    nfft: Option<usize>,
    window: WindowType,
    averaging: Option<AveragingMethod>,
) -> (Vec<T>, Vec<T>)
where
    T: Float + std::fmt::Debug + std::ops::RemAssign + std::ops::DivAssign + Send + Sync + num_traits::FromPrimitive + num_traits::Signed + 'static,
{
    // Compute defaults
    let (overlap, fft_len) = compute_defaults(nperseg, noverlap, nfft);
    let avg_method = averaging.unwrap_or(AveragingMethod::Mean);

    // Generate window
    let window_vec = generate_window(window, nperseg);

    // Extract overlapping segments
    let segments = extract_segments(signal, nperseg, overlap);

    if segments.is_empty() {
        // Return empty result if no segments
        return (vec![], vec![]);
    }

    // Compute PSD for each segment
    let psds: Vec<Vec<T>> = segments
        .iter()
        .map(|seg| compute_segment_psd(seg, &window_vec, fft_len))
       

fn compute_defaults(
    nperseg: usize,
    noverlap: Option<usize>,
    nfft: Option<usize>,
) -> (usize, usize) {
    let overlap = noverlap.unwrap_or(nperseg / 2);
    let fft_len = nfft.unwrap_or(nperseg);
    (overlap, fft_len)
}

fn extract_segments<T: Float>(
    signal: &[Complex<T>],
    nperseg: usize,
    noverlap: usize,
) -> Vec<Vec<Complex<T>>> {
    let step = nperseg - noverlap;
    let mut segments = Vec::new();

    let mut start = 0;
    while start + nperseg <= signal.len() {
        let segment = signal[start..start + nperseg].to_vec();
        segments.push(segment);
        start += step;
    }

    segments
}

fn compute_segment_psd<T>(
    segment: &[Complex<T>],
    window: &[T],
    nfft: usize,
) -> Vec<T>
where
    T: Float + std::fmt::Debug + std::ops::RemAssign + std::ops::DivAssign + Send + Sync + num_traits::FromPrimitive + num_traits::Signed + 'static,
{
    // Apply window
    let windowed = apply_window(segment, window);

    // Pad or truncate to nfft length
    let mut fft_input = prepare_fft_input(windowed, nfft);

    // Compute FFT
    fft(&mut fft_input);

    // Compute power: |FFT|^2
    fft_input.iter().map(|c| c.norm_sqr()).collect()
}

fn apply_window<T: Float>(segment: &[Complex<T>], window: &[T]) -> Vec<Complex<T>> {
    segment
        .iter()
        .zip(window.iter())
        .map(|(s, &w)| Complex::new(s.re * w, s.im * w))
        .collect()
}

fn prepare_fft_input<T: Float>(windowed: Vec<Complex<T>>, nfft: usize) -> Vec<Complex<T>> {
    let mut fft_input = windowed;

    if nfft > fft_input.len() {
        // Zero-pad
        fft_input.resize(nfft, Complex::new(T::zero(), T::zero()));
    } else if nfft < fft_input.len() {
        // Truncate
        fft_input.truncate(nfft);
    }

    fft_input
}

fn average_spectra<T: Float>(spectra: Vec<Vec<T>>, method: AveragingMethod) -> Vec<T> {
    match method {
        AveragingMethod::Mean => compute_mean(&spectra),
        AveragingMethod::Median => compute_median(&spectra),
        AveragingMethod::Max => compute_max(&spectra),
        AveragingMethod::Min => compute_min(&spectra),
    }
}

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

    let n_bins = spectra[0].len();
    let n_spectra = T::from(spectra.len()).unwrap();

    (0..n_bins)
        .map(|i| {
            let sum = spectra.iter().fold(T::zero(), |acc, spec| acc + spec[i]);
            sum / n_spectra
        })
        .collect()
}

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

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

    (0..n_bins)
        .map(|i| {
            let mut values: Vec<T> = spectra.iter().map(|spec| spec[i]).collect();
            median_of_slice(&mut values)
        })
        .collect()
}

fn compute_max<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::neg_infinity(), |acc, val| {
                    if val > acc { val } else { acc }
                })
        })
        .collect()
}