Function bodies 301 total

def combine_carriers(*carrier_arrays):
    """
    Combine multiple carrier IQ samples by element-wise addition.

    Parameters:
        *carrier_arrays: Variable number of numpy arrays to combine

    Returns:
        numpy.ndarray: Combined signal

    Example:
        >>> import signal_kit as sk
        >>> carrier1 = sk.Carrier(...).generate(1000)
        >>> carrier2 = sk.Carrier(...).generate(1000)
        >>> combined = sk.combine_carriers(carrier1, carrier2)
    """
    import numpy as np
    return np.sum(carrier_arrays, axis=0)

    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
        match self {
            CafError::InvalidDopplerRange { min, max } => {
                write!(f, "Invalid Doppler range: min ({}) >= max ({})", min, max)
            }
            CafError::InvalidDelayRange { min, max } => {
                write!(f, "Invalid delay range: min ({}) >= max ({})", min, max)
            }
            CafError::SignalLengthMismatch {
                signal1_len,
                signal2_len,
            } => write!(
                f,
                "Signal length mismatch: {} != {}",
                signal1_len, signal2_len
            ),
        }
    }

pub fn compute_caf<T>(
    signal1: &[Complex<T>],
    signal2: &[Complex<T>],
    params: &CafParams,
) -> Result<CafSurface<T>, CafError>
where
    T: Float + RemAssign + DivAssign + Send + Sync + FromPrimitive + Signed + Debug + 'static,
{
    // Validate signal lengths
    if signal1.len() != signal2.len() {
        return Err(CafError::SignalLengthMismatch {
            signal1_len: signal1.len(),
            signal2_len: signal2.len(),
        });
    }

    // Validate Doppler range
    if params.doppler_range_hz.0 >= params.doppler_range_hz.1 {
        return Err(CafError::InvalidDopplerRange {
            min: params.doppler_range_hz.0,
            max: params.doppler_range_hz.1,
        });
    }

    // Validate delay range
    if params.delay_range_samples.0 >= params.delay_range_samples.1 {
        return Err(CafError::InvalidDelayRange {
            min: params.delay_range_samples.0,
            max: params.delay_range_samples.1,
        });
    }

    let n = signal1.len

fn compute_correlation_at_doppler<T>(
    signal1: &[Complex<T>],
    signal2: &[Complex<T>],
    doppler_hz: f64,
    sample_rate_hz: f64,
) -> Vec<Complex<T>>
where
    T: Float + RemAssign + DivAssign + Send + Sync + FromPrimitive + Signed + Debug + 'static,
{
    // Apply frequency shift to signal1
    let mut signal1_shifted: Vec<Complex<T>> = signal1.to_vec();
    freq_shift_inplace(&mut signal1_shifted, doppler_hz, sample_rate_hz);

    // FFT of shifted signal1
    let mut s1_freq = signal1_shifted.clone();
    fft(&mut s1_freq);

    // FFT of signal2
    let mut s2_freq = signal2.to_vec();
    fft(&mut s2_freq);

    // Multiply conjugate: S1* × S2
    let mut correlation_freq: Vec<Complex<T>> = s1_freq
        .iter()
        .zip(s2_freq.iter())
        .map(|(s1, s2)| s1.conj() * s2)
        .collect();

    // IFFT to get correlation in time domain
    ifft(&mut correlation_freq);

    correlation_freq
}

fn freq_shift_inplace<T>(signal: &mut [Complex<T>], freq_hz: f64, sample_rate_hz: f64)
where
    T: Float + FromPrimitive,
{
    let two_pi = T::from(2.0 * std::f64::consts::PI).unwrap();
    let fs = T::from(sample_rate_hz).unwrap();
    let f = T::from(freq_hz).unwrap();

    for (i, sample) in signal.iter_mut().enumerate() {
        let t = T::from(i).unwrap() / fs;
        let phase = two_pi * f * t;
        let phasor = Complex::new(phase.cos(), phase.sin());
        *sample = *sample * phasor;
    }
}

pub fn estimate_noise_floor<T>(surface: &CafSurface<T>) -> T
where
    T: Float + Debug,
{
    let mut all_values: Vec<T> = surface
        .surface
        .iter()
        .flat_map(|row| row.iter().copied())
        .filter(|&val| val > T::zero()) // Only positive values
        .collect();

    if all_values.is_empty() {
        return T::from(1e-10).unwrap(); // Fallback for empty surface
    }

    // Sort to find median
    all_values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    let median_idx = all_values.len() / 2;
    all_values[median_idx]
}

pub fn find_peak<T>(surface: &CafSurface<T>) -> Peak
where
    T: Float + Debug,
{
    // Estimate noise floor using median
    let noise_floor = estimate_noise_floor(surface);

    let mut max_magnitude = T::zero();
    let mut max_doppler_idx = 0;
    let mut max_delay_idx = 0;

    for (doppler_idx, row) in surface.surface.iter().enumerate() {
        for (delay_idx, &magnitude) in row.iter().enumerate() {
            if magnitude > max_magnitude {
                max_magnitude = magnitude;
                max_doppler_idx = doppler_idx;
                max_delay_idx = delay_idx;
            }
        }
    }

    // Calculate SNR in dB
    let magnitude_f64 = max_magnitude.to_f64().unwrap();
    let noise_floor_f64 = noise_floor.to_f64().unwrap();
    let snr_db = 10.0 * (magnitude_f64 / noise_floor_f64).log10();

    Peak {
        magnitude: magnitude_f64,
        delay_idx: max_delay_idx,
        doppler_idx: max_doppler_idx,
        delay_samples: surface.time_delays[max_delay_i

pub fn interpolate_peak<T>(
    surface: &CafSurface<T>,
    coarse_peak: &Peak,
    method: InterpolationMethod,
) -> RefinedPeak
where
    T: Float + Debug,
{
    match method {
        InterpolationMethod::Parabolic2D => interpolate_parabolic_2d(surface, coarse_peak),
        InterpolationMethod::Sinc2D => interpolate_sinc_2d(surface, coarse_peak),
        InterpolationMethod::None => RefinedPeak {
            magnitude: coarse_peak.magnitude,
            delay_samples: coarse_peak.delay_samples,
            doppler_hz: coarse_peak.doppler_hz,
            snr_db: coarse_peak.snr_db,
            noise_floor: coarse_peak.noise_floor,
        },
    }
}

fn interpolate_parabolic_2d<T>(surface: &CafSurface<T>, peak: &Peak) -> RefinedPeak
where
    T: Float + Debug,
{
    let d_idx = peak.doppler_idx;
    let t_idx = peak.delay_idx;

    let n_doppler = surface.surface.len();
    let n_time = surface.surface[0].len();

    // Check if we have enough neighbors for interpolation
    if d_idx == 0 || d_idx >= n_doppler - 1 || t_idx == 0 || t_idx >= n_time - 1 {
        // Peak is on the edge, can't interpolate
        return RefinedPeak {
            magnitude: peak.magnitude,
            delay_samples: peak.delay_samples,
            doppler_hz: peak.doppler_hz,
            snr_db: peak.snr_db,
            noise_floor: peak.noise_floor,
        };
    }

    // Extract 3x3 neighborhood (doppler, time)
    let z00 = surface.surface[d_idx - 1][t_idx - 1].to_f64().unwrap();
    let z01 = surface.surface[d_idx - 1][t_idx].to_f64().unwrap();
    let z02 = surface.surface[d_idx - 1][t_idx + 1].to_f64().unwrap();
    let z10 = surface.surface[d_i

fn interpolate_sinc_2d<T>(surface: &CafSurface<T>, peak: &Peak) -> RefinedPeak
where
    T: Float + Debug,
{
    // Extract region around peak
    let d_idx = peak.doppler_idx;
    let t_idx = peak.delay_idx;

    let n_doppler = surface.surface.len();
    let n_time = surface.surface[0].len();

    // Define interpolation window size
    let window_size = 16.min(n_time / 2).min(n_doppler / 2);

    // Extract window around peak
    let d_start = d_idx.saturating_sub(window_size / 2);
    let d_end = (d_idx + window_size / 2).min(n_doppler);
    let t_start = t_idx.saturating_sub(window_size / 2);
    let t_end = (t_idx + window_size / 2).min(n_time);

    let d_size = d_end - d_start;
    let t_size = t_end - t_start;

    // Extract local surface
    let mut local_surface: Vec<Vec<f64>> = vec![vec![0.0; t_size]; d_size];
    for (i, d) in (d_start..d_end).enumerate() {
        for (j, t) in (t_start..t_end).enumerate() {
            local_surface[i][j] = surface.surface[d][t].to_f64(

fn upsample_2d(surface: &[Vec<f64>], factor: usize) -> Vec<Vec<f64>> {
    let n_doppler = surface.len();
    let n_time = surface[0].len();

    let new_n_doppler = n_doppler * factor;
    let new_n_time = n_time * factor;

    // For simplicity, upsample each dimension separately
    // First upsample time dimension
    let mut time_upsampled: Vec<Vec<f64>> = Vec::new();
    for row in surface.iter() {
        time_upsampled.push(upsample_1d(row, factor));
    }

    // Then upsample doppler dimension (transpose, upsample, transpose back)
    let mut doppler_upsampled: Vec<Vec<f64>> = vec![vec![0.0; new_n_time]; new_n_doppler];

    for t in 0..new_n_time {
        let column: Vec<f64> = time_upsampled.iter().map(|row| row[t]).collect();
        let upsampled_column = upsample_1d(&column, factor);
        for (d, &val) in upsampled_column.iter().enumerate() {
            doppler_upsampled[d][t] = val;
        }
    }

    doppler_upsampled
}

fn upsample_1d(signal: &[f64], factor: usize) -> Vec<f64> {
    let n = signal.len();
    let new_n = n * factor;

    // Convert to complex
    let mut signal_complex: Vec<Complex<f64>> = signal.iter().map(|&x| Complex::new(x, 0.0)).collect();

    // FFT
    fft(&mut signal_complex);

    // Zero-pad in frequency domain
    let mut padded = vec![Complex::new(0.0, 0.0); new_n];

    // Copy positive frequencies
    let half = n / 2;
    for i in 0..=half {
        padded[i] = signal_complex[i];
    }

    // Copy negative frequencies
    for i in 1..half {
        padded[new_n - i] = signal_complex[n - i];
    }

    // IFFT
    ifft(&mut padded);

    // Extract real part and scale
    padded.iter().map(|c| c.re * factor as f64).collect()
}

    fn circular_shift<T: Clone>(signal: &[T], shift: isize) -> Vec<T> {
        let n = signal.len();
        let shift = ((shift % n as isize) + n as isize) as usize % n;

        let mut shifted = signal.to_vec();
        shifted.rotate_right(shift);
        shifted
    }

    fn test_caf_autocorrelation() {
        if !should_plot() {
            println!("Skipping test_caf_autocorrelation plot (set PLOT=true to enable)");
            return;
        }

        let sample_rate = 1e6;
        let length = 1024;

        // Generate signal
        let signal = generate_test_signal(length, sample_rate, 42);

        // CAF parameters
        let params = CafParams {
            time_step: 1,
            freq_step_hz: 50.0,
            doppler_range_hz: (-500.0, 500.0),
            delay_range_samples: (-50.0, 50.0),
            sample_rate_hz: sample_rate,
        };

        // Compute CAF
        let surface = compute_caf(&signal, &signal, &params).unwrap();

        // Find peak
        let peak = find_peak(&surface);

        println!("Autocorrelation peak:");
        println!("  Delay: {} samples", peak.delay_samples);
        println!("  FDOA: {} Hz", peak.doppler_hz);
        println!("  Magnitude: {}", peak.magnitude);
        println!("  SNR: {:.1

    fn test_caf_time_delay() {
        let sample_rate = 1e6;
        let length = 1024;
        let time_delay = 50; // samples

        // Generate common signal
        let signal_common = generate_test_signal(length, sample_rate, 0);
        let noise1 = generate_test_signal(length, sample_rate, 1);
        let noise2 = generate_test_signal(length, sample_rate, 2);

        // signal1 = common + noise1 * 0.1
        let signal1 = &signal_common + &(&noise1 * 0.1);

        // signal2 = shifted(common) + noise2 * 0.1
        let signal_common_shifted = ComplexVec::from_vec(circular_shift(&signal_common, time_delay as isize));
        let signal2 = &signal_common_shifted + &(&noise2 * 0.1);

        // CAF parameters
        let params = CafParams {
            time_step: 1,
            freq_step_hz: 50.0,
            doppler_range_hz: (-500.0, 500.0),
            delay_range_samples: (0.0, 100.0),
            sample_rate_hz: sample_rate,
        };

        // Compute CAF
        le

    fn test_caf_frequency_shift() {
        if !should_plot() {
            println!("Skipping test_caf_frequency_shift plot (set PLOT=true to enable)");
            return;
        }

        let sample_rate = 1e6;
        let length = 1048576; //2**20
        let freq_shift_hz = 500.0;

        // Generate common signal
        let signal_common = generate_test_signal(length, sample_rate, 42);
        let noise1 = generate_test_signal(length, sample_rate, 100);
        let noise2 = generate_test_signal(length, sample_rate, 200);

        // signal1 = common + noise1 * 0.1
        let signal1 = &signal_common + &(&noise1 * 0.1);

        // signal2 = freq_shift(common) + noise2 * 0.1
        let mut signal_common_shifted = signal_common.clone();
        signal_common_shifted.freq_shift(freq_shift_hz, sample_rate);
        let signal2 = &signal_common_shifted + &(&noise2 * 0.1);

        // CAF parameters
        let params = CafParams {
            time_step: 1,
            freq_step_

    fn test_caf_tdoa_doppler_combined() {
        if !should_plot() {
            println!("Skipping test_caf_tdoa_doppler_combined plot (set PLOT=true to enable)");
            return;
        }

        let sample_rate = 1e6;
        let length = 1024;
        let time_delay = 40;
        let freq_shift_hz = 150.0;

        // Generate common signal
        let signal_common = generate_test_signal(length, sample_rate, 42);
        let noise1 = generate_test_signal(length, sample_rate, 100);
        let noise2 = generate_test_signal(length, sample_rate, 200);

        // signal1 = common + noise1 * 0.1
        let signal1 = &signal_common + &(&noise1 * 0.1);

        // signal2 = freq_shift(time_shift(common)) + noise2 * 0.1
        let mut signal_common_shifted = ComplexVec::from_vec(circular_shift(&signal_common, time_delay as isize));
        signal_common_shifted.freq_shift(freq_shift_hz, sample_rate);
        let signal2 = &signal_common_shifted + &(&noise2 * 0.1);

        // CA

    fn test_caf_interpolation() {
        if !should_plot() {
            println!("Skipping test_caf_interpolation plot (set PLOT=true to enable)");
            return;
        }

        let sample_rate = 1e6;
        let length = 2048;
        let time_delay = 37;
        let freq_shift_hz = 123.0;

        // Generate common signal
        let signal_common = generate_test_signal(length, sample_rate, 42);

        // Apply shifts
        let signal1 = signal_common.clone();
        let mut signal2 = ComplexVec::from_vec(circular_shift(&signal_common, time_delay as isize));
        signal2.freq_shift(freq_shift_hz, sample_rate);

        // CAF parameters (coarser grid to test interpolation)
        let params = CafParams {
            time_step: 2, // Coarse mode
            freq_step_hz: 50.0,
            doppler_range_hz: (-500.0, 500.0),
            delay_range_samples: (0.0, 100.0),
            sample_rate_hz: sample_rate,
        };

        // Compute CAF
        let surface 

    fn test_caf_invalid_doppler_range() {
        let sample_rate = 1e6;
        let length = 1024;

        let signal1 = generate_test_signal(length, sample_rate, 42);
        let signal2 = generate_test_signal(length, sample_rate, 43);

        // Invalid Doppler range (min >= max)
        let params = CafParams {
            time_step: 1,
            freq_step_hz: 50.0,
            doppler_range_hz: (500.0, -500.0), // min > max
            delay_range_samples: (0.0, 100.0),
            sample_rate_hz: sample_rate,
        };

        let result = compute_caf(&signal1, &signal2, &params);
        assert!(result.is_err());

        if let Err(CafError::InvalidDopplerRange { min, max }) = result {
            assert_eq!(min, 500.0);
            assert_eq!(max, -500.0);
        } else {
            panic!("Expected InvalidDopplerRange error");
        }
    }

    fn test_caf_invalid_delay_range() {
        let sample_rate = 1e6;
        let length = 1024;

        let signal1 = generate_test_signal(length, sample_rate, 42);
        let signal2 = generate_test_signal(length, sample_rate, 43);

        // Invalid delay range (min >= max)
        let params = CafParams {
            time_step: 1,
            freq_step_hz: 50.0,
            doppler_range_hz: (-500.0, 500.0),
            delay_range_samples: (100.0, 0.0), // min > max
            sample_rate_hz: sample_rate,
        };

        let result = compute_caf(&signal1, &signal2, &params);
        assert!(result.is_err());

        if let Err(CafError::InvalidDelayRange { min, max }) = result {
            assert_eq!(min, 100.0);
            assert_eq!(max, 0.0);
        } else {
            panic!("Expected InvalidDelayRange error");
        }
    }

    fn test_caf_signal_length_mismatch() {
        let sample_rate = 1e6;

        let signal1 = generate_test_signal(1024, sample_rate, 42);
        let signal2 = generate_test_signal(2048, sample_rate, 43); // Different length

        let params = CafParams {
            time_step: 1,
            freq_step_hz: 50.0,
            doppler_range_hz: (-500.0, 500.0),
            delay_range_samples: (0.0, 100.0),
            sample_rate_hz: sample_rate,
        };

        let result = compute_caf(&signal1, &signal2, &params);
        assert!(result.is_err());

        if let Err(CafError::SignalLengthMismatch { signal1_len, signal2_len }) = result {
            assert_eq!(signal1_len, 1024);
            assert_eq!(signal2_len, 2048);
        } else {
            panic!("Expected SignalLengthMismatch error");
        }
    }

    fn test_auto_compute_caf() {
        let sample_rate = 1e6;
        let bandwidth = 10e3; // 10 kHz bandwidth
        let length = 65536; // Longer signal for better frequency resolution
        let time_delay = 50; // samples
        let freq_shift_hz = 100.0;

        // Generate common signal
        let signal_common = generate_test_signal(length, sample_rate, 42);
        let noise1 = generate_test_signal(length, sample_rate, 100);
        let noise2 = generate_test_signal(length, sample_rate, 200);

        // signal1 = common + noise1 * 0.1
        let signal1 = &signal_common + &(&noise1 * 0.1);

        // signal2 = freq_shift(time_shift(common)) + noise2 * 0.1
        let mut signal_common_shifted = ComplexVec::from_vec(circular_shift(&signal_common, time_delay as isize));
        signal_common_shifted.freq_shift(freq_shift_hz, sample_rate);
        let signal2 = &signal_common_shifted + &(&noise2 * 0.1);

        // Use auto_compute_caf
        let surface = auto_compute

    pub fn plot_caf_surface_3d(surface: &CafSurface<f64>, peak: Option<&Peak>, title: &str) {
        let mut plot = Plot::new();

        // Estimate noise floor using median
        let noise_floor = estimate_noise_floor(surface);
        let noise_floor_db = 10.0 * noise_floor.log10();

        // Convert to dB relative to noise floor
        let surface_db: Vec<Vec<f64>> = surface
            .surface
            .iter()
            .map(|row| {
                row.iter()
                    .map(|&val| {
                        if val > 0.0 {
                            let val_db = 10.0 * val.log10();
                            // Relative to noise floor, clamped at -100 dB
                            (val_db - noise_floor_db).max(-100.0)
                        } else {
                            -100.0
                        }
                    })
                    .collect()
            })
            .collect();

        let trace = Surface::new(surface_db)
           

    pub fn plot_caf_heatmap(surface: &CafSurface<f64>, peak: Option<&Peak>, title: &str) {
        let mut plot = Plot::new();

        // Estimate noise floor using median
        let noise_floor = estimate_noise_floor(surface);
        let noise_floor_db = 10.0 * noise_floor.log10();

        // Convert to dB relative to noise floor
        let surface_db: Vec<Vec<f64>> = surface
            .surface
            .iter()
            .map(|row| {
                row.iter()
                    .map(|&val| {
                        if val > 0.0 {
                            let val_db = 10.0 * val.log10();
                            // Relative to noise floor, clamped at -100 dB
                            (val_db - noise_floor_db).max(-100.0)
                        } else {
                            -100.0
                        }
                    })
                    .collect()
            })
            .collect();

        let trace = HeatMap::new(
            surface.time_

    pub fn plot_caf_slices(surface: &CafSurface<f64>, peak: &Peak, title: &str) {
        use plotly::layout::{GridPattern, LayoutGrid, RowOrder};

        let mut plot = Plot::new();

        // Estimate noise floor using median
        let noise_floor = estimate_noise_floor(surface);
        let noise_floor_db = 10.0 * noise_floor.log10();

        // Time slice at peak Doppler (relative to noise floor)
        let time_slice: Vec<f64> = surface.surface[peak.doppler_idx]
            .iter()
            .map(|&val| {
                if val > 0.0 {
                    let val_db = 10.0 * val.log10();
                    (val_db - noise_floor_db).max(-100.0)
                } else {
                    -100.0
                }
            })
            .collect();

        let time_trace = Scatter::new(
            surface.time_delays.iter().map(|&x| x).collect(),
            time_slice,
        )
        .mode(Mode::Lines)
        .name("Time Slice")
        .x_axis("x1")
        .y_a

    pub fn normalize(&mut self) -> ComplexVec<T> {
        ComplexVec::from_vec(
            self.vector.iter()
            .map(|c| {
                let mag = c.norm();
                if mag > T::zero() {
                    *c / mag
                } else {
                    *c
                }
            }).collect()
        )
    }

    pub fn normalize_inplace(&mut self) {
        for c in self.vector.iter_mut() {
            let mag = c.norm();
            if mag > T::zero() {
                *c = *c / mag;
            }
        }
    }

    pub fn convolve(&self, kernel: &ComplexVec<T>, mode: ConvMode) -> ComplexVec<T> {
        let input_len = self.vector.len();
        let kernel_len = kernel.vector.len();

        // Compute convolution in full mode first
        let full_size = input_len + kernel_len - 1;
        let mut full_result = vec![Complex::new(T::zero(), T::zero()); full_size];

        // Full convolution computation
        // Note: This implements cross-correlation (not true convolution with kernel flip)
        // This matches the behavior of the original implementation
        for i in 0..full_size {
            let mut sum = Complex::new(T::zero(), T::zero());
            for j in 0..kernel_len {
                // Correlation: y[i] = sum(x[i+j] * h[j])  for valid indices
                // For full mode, we need to account for zero-padding
                if i + j >= kernel_len - 1 && i + j < input_len + kernel_len - 1 {
                    let signal_idx = i + j - (kernel_len - 1);
               

    pub fn freq_shift(&mut self, freq_offset: f64, sample_rate: f64) {
        assert!(sample_rate > 0.0, "sample_rate must be positive");

        let normalized_freq = freq_offset / sample_rate;
        let two_pi = T::from(std::f64::consts::PI).unwrap() * T::from(2.0).unwrap();

        for (n, sample) in self.vector.iter_mut().enumerate() {
            let phase = two_pi * T::from(normalized_freq).unwrap() * T::from(n as f64).unwrap();
            let phase_f64 = phase.to_f64().unwrap();

            // e^(j*phase) = cos(phase) + j*sin(phase)
            let cos_phase = T::from(phase_f64.cos()).unwrap();
            let sin_phase = T::from(phase_f64.sin()).unwrap();
            let rotator = Complex::new(cos_phase, sin_phase);

            *sample = *sample * rotator;
        }
    }

    pub fn measure_power(&self, oversample_rate: Option<f64>) -> f64 {
        let mut sum_power = T::from(0.0).unwrap();
        for sample in self.vector.iter() {
            let power = (sample.re * sample.re) + (sample.im * sample.im);
            sum_power = sum_power + power;
        }
        let factor = oversample_rate.unwrap_or(1.0_f64);
        sum_power.to_f64().unwrap() / ((self.vector.len() as f64) / factor)
    }

    pub fn scale_to_power(&self, target_power: f64, oversample_rate: Option<f64>) -> ComplexVec<T> {
        if target_power <= 0.0 {
            panic!("target_power must be positive");
        }

        let current_power = self.measure_power(oversample_rate);

        if current_power == 0.0 {
            // Signal is silent, return a copy as-is
            return ComplexVec::from_vec(self.vector.clone());
        }

        let scale = (target_power / current_power).sqrt();
        let scale_t = T::from(scale).unwrap();

        let result = self.vector.iter().map(|sample| sample * scale_t).collect();
        ComplexVec::from_vec(result)
    }

    fn add(self, other: ComplexVec<T>) -> Self {
        assert_eq!(
            self.vector.len(),
            other.vector.len(),
            "ComplexVec addition requires equal lengths"
        );
        let result = self
            .vector
            .iter()
            .zip(other.vector.iter())
            .map(|(a, b)| a + b)
            .collect();
        ComplexVec::from_vec(result)
    }

    fn add(self, other: &ComplexVec<T>) -> ComplexVec<T> {
        assert_eq!(
            self.vector.len(),
            other.vector.len(),
            "ComplexVec addition requires equal lengths"
        );
        let result = self
            .vector
            .iter()
            .zip(other.vector.iter())
            .map(|(a, b)| a + b)
            .collect();
        ComplexVec::from_vec(result)
    }

    fn sub(self, other: ComplexVec<T>) -> Self {
        assert_eq!(
            self.vector.len(),
            other.vector.len(),
            "ComplexVec subtraction requires equal lengths"
        );
        let result = self
            .vector
            .iter()
            .zip(other.vector.iter())
            .map(|(a, b)| a - b)
            .collect();
        ComplexVec::from_vec(result)
    }

    fn sub(self, other: &ComplexVec<T>) -> ComplexVec<T> {
        assert_eq!(
            self.vector.len(),
            other.vector.len(),
            "ComplexVec subtraction requires equal lengths"
        );
        let result = self
            .vector
            .iter()
            .zip(other.vector.iter())
            .map(|(a, b)| a - b)
            .collect();
        ComplexVec::from_vec(result)
    }

    fn mul(self, other: ComplexVec<T>) -> Self {
        assert_eq!(
            self.vector.len(),
            other.vector.len(),
            "ComplexVec element-wise multiplication requires equal lengths"
        );
        let result = self
            .vector
            .iter()
            .zip(other.vector.iter())
            .map(|(a, b)| a * b)
            .collect();
        ComplexVec::from_vec(result)
    }

    fn mul(self, other: &ComplexVec<T>) -> ComplexVec<T> {
        assert_eq!(
            self.vector.len(),
            other.vector.len(),
            "ComplexVec element-wise multiplication requires equal lengths"
        );
        let result = self
            .vector
            .iter()
            .zip(other.vector.iter())
            .map(|(a, b)| a * b)
            .collect();
        ComplexVec::from_vec(result)
    }

    fn test_replace_vec() {
        let mut cv = ComplexVec::new();
        let data = vec![Complex::new(1.0, 0.0), Complex::new(0.0, 1.0)];
        cv.replace_vec(data);
        assert_eq!(cv.len(), 2);
    }

    fn test_abs() {
        let data = vec![
            Complex::new(3.0, 4.0),  // mag = 5.0
            Complex::new(0.0, 1.0),  // mag = 1.0
        ];
        let mut cv = ComplexVec::from_vec(data);
        let magnitudes = cv.abs();

        assert!((magnitudes[0] - 5.0).abs() < 1e-10);
        assert!((magnitudes[1] - 1.0).abs() < 1e-10);
    }

    fn test_normalize() {
        let data = vec![
            Complex::new(3.0, 4.0),
            Complex::new(5.0, 12.0),
        ];
        let mut cv = ComplexVec::from_vec(data);
        let normalized = cv.normalize();
        let mags = normalized.vector;

        assert!((mags[0].norm() - 1.0).abs() < 1e-10);
        assert!((mags[1].norm() - 1.0).abs() < 1e-10);
    }

    fn test_normalize_inplace() {
        let data = vec![
            Complex::new(3.0, 4.0),
            Complex::new(5.0, 12.0),
        ];
        let mut cv = ComplexVec::from_vec(data);
        cv.normalize_inplace();

        assert!((cv.vector[0].norm() - 1.0).abs() < 1e-10);
        assert!((cv.vector[1].norm() - 1.0).abs() < 1e-10);
    }

    fn test_normalize_zero_magnitude() {
        let data = vec![
            Complex::new(0.0, 0.0),
            Complex::new(3.0, 4.0),
        ];
        let mut cv = ComplexVec::from_vec(data);
        cv.normalize_inplace();

        assert_eq!(cv.vector[0], Complex::new(0.0, 0.0));
        assert!((cv.vector[1].norm() - 1.0).abs() < 1e-10);
    }

    fn test_convolve_valid() {
        let signal = vec![
            Complex::new(1.0, 0.0),
            Complex::new(2.0, 0.0),
            Complex::new(3.0, 0.0),
            Complex::new(4.0, 0.0),
            Complex::new(5.0, 0.0),
        ];
        let kernel = vec![
            Complex::new(2.0, 0.0),
            Complex::new(1.0, 0.0),
        ];

        let sig = ComplexVec::from_vec(signal);
        let ker = ComplexVec::from_vec(kernel);
        let result = sig.convolve(&ker, ConvMode::Valid);

        assert_eq!(result.len(), 4); // 5 - 2 + 1 = 4
        assert_eq!(result.vector[0], Complex::new(4.0, 0.0));
        assert_eq!(result.vector[1], Complex::new(7.0, 0.0));
        assert_eq!(result.vector[2], Complex::new(10.0, 0.0));
        assert_eq!(result.vector[3], Complex::new(13.0, 0.0));
    }

    fn test_convolve_full() {
        let signal = vec![
            Complex::new(1.0, 0.0),
            Complex::new(2.0, 0.0),
            Complex::new(3.0, 0.0),
        ];
        let kernel = vec![
            Complex::new(1.0, 0.0),
            Complex::new(1.0, 0.0),
        ];

        let sig = ComplexVec::from_vec(signal);
        let ker = ComplexVec::from_vec(kernel);
        let result = sig.convolve(&ker, ConvMode::Full);

        assert_eq!(result.len(), 4); // 3 + 2 - 1 = 4
        assert_eq!(result.vector[0], Complex::new(1.0, 0.0));
        assert_eq!(result.vector[1], Complex::new(3.0, 0.0));
        assert_eq!(result.vector[2], Complex::new(5.0, 0.0));
        assert_eq!(result.vector[3], Complex::new(3.0, 0.0));
    }

    fn test_convolve_same() {
        let signal = vec![
            Complex::new(1.0, 0.0),
            Complex::new(2.0, 0.0),
            Complex::new(3.0, 0.0),
        ];
        let kernel = vec![
            Complex::new(1.0, 0.0),
            Complex::new(1.0, 0.0),
        ];

        let sig = ComplexVec::from_vec(signal);
        let ker = ComplexVec::from_vec(kernel);
        let result = sig.convolve(&ker, ConvMode::Same);

        assert_eq!(result.len(), 3); // Same as input
        // In same mode with kernel size 2, output is shifted by (2-1)/2 = 0
        assert_eq!(result.vector[0], Complex::new(1.0, 0.0));
        assert_eq!(result.vector[1], Complex::new(3.0, 0.0));
        assert_eq!(result.vector[2], Complex::new(5.0, 0.0));
    }

    fn test_convolve_inplace() {
        let signal = vec![
            Complex::new(1.0, 0.0),
            Complex::new(2.0, 0.0),
            Complex::new(3.0, 0.0),
        ];
        let kernel = vec![
            Complex::new(1.0, 0.0),
            Complex::new(1.0, 0.0),
        ];

        let mut sig = ComplexVec::from_vec(signal);
        let ker = ComplexVec::from_vec(kernel);
        sig.convolve_inplace(&ker, ConvMode::Valid);

        assert_eq!(sig.len(), 2);
        assert_eq!(sig.vector[0], Complex::new(3.0, 0.0));
        assert_eq!(sig.vector[1], Complex::new(5.0, 0.0));
    }

    fn test_convolve_impulse() {
        let signal = vec![
            Complex::new(1.0, 2.0),
            Complex::new(3.0, 4.0),
        ];
        let impulse = vec![Complex::new(1.0, 0.0)];

        let sig = ComplexVec::from_vec(signal.clone());
        let imp = ComplexVec::from_vec(impulse);
        let result = sig.convolve(&imp, ConvMode::Valid);

        for i in 0..signal.len() {
            println!("Signal: {}, Convolved Result: {}", signal[i], result[i]);
        }
        assert_eq!(result.len(), 2);
        assert_eq!(result.vector[0], signal[0]);
        assert_eq!(result.vector[1], signal[1]);
    }

    fn test_indexing() {
        let data = vec![
            Complex::new(1.0, 2.0),
            Complex::new(3.0, 4.0),
            Complex::new(5.0, 6.0),
        ];
        let mut cv = ComplexVec::from_vec(data);

        // Test read indexing
        assert_eq!(cv[0], Complex::new(1.0, 2.0));
        assert_eq!(cv[1], Complex::new(3.0, 4.0));
        assert_eq!(cv[2], Complex::new(5.0, 6.0));

        // Test write indexing
        cv[1] = Complex::new(7.0, 8.0);
        assert_eq!(cv[1], Complex::new(7.0, 8.0));
    }

    fn test_mul_element_wise() {
        let a = ComplexVec::from_vec(vec![
            Complex::new(2.0, 1.0),
            Complex::new(3.0, 2.0),
        ]);
        let b = ComplexVec::from_vec(vec![
            Complex::new(1.0, 1.0),
            Complex::new(2.0, 0.0),
        ]);
        
        // Test owned multiplication
        let c = a.clone() * b.clone();
        // (2+i)(1+i) = 2 + 2i + i + i^2 = 2 + 3i - 1 = 1 + 3i
        assert_eq!(c[0], Complex::new(1.0, 3.0));
        // (3+2i)(2+0i) = 6 + 4i
        assert_eq!(c[1], Complex::new(6.0, 4.0));
        
        // Test borrowed multiplication
        let d = &a * &b;
        assert_eq!(d[0], Complex::new(1.0, 3.0));
        assert_eq!(d[1], Complex::new(6.0, 4.0));
    }

    fn test_mul_complex_scalar() {
        let vec = ComplexVec::from_vec(vec![
            Complex::new(1.0, 2.0),
            Complex::new(3.0, 4.0),
        ]);
        let scalar = Complex::new(2.0, 1.0);
        
        // Test vec * scalar (owned)
        let result1 = vec.clone() * scalar;
        // (1+2i)(2+i) = 2 + i + 4i + 2i^2 = 2 + 5i - 2 = 0 + 5i
        assert_eq!(result1[0], Complex::new(0.0, 5.0));
        // (3+4i)(2+i) = 6 + 3i + 8i + 4i^2 = 6 + 11i - 4 = 2 + 11i
        assert_eq!(result1[1], Complex::new(2.0, 11.0));
        
        // Test &vec * scalar (borrowed)
        let result2 = &vec * scalar;
        assert_eq!(result2[0], Complex::new(0.0, 5.0));
        assert_eq!(result2[1], Complex::new(2.0, 11.0));
        
        // Test scalar * vec (commutative, owned)
        let result3 = scalar * vec.clone();
        assert_eq!(result3[0], Complex::new(0.0, 5.0));
        assert_eq!(result3[1], Complex::new(2.0, 11.0));
        
        // Test scalar * &vec