Function bodies 301 total

    fn test_mul_real_scalar_f64() {
        let vec = ComplexVec::from_vec(vec![
            Complex::new(1.0_f64, 2.0_f64),
            Complex::new(3.0_f64, 4.0_f64),
        ]);
        let scalar = 2.5_f64;
        
        // Test vec * scalar (owned)
        let result1 = vec.clone() * scalar;
        assert_eq!(result1[0], Complex::new(2.5, 5.0));
        assert_eq!(result1[1], Complex::new(7.5, 10.0));
        
        // Test &vec * scalar (borrowed)
        let result2 = &vec * scalar;
        assert_eq!(result2[0], Complex::new(2.5, 5.0));
        assert_eq!(result2[1], Complex::new(7.5, 10.0));
        
        // Test scalar * vec (commutative, owned)
        let result3 = scalar * vec.clone();
        assert_eq!(result3[0], Complex::new(2.5, 5.0));
        assert_eq!(result3[1], Complex::new(7.5, 10.0));
        
        // Test scalar * &vec (commutative, borrowed)
        let result4 = scalar * &vec;
        assert_eq!(result4[0], Complex::new(2.5, 5.0));
        assert

    fn test_mul_real_scalar_f32() {
        let vec = ComplexVec::from_vec(vec![
            Complex::new(1.0_f32, 2.0_f32),
            Complex::new(3.0_f32, 4.0_f32),
        ]);
        let scalar = 2.5_f32;
        
        // Test vec * scalar (owned)
        let result1 = vec.clone() * scalar;
        assert_eq!(result1[0], Complex::new(2.5, 5.0));
        assert_eq!(result1[1], Complex::new(7.5, 10.0));
        
        // Test &vec * scalar (borrowed)
        let result2 = &vec * scalar;
        assert_eq!(result2[0], Complex::new(2.5, 5.0));
        assert_eq!(result2[1], Complex::new(7.5, 10.0));
        
        // Test scalar * vec (commutative, owned)
        let result3 = scalar * vec.clone();
        assert_eq!(result3[0], Complex::new(2.5, 5.0));
        assert_eq!(result3[1], Complex::new(7.5, 10.0));
        
        // Test scalar * &vec (commutative, borrowed)
        let result4 = scalar * &vec;
        assert_eq!(result4[0], Complex::new(2.5, 5.0));
        assert

pub fn fft<T>(input : &mut [Complex<T>])
where
    T: Float + RemAssign + DivAssign + Send + Sync + FromPrimitive + Signed + Debug + 'static,
{
    let mut planner = FftPlanner::<T>::new();
    let fft_forward = planner.plan_fft_forward(input.len());
    fft_forward.process(input);
    scale::<T>(input);
}

pub fn ifft<T>(input: &mut [Complex<T>])
where
    T: Float + Send + Sync + FromPrimitive + Signed + Debug + 'static,
{
    let mut planner = FftPlanner::<T>::new();
    let fft_inverse = planner.plan_fft_inverse(input.len());
    fft_inverse.process(input);
}

pub fn scale<T>(input: &mut [Complex<T>])
where
    T: Float + RemAssign + DivAssign,
{
    let nfft = T::from(input.len()).unwrap();
    let scale_val = Complex::<T>::new(
        T::from(1.0).unwrap() / nfft,
        T::from(0.0).unwrap()
    );
    input.iter_mut().for_each(|x| *x = *x * scale_val);
}

pub fn fftshift<T>(input_vec: &mut [T]) {
    let n = input_vec.len();
    let mid = (n + 1) / 2;
    input_vec.rotate_left(mid);
}

pub fn fftfreqs<T: Float>(start: T, stop: T, num_points: usize) -> Vec<T> {
    let mut v = Vec::with_capacity(num_points);
    let step: T = (stop - start) / (T::from(num_points).unwrap() - T::from(1.0).unwrap());
    for i in 0..num_points {
        v.push(start + (T::from(i).unwrap() * step));
    }
    v
}

    fn test_fft() {
        let blocksize: usize = 2000;
        let freq_hz = 100.0f64;
        let sample_rate_hz = 1000.0f64;

        let mut cw_generator = CW::new(freq_hz, sample_rate_hz, blocksize);
        let mut cw: Vec<Complex<f32>> = cw_generator.generate_block::<f32>();
        fft::<f32>(&mut cw);
        let mut cw_fft = ComplexVec::from_vec(cw);
        let mut cw_fft_abs: Vec<f32> = cw_fft.abs();
        fftshift::<f32>(&mut cw_fft_abs);

        let freqs: Vec<f32> = fftfreqs::<f32>((-sample_rate_hz/2_f64) as f32, (sample_rate_hz/2_f64) as f32, blocksize);

        if should_plot() {
            use crate::plot::plot_spectrum;
            plot_spectrum(&freqs, &cw_fft_abs, "CW Signal Spectrum");
        }

        let (max_ind, max_val) = vector_ops::max(&cw_fft_abs);
        let max_freq = freqs[max_ind].clone();
        println!("Max (ind,val) is: {max_freq},{max_val}");
        assert_eq!(((max_freq as f64) - freq_hz).abs() < 1.0, true);
    }

    fn test_ifft() {
        let signal_original = vec![
            Complex::new(1.0, 1.0),
            Complex::new(2.0, 2.0),
            Complex::new(3.0, 3.0),
            Complex::new(4.0, 4.0),
            Complex::new(5.0, 5.0),
            Complex::new(6.0, 6.0),
            Complex::new(7.0, 7.0),
            Complex::new(8.0, 8.0)
        ];
        let mut signal_modified = signal_original.clone();
        fft(&mut signal_modified);
        ifft(&mut signal_modified);

        for i in 0..signal_original.len() {
            println!("i: {i} -- Original value: {} vs Modified value {}", signal_original[i], signal_modified[i]);
            assert_near::<f32>(signal_original[i].norm(), signal_modified[i].norm(), 1e-3);
        }
    }

    fn test_fftshift_even() {
        let vec_initial = vec![0, 1, 2, 3, 4, 5, 6, 7];
        let mut vec_shifted = vec_initial.clone();
        fftshift(&mut vec_shifted);
        assert_eq!(vec_shifted, vec![4, 5, 6, 7, 0, 1, 2, 3]);
    }

    fn test_fftshift_odd() {
        let vec_initial = vec![0, 1, 2, 3, 4, 5, 6, 7, 8];
        let mut vec_shifted = vec_initial.clone();
        fftshift(&mut vec_shifted);
        assert_eq!(vec_shifted, vec![5, 6, 7, 8, 0, 1, 2, 3, 4]);
    }

    fn test_fftshift_cw_even() {
        let vec_cw = vec![Complex::new(1.0_f32, 0.0_f32); 10];
        let mut vec_cv = ComplexVec::<f32>::from_vec(vec_cw);
        fft(&mut vec_cv);
        fftshift(&mut vec_cv);
        vec_cv.iter().for_each(|c| println!("{c}"));
        assert_eq!(vec_cv[5].norm(), 1.0);
    }

    fn test_fftshift_cw_odd() {
        let vec_cw = vec![Complex::new(1.0_f32, 0.0_f32); 9];
        let mut vec_cv = ComplexVec::<f32>::from_vec(vec_cw);
        fft(&mut vec_cv);
        fftshift(&mut vec_cv);
        vec_cv.iter().for_each(|c| println!("{c}"));
        assert_eq!(vec_cv[4].norm(), 1.0);
    }

    fn test_fftfreq() {
        let start = -100.0;
        let stop = 100.0;
        let num_points = 1000;
        let freqs = fftfreqs::<f64>(start, stop, num_points);
        println!("start: {start}, stop: {stop}, num_points: {num_points}");
        assert_eq!(*freqs.first().unwrap(), start);
        assert_eq!(*freqs.last().unwrap(), stop);
    }

pub fn create_butterworth_filter<T: Float>(num_samples: usize, order: i32, cutoff: T) -> Vec<T> {
    let mut filter_response = vec![T::zero(); num_samples];
    let one = T::one();

    for i in 0..num_samples {
        // Convert index to normalized frequency [-0.5, 0.5]
        let f = if i < num_samples / 2 {
            T::from(i).unwrap() / T::from(num_samples).unwrap()
        } else {
            (T::from(i).unwrap() - T::from(num_samples).unwrap()) / T::from(num_samples).unwrap()
        };

        let norm_freq = f.abs();

        // Butterworth magnitude response: H(f) = 1 / sqrt(1 + (f/fc)^(2n))
        let ratio = norm_freq / cutoff;
        let power = ratio.powi(2 * order);
        filter_response[i] = one / (one + power).sqrt();
    }

    filter_response
}

pub fn apply_butterworth_filter<T>(signal: &mut [Complex<T>], order: i32, cutoff: T)
where
    T: Float + Signed + Send + Sync + std::fmt::Debug + num_traits::cast::FromPrimitive + std::ops::RemAssign + std::ops::DivAssign + 'static,
{
    let n = signal.len();

    // Create filter in frequency domain
    let filter = create_butterworth_filter(n, order, cutoff);

    // FFT to frequency domain
    let mut planner = FftPlanner::new();
    let fft = planner.plan_fft_forward(n);
    fft.process(signal);

    // Apply filter
    for (i, sample) in signal.iter_mut().enumerate() {
        *sample = *sample * filter[i];
    }

    // IFFT back to time domain
    let ifft = planner.plan_fft_inverse(n);
    ifft.process(signal);

    // Normalize (rustfft doesn't normalize)
    let scale = T::one() / T::from(n).unwrap();
    for sample in signal.iter_mut() {
        *sample = *sample * scale;
    }
}

    fn generate_iq_noise(n_samples: usize, noise_floor_db: f64) -> Vec<Complex<f64>> {
        use rand::thread_rng;
        use rand_distr::{Distribution, Normal};

        let mut rng = thread_rng();
        let noise_power = 10_f64.powf(noise_floor_db / 10.0);
        let sigma = (noise_power / 2.0).sqrt();
        let normal = Normal::new(0.0, sigma).unwrap();

        (0..n_samples)
            .map(|_| Complex::new(normal.sample(&mut rng), normal.sample(&mut rng)))
            .collect()
    }

    fn test_filter_creation() {
        let filter = create_butterworth_filter(1024, 3, 0.45);
        assert_eq!(filter.len(), 1024);

        // DC should be 1.0 (unity gain)
        assert!((filter[0] - 1.0).abs() < 1e-10);

        // Frequencies beyond cutoff should be attenuated
        // Bin 512 is Nyquist (0.5), cutoff is 0.45, so it should be somewhat attenuated
        // For 3rd order Butterworth: H(0.5/0.45) = 1/sqrt(1 + (0.5/0.45)^6) ≈ 0.64
        assert!(filter[512] < 0.7);

        // Very high frequency (bin 256 = 0.25, which is < cutoff) should pass more
        assert!(filter[256] > 0.9);
    }

    fn test_dc_response() {
        // Create DC signal
        let mut signal: Vec<Complex<f64>> = vec![Complex::new(1.0, 1.0); 1024];

        // Apply filter
        apply_butterworth_filter(&mut signal, 3, 0.45);

        // DC should pass through with unity gain
        for sample in signal.iter() {
            assert!((sample.re - 1.0).abs() < 0.01);
            assert!((sample.im - 1.0).abs() < 0.01);
        }
    }

    fn test_butterworth_spectrum() {
        use crate::spectrum::welch::welch;
        use crate::spectrum::window::WindowType;
        use crate::vector_ops;
        use crate::plot::plot_spectrum;
        use crate::test_utils::should_plot;

        if !should_plot() {
            println!("Skipping Butterworth filter spectrum plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Butterworth Filter Frequency Response (FFT-based) ===\n");

        let sample_rate = 1e6;
        let n_samples = (2.0_f64).powf(20.0) as usize;

        // Generate white noise
        let noise = generate_iq_noise(n_samples, -80.0);

        // Apply 3rd order Butterworth filter
        println!("Filtering with 3rd order Butterworth (cutoff = 0.45*Nyquist = 225 kHz)");
        let mut filtered_signal = noise.clone();
        apply_butterworth_filter(&mut filtered_signal, 4, 0.55);

        // Compute Welch PSD of filtered signal
        let (freqs, psd) = welch(
       

    fn test_filter_response() {
        use crate::spectrum::welch::welch;
        use crate::spectrum::window::WindowType;
        use crate::vector_ops;
        use crate::plot::plot_spectrum;
        use crate::test_utils::should_plot;

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

        println!("\n=== Butterworth Filter Response: Before vs After ===\n");

        let sample_rate = 1e6;
        let n_samples = (2.0_f64).powf(20.0) as usize;

        // Generate white noise
        println!("Generating AWGN...");
        let noise = generate_iq_noise(n_samples, -80.0);

        // Compute PSD of UNFILTERED noise
        let (freqs, psd_unfiltered) = welch(
            &noise,
            sample_rate,
            2048,
            None,
            None,
            WindowType::Hann,
            None,
        );
        let psd_unfiltered_db: Vec<f64> = vector_ops::to_db(&psd_unfiltered);

    fn test_impulse_response() {
        use crate::plot::plot_spectrum;
        use crate::test_utils::should_plot;

        if !should_plot() {
            println!("Skipping impulse response test (set PLOT=true to enable)");
            return;
        }

        println!("\n=== Butterworth Filter Impulse Response ===\n");

        let sample_rate = 1e6;
        let n_samples = 4096;

        // Create impulse signal (1.0 at t=0, 0.0 elsewhere)
        let mut impulse: Vec<Complex<f64>> = vec![Complex::new(0.0, 0.0); n_samples];
        impulse[0] = Complex::new(1.0, 0.0);

        // Filter the impulse
        apply_butterworth_filter(&mut impulse, 3, 0.45);

        // Compute FFT to get frequency response
        use rustfft::FftPlanner;
        let mut planner = FftPlanner::new();
        let fft = planner.plan_fft_forward(n_samples);
        let mut freq_response = impulse.clone();
        fft.process(&mut freq_response);

        // Compute magnitude in dB and shift
        le

    fn test_dc_suppression_diagnostic() {
        use crate::spectrum::welch::welch;
        use crate::spectrum::window::WindowType;
        use crate::vector_ops;
        use crate::plot::plot_spectrum;
        use crate::fft::{fft, fftshift, fftfreqs};
        use crate::test_utils::should_plot;

        if !should_plot() {
            println!("Skipping DC suppression diagnostic (set PLOT=true to enable)");
            return;
        }

        println!("\n=== DC SUPPRESSION DIAGNOSTIC ===");
        println!("Testing 3 methods to determine source of DC drop:\n");

        let sample_rate = 1e6;
        let n_samples = (2.0_f64).powf(18.0) as usize; // 262144 samples
        let order = 3;
        let cutoff = 0.45;

        // METHOD 1: Direct filter coefficients (GROUND TRUTH)
        println!("METHOD 1: Direct Filter Coefficients");
        let filter = create_butterworth_filter(n_samples, order, cutoff);
        println!("  Filter[0] (DC bin) = {:.15}", filter[0]);
        pri

    fn test_dc_component_in_noise() {
        use crate::generate::awgn::AWGN;

        println!("\n=== CHECKING DC COMPONENT IN AWGN ===\n");

        let n_samples = 262144;
        let sample_rate = 1e6;

        let mut awgn = AWGN::new_from_seed(sample_rate, n_samples, 1.0, 42);
        let noise: Vec<Complex<f64>> = awgn.generate_block().to_vec();

        // Compute DC (mean) of the noise
        let sum: Complex<f64> = noise.iter().sum();
        let dc_mean = sum / (n_samples as f64);

        println!("Noise samples: {}", n_samples);
        println!("DC (mean) = {} + {}i", dc_mean.re, dc_mean.im);
        println!("DC magnitude = {}", dc_mean.norm());
        println!("DC power (magnitude^2) = {}", dc_mean.norm_sqr());
        println!("DC in dB = {} dB\n", 10.0 * dc_mean.norm_sqr().log10());

        println!("✓ AWGN has near-zero DC by design (zero mean Gaussian)!");
        println!("✓ This is why you see DC drop in spectrograms of filtered noise");
        println!("✓ Th

    fn test_butterworth_preserves_dc() {
        use crate::generate::awgn::AWGN;

        println!("\n=== PROOF: Butterworth Preserves DC ===\n");

        let n_samples = 262144;
        let sample_rate = 1e6;

        // Generate noise and ADD a DC offset
        let mut awgn = AWGN::new_from_seed(sample_rate, n_samples, 1.0, 42);
        let mut signal_with_dc: Vec<Complex<f64>> = awgn.generate_block().to_vec();

        // Add DC offset
        let dc_offset = Complex::new(10.0, 5.0);
        for sample in signal_with_dc.iter_mut() {
            *sample += dc_offset;
        }

        // Measure DC before filtering
        let dc_before: Complex<f64> = signal_with_dc.iter().sum::<Complex<f64>>() / (n_samples as f64);
        println!("DC BEFORE filtering: {} + {}i", dc_before.re, dc_before.im);
        println!("DC magnitude BEFORE: {}", dc_before.norm());
        println!("DC in dB BEFORE: {} dB\n", 20.0 * dc_before.norm().log10());

        // Apply Butterworth filter
        a

pub fn create_cosine_taper_filter(num_samples: usize, passband_end: f64, stopband_start: f64) -> Vec<f64> {
    let mut filter_response = vec![1.0; num_samples];
    
    for i in 0..num_samples {
        // Convert index to normalized frequency [-0.5, 0.5]
        let f = if i < num_samples / 2 {
            i as f64 / num_samples as f64
        } else {
            (i as f64 - num_samples as f64) / num_samples as f64
        };
        
        let norm_freq = f.abs();
        
        if norm_freq <= passband_end {
            // Passband: unity gain (including DC at i=0)
            filter_response[i] = 1.0;
        } else if norm_freq >= stopband_start {
            // Stopband: zero gain
            filter_response[i] = 0.0;
        } else {
            // Transition: raised cosine taper
            // Smoothly transitions from 1.0 to 0.0 using cosine
            let transition_width = stopband_start - passband_end;
            let position = (norm_freq - passband_end) / transiti

pub fn apply_cosine_taper_filter(signal: &mut [Complex<f64>], passband_end: f64, stopband_start: f64) {
    let n = signal.len();
    let filter = create_cosine_taper_filter(n, passband_end, stopband_start);
    
    // FFT to frequency domain
    let mut planner = FftPlanner::new();
    let fft = planner.plan_fft_forward(n);
    fft.process(signal);
    
    // Apply filter in frequency domain
    for (i, sample) in signal.iter_mut().enumerate() {
        *sample = *sample * filter[i];
    }
    
    // IFFT back to time domain
    let ifft = planner.plan_fft_inverse(n);
    ifft.process(signal);
    
    // Normalize (rustfft doesn't normalize IFFT)
    let scale = 1.0 / (n as f64);
    for sample in signal.iter_mut() {
        *sample = *sample * scale;
    }
}

    fn test_filter_creation() {
        let filter = create_cosine_taper_filter(1024, 0.40, 0.48);
        assert_eq!(filter.len(), 1024);

        // DC should be 1.0 (unity gain, NO suppression)
        assert_eq!(filter[0], 1.0);

        // Low frequencies in passband should be 1.0
        for i in 1..100 {
            let f = i as f64 / 1024.0;
            if f < 0.40 {
                assert_eq!(filter[i], 1.0, "Passband bin {} should be 1.0", i);
            }
        }

        // High frequencies in stopband should be 0.0
        let stopband_bin = (0.48 * 1024.0) as usize;
        assert!(filter[stopband_bin] < 0.01, "Stopband should be near 0");
    }

    fn test_dc_preserved() {
        // Create DC signal
        let mut signal: Vec<Complex<f64>> = vec![Complex::new(1.0, 1.0); 1024];

        // Apply filter
        apply_cosine_taper_filter(&mut signal, 0.40, 0.48);

        // DC should pass through with unity gain
        for sample in signal.iter() {
            assert!((sample.re - 1.0).abs() < 0.01, "Real part should be preserved");
            assert!((sample.im - 1.0).abs() < 0.01, "Imag part should be preserved");
        }
    }

    fn test_transition_band() {
        let filter = create_cosine_taper_filter(10000, 0.40, 0.48);
        
        // Check that transition is smooth (monotonic decreasing)
        let start_bin = (0.40 * 10000.0) as usize;
        let end_bin = (0.48 * 10000.0) as usize;
        
        for i in start_bin..(end_bin - 1) {
            assert!(filter[i] >= filter[i + 1], 
                "Transition should be monotonic: bin {} = {}, bin {} = {}", 
                i, filter[i], i+1, filter[i+1]);
        }
    }

pub fn resample<T>(signal: &ComplexVec<T>, factor: f64) -> ComplexVec<T>
where
    T: Float + RemAssign + DivAssign + Send + Sync + FromPrimitive + Signed + Debug + 'static,
{
    let input_len = signal.len();
    let output_len = ((input_len as f64) * factor).round() as usize;

    // Handle edge cases
    if output_len == 0 {
        return ComplexVec::new();
    }
    if output_len == input_len {
        return ComplexVec::from_vec(signal.iter().cloned().collect());
    }

    if output_len > input_len {
        upsample(signal, output_len)
    } else {
        downsample(signal, output_len)
    }
}

fn upsample<T>(signal: &ComplexVec<T>, output_len: usize) -> ComplexVec<T>
where
    T: Float + RemAssign + DivAssign + Send + Sync + FromPrimitive + Signed + Debug + 'static,
{
    let input_len = signal.len();

    // Clone input data for FFT processing
    let mut freq_data: Vec<Complex<T>> = signal.iter().cloned().collect();

    // Forward FFT (already scaled by 1/N)
    fft(&mut freq_data[..]);

    // Create zero-padded frequency domain data
    let mut padded_freq = vec![Complex::new(T::zero(), T::zero()); output_len];

    if input_len % 2 == 0 {
        // Even length: split Nyquist bin
        let half_input = input_len / 2;
        let nyquist_bin = freq_data[half_input];

        // Copy DC and positive frequencies
        for i in 0..=half_input {
            padded_freq[i] = freq_data[i];
        }

        // Split Nyquist bin (multiply by 0.5 and place at both ends)
        padded_freq[half_input] = nyquist_bin * T::from(0.5).unwrap();
        padded_freq[output_len - 

fn downsample<T>(signal: &ComplexVec<T>, output_len: usize) -> ComplexVec<T>
where
    T: Float + RemAssign + DivAssign + Send + Sync + FromPrimitive + Signed + Debug + 'static,
{
    let input_len = signal.len();

    // Clone input data for FFT processing
    let mut freq_data: Vec<Complex<T>> = signal.iter().cloned().collect();

    // Forward FFT (already scaled by 1/N)
    fft(&mut freq_data[..]);

    // Create truncated frequency domain data
    let mut truncated_freq = vec![Complex::new(T::zero(), T::zero()); output_len];

    if output_len % 2 == 0 {
        // Even output length: include Nyquist bin
        let half_output = output_len / 2;

        // Copy DC and positive frequencies
        for i in 0..=half_output {
            truncated_freq[i] = freq_data[i];
        }

        // Copy negative frequencies
        for i in 1..half_output {
            truncated_freq[output_len - i] = freq_data[input_len - i];
        }
    } else {
        // Odd output length: no Nyquis

    fn test_upsample_by_2() {
        // Create a simple DC signal
        let input: Vec<Complex<f64>> = vec![
            Complex::new(1.0, 0.0),
            Complex::new(1.0, 0.0),
            Complex::new(1.0, 0.0),
            Complex::new(1.0, 0.0),
        ];
        let signal = ComplexVec::from_vec(input.clone());

        // Upsample by factor of 2
        let resampled = resample(&signal, 2.0);

        // Check output length
        assert_eq!(resampled.len(), 8);

        // For DC signal, all values should remain approximately 1.0
        for i in 0..resampled.len() {
            assert!((resampled[i].re - 1.0).abs() < 1e-10,
                "Sample {} has value {}, expected 1.0", i, resampled[i].re);
            assert!(resampled[i].im.abs() < 1e-10,
                "Sample {} has imaginary part {}, expected 0.0", i, resampled[i].im);
        }
    }

    fn test_downsample_by_2() {
        // Create a simple DC signal
        let input: Vec<Complex<f64>> = vec![
            Complex::new(2.0, 0.0),
            Complex::new(2.0, 0.0),
            Complex::new(2.0, 0.0),
            Complex::new(2.0, 0.0),
            Complex::new(2.0, 0.0),
            Complex::new(2.0, 0.0),
            Complex::new(2.0, 0.0),
            Complex::new(2.0, 0.0),
        ];
        let signal = ComplexVec::from_vec(input.clone());

        // Downsample by factor of 0.5
        let resampled = resample(&signal, 0.5);

        // Check output length
        assert_eq!(resampled.len(), 4);

        // For DC signal, all values should remain approximately 2.0
        for i in 0..resampled.len() {
            assert!((resampled[i].re - 2.0).abs() < 1e-10,
                "Sample {} has value {}, expected 2.0", i, resampled[i].re);
            assert!(resampled[i].im.abs() < 1e-10,
                "Sample {} has imaginary part {}, expected 0.0", i, resamp

    fn test_bpsk_upsample_spectrum() {
        use crate::generate::psk_carrier::PskCarrier;
        use crate::mod_type::ModType;
        use crate::fft::{fft, fftshift, fftfreqs};
        use crate::vector_ops;
        use crate::test_utils::should_plot;

        if !should_plot() {
            println!("Skipping BPSK upsample spectrum plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== BPSK Upsample Spectrum Test ===");

        // Create BPSK carrier parameters
        let sample_rate_hz = 10e6_f64;
        let symbol_rate_hz = 1e6_f64;
        let rolloff_factor = 0.35_f64;
        let block_size = 4096;
        let filter_taps = 51;
        let upsample_rate = 1.5;

        // Generate BPSK carrier signal
        let mut carrier = PskCarrier::new(
            sample_rate_hz,
            symbol_rate_hz,
            ModType::_BPSK,
            rolloff_factor,
            block_size,
            filter_taps,
            Some(42), // seed for repr

    pub fn new(num_filter_taps: usize, sample_rate: f64, symbol_rate: f64, beta: f64) -> Self {
        let sps: f64 = sample_rate / symbol_rate;
        let scale = 1.0f64 / sample_rate;

        RRCFilter {
            num_filter_taps,
            beta,
            sps,
            symbol_rate,
            scale
        }
    }

    pub fn build_filter<T: Float>(&self) -> ComplexVec<T> {
        let pi = std::f64::consts::PI;
        let zero = 0.0f64;
        let one = 1.0f64;
        let two = 2.0f64;
        let four = 4.0f64;

        let mut taps: Vec<f64> = (0..self.num_filter_taps).map(|i| {
            // Use integer center to ensure t=0 at center tap
            let center = (self.num_filter_taps - 1) as f64 / two;
            let t = (i as f64 - center) * self.scale;

            // Normalize time by symbol period: t_norm = t * symbol_rate = t / T
            let t_norm = t * self.symbol_rate;

            if t == zero {
                // At t=0: h(0) = (1 + β*(4/π - 1)) / T
                (one + self.beta * (four / pi - one)) * self.symbol_rate
            }
            else if (four * self.beta * t_norm).abs() == one {
                // At discontinuity: t = ±T/(4β)
                let sqrt2 = two.sqrt();
                (self.beta / sqrt2) *
                ((one + two / pi) * (pi / (four * se

    fn test_build_filter_f32() {
        let filter32 = RRCFilter::new(101, 1e6f64, 100e3f64, 0.35f64);
        let taps32 = filter32.build_filter::<f32>();
        assert_eq!(taps32.len(), 101);
        println!("f32 Filter taps count: {}", taps32.len());
    }

    fn test_build_filter_f64() {
        let filter64 = RRCFilter::new(51, 1e6_f64, 500e3_f64, 0.25_f64);
        let taps64 = filter64.build_filter::<f64>();
        assert_eq!(taps64.len(), 51);
        println!("f64 Filter taps count: {}", taps64.len());
    }

    fn test_build_filter_complex32() {
        let filter32 = RRCFilter::new(101, 1.5, 1.0, 0.20);
        let taps32 = filter32.build_filter::<f32>();
        assert_eq!(taps32.len(), 101);

        println!("Complex32 Filter taps count: {}", taps32.len());

        // Verify all imaginary parts are near zero (filter taps should be real)
        for i in 0..5 {
            println!("tap[{}] = {} + {}i", i, taps32[i].re, taps32[i].im);
            assert!(taps32[i].im.abs() < 1e-10);
        }
    }

    pub fn new_from_seed(
        sample_rate_hz: f64,
        block_size: usize,
        noise_power: f64,
        seed: u64,
    ) -> Self {
        AWGN {
            sample_rate_hz,
            block_size,
            noise_power,
            rng: StdRng::seed_from_u64(seed),
        }
    }

    pub fn new_from_entropy(
        sample_rate_hz: f64,
        block_size: usize,
        noise_power: f64,
    ) -> Self {
        AWGN {
            sample_rate_hz,
            block_size,
            noise_power,
            rng: StdRng::from_entropy(),
        }
    }

    pub fn generate_block<T: Float>(&mut self) -> ComplexVec<T> {
        let mut samples = Vec::with_capacity(self.block_size);

        // Standard deviation per component for desired total power
        // Total power = I^2 + Q^2 = sigma^2 + sigma^2 = 2*sigma^2 = noise_power
        // Therefore: sigma = sqrt(noise_power/2)
        let std_dev = (self.noise_power / 2.0).sqrt();
        let std_dev_t = T::from(std_dev).unwrap();

        // Generate block_size complex samples
        for _ in 0..self.block_size {
            // Generate I and Q components from standard normal
            let i_std: f64 = StandardNormal.sample(&mut self.rng);
            let q_std: f64 = StandardNormal.sample(&mut self.rng);

            // Scale to desired standard deviation and convert to type T
            let i = T::from(i_std).unwrap() * std_dev_t;
            let q = T::from(q_std).unwrap() * std_dev_t;

            samples.push(Complex::new(i, q));
        }

        ComplexVec::from_vec(sample

    fn test_awgn_from_seed() {
        let mut awgn = AWGN::new_from_seed(1e6, 1000, 1.0, 42);
        let block1 = awgn.generate_block::<f64>();

        assert_eq!(block1.len(), 1000);

        // Generate another block and verify it's different (stateful)
        let block2 = awgn.generate_block::<f64>();
        assert_eq!(block2.len(), 1000);

        // Blocks should be different due to RNG state progression
        assert_ne!(block1[0], block2[0]);
    }

    fn test_awgn_reproducibility() {
        // Same seed should produce identical sequences
        let mut awgn1 = AWGN::new_from_seed(1e6, 100, 1.0, 123);
        let mut awgn2 = AWGN::new_from_seed(1e6, 100, 1.0, 123);

        let block1 = awgn1.generate_block::<f64>();
        let block2 = awgn2.generate_block::<f64>();

        // Should be identical
        for i in 0..block1.len() {
            assert_eq!(block1[i], block2[i]);
        }
    }

    fn test_awgn_from_entropy() {
        let mut awgn = AWGN::new_from_entropy(1e6, 500, 2.0);
        let block = awgn.generate_block::<f32>();

        assert_eq!(block.len(), 500);

        // Basic sanity check: samples should not all be zero
        let mut has_nonzero = false;
        for i in 0..block.len() {
            if block[i].re != 0.0 || block[i].im != 0.0 {
                has_nonzero = true;
                break;
            }
        }
        assert!(has_nonzero);
    }

    fn test_awgn_f32_and_f64() {
        let mut awgn = AWGN::new_from_seed(1e6, 100, 1.0, 999);

        // Test with f32
        let block_f32 = awgn.generate_block::<f32>();
        assert_eq!(block_f32.len(), 100);

        // Reset RNG
        awgn = AWGN::new_from_seed(1e6, 100, 1.0, 999);

        // Test with f64
        let block_f64 = awgn.generate_block::<f64>();
        assert_eq!(block_f64.len(), 100);

        // Values should be approximately equal (within f32 precision)
        for i in 0..block_f32.len() {
            let diff_re = (block_f32[i].re as f64 - block_f64[i].re).abs();
            let diff_im = (block_f32[i].im as f64 - block_f64[i].im).abs();
            assert!(diff_re < 1e-6, "Real part differs too much at {}: {} vs {}", i, block_f32[i].re, block_f64[i].re);
            assert!(diff_im < 1e-6, "Imag part differs too much at {}: {} vs {}", i, block_f32[i].im, block_f64[i].im);
        }
    }

    fn test_awgn_statistical_properties() {
        // Generate many samples to check statistical properties
        let mut awgn = AWGN::new_from_seed(1e6, 10000, 1.0, 777);
        let block = awgn.generate_block::<f64>();

        // Calculate mean (should be close to 0)
        let mut mean_i = 0.0;
        let mut mean_q = 0.0;
        for i in 0..block.len() {
            mean_i += block[i].re;
            mean_q += block[i].im;
        }
        mean_i /= block.len() as f64;
        mean_q /= block.len() as f64;

        // Mean should be close to 0 (within 3*sigma/sqrt(N) with high probability)
        let expected_mean_error = 3.0 * (0.5_f64.sqrt()) / (block.len() as f64).sqrt();
        assert!(mean_i.abs() < expected_mean_error, "Mean I = {}, expected < {}", mean_i, expected_mean_error);
        assert!(mean_q.abs() < expected_mean_error, "Mean Q = {}, expected < {}", mean_q, expected_mean_error);

        // Calculate variance (should be close to noise_power)
        let mu

    pub fn new(
        modulation: ModType,
        bandwidth: f64,
        center_freq: f64,
        snr_db: f64,
        rolloff: f64,
        sample_rate_hz: f64,
        seed: Option<u64>,
    ) -> Self {
        assert!(
            bandwidth > 0.0 && bandwidth <= 1.0,
            "bandwidth must be in range (0.0, 1.0]"
        );
        assert!(
            center_freq >= -0.5 && center_freq <= 0.5,
            "center_freq must be in range [-0.5, 0.5]"
        );
        assert!(
            rolloff >= 0.0 && rolloff <= 1.0,
            "rolloff must be in range [0.0, 1.0]"
        );
        assert!(sample_rate_hz > 0.0, "sample_rate_hz must be positive");

        Carrier {
            modulation,
            bandwidth,
            center_freq,
            snr_db,
            rolloff,
            sample_rate_hz,
            seed,
        }
    }