Function bodies 301 total

    pub fn generate_clean<T: Float>(&self, num_samples: usize) -> ComplexVec<T> {
        // Calculate symbol rate from bandwidth
        // bandwidth = (symbol_rate * (1 + rolloff)) / sample_rate
        // symbol_rate = (bandwidth * sample_rate) / (1 + rolloff)
        let symbol_rate = (self.bandwidth * self.sample_rate_hz) / (1.0 + self.rolloff);

        // Generate baseband signal using appropriate carrier type
        let mut iq_samples = match self.modulation {
            ModType::_BPSK | ModType::_QPSK | ModType::_8PSK | ModType::_16APSK
            | ModType::_16QAM | ModType::_32QAM | ModType::_64QAM => {
                self.generate_psk_based(symbol_rate, num_samples)
            }
            // FSK handling would go here if we support it
            ModType::_CW => {
                self.generate_cw(num_samples)
            }
        };

        // Frequency shift to center frequency (if not baseband)
        if self.center_freq.abs() > 1e-9 {
            let freq_offse

    pub fn generate<T: Float>(&self, num_samples: usize) -> ComplexVec<T> {
        let iq_samples = self.generate_clean::<T>(num_samples);

        // Measure signal power
        let oversample_rate = 1.0_f64 / self.bandwidth;
        let signal_power = iq_samples.measure_power(Some(oversample_rate));

        // Add AWGN to achieve target SNR
        // SNR_dB = 10 * log10(signal_power / noise_power)
        // noise_power = signal_power / 10^(SNR_dB/10)
        let snr_linear = 10.0_f64.powf(self.snr_db / 10.0);
        let noise_power = signal_power / snr_linear;

        let mut awgn = match self.seed {
            Some(s) => AWGN::new_from_seed(self.sample_rate_hz, num_samples, noise_power, s),
            None => AWGN::new_from_entropy(self.sample_rate_hz, num_samples, noise_power),
        };

        let noise = awgn.generate_block::<T>();
        let result = iq_samples + noise;

        result
    }

    fn generate_psk_based<T: Float>(&self, symbol_rate: f64, num_samples: usize) -> ComplexVec<T> {
        // Estimate block size (will be refined by PSK carrier)
        let block_size = (num_samples as f64 * symbol_rate / self.sample_rate_hz).ceil() as usize;

        let mut psk = PskCarrier::<T>::new(
            T::from(self.sample_rate_hz).unwrap(),
            T::from(symbol_rate).unwrap(),
            self.modulation,
            T::from(self.rolloff).unwrap(),
            block_size,
            91, // Default filter taps
            self.seed,
        );

        // Generate blocks until we have enough samples
        let mut result = ComplexVec::new();
        while result.len() < num_samples {
            let block = psk.generate_block();
            result.extend(block.iter().cloned());
        }

        // Trim to exact sample count
        let vec: Vec<_> = result
            .iter()
            .take(num_samples)
            .cloned()
            .collect();
        C

    fn generate_cw<T: Float>(&self, num_samples: usize) -> ComplexVec<T> {
        // CW is just a constant magnitude signal at baseband (0 Hz)
        // Generate as if it's a CW at frequency 0
        let mut cw = CW::new(0.0, self.sample_rate_hz, num_samples);
        let samples = cw.generate_block::<T>();
        ComplexVec::from_vec(samples)
    }

    fn test_carrier_qpsk_creation() {
        let carrier = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.1,
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        assert_eq!(carrier.modulation, ModType::_QPSK);
        assert_eq!(carrier.snr_db, 10.0);
    }

    fn test_carrier_generate_qpsk() {
        let carrier = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.0,
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        let iq = carrier.generate::<f64>(1000);
        assert_eq!(iq.len(), 1000);
    }

    fn test_carrier_generate_cw() {
        let carrier = Carrier::new(
            ModType::_CW,
            0.1,
            0.0,
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        let iq = carrier.generate::<f64>(1000);
        assert_eq!(iq.len(), 1000);
    }

    fn test_carrier_with_frequency_shift() {
        let carrier = Carrier::new(
            ModType::_CW,
            0.1,
            0.25, // Shift to 250 kHz
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        let iq = carrier.generate::<f64>(100);
        assert_eq!(iq.len(), 100);
    }

    fn test_invalid_bandwidth() {
        let result = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            Carrier::new(
                ModType::_QPSK,
                1.5, // Invalid: > 1.0
                0.0,
                10.0,
                0.35,
                1e6,
                None,
            );
        }));
        assert!(result.is_err());
    }

    fn test_invalid_center_freq() {
        let result = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
            Carrier::new(
                ModType::_QPSK,
                0.1,
                0.6, // Invalid: > 0.5
                10.0,
                0.35,
                1e6,
                None,
            );
        }));
        assert!(result.is_err());
    }

    fn test_carrier_combination_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;

        // Check if plotting is enabled
        if !should_plot() {
            println!("Skipping carrier combination spectrum plot (set PLOT=true to enable)");
            return;
        }

        let sample_rate = 1e6;

        // Create two carriers with contrasting characteristics
        // Carrier 1: Strong, narrow bandwidth signal
        let carrier1 = Carrier::new(
            ModType::_QPSK,
            0.05,     // 50 kHz bandwidth (narrow)
            0.15,     // 150 kHz center freq
            15.0,     // 25 dB SNR (high power)
            0.35,     // RRC rolloff
            sample_rate,
            Some(42),
        );

        // Carrier 2: Weak, wide bandwidth signal (interference/noise)
        let carrier2 = Carri

    fn test_carrier_with_awgn() {
        use num_complex::Complex;
        use crate::vector_ops;

        let carrier = Carrier::new(ModType::_QPSK, 0.25, 0.0, 10.0, 0.15, 1e6, Some(0));
        let _oversample_rate = 1.0f64 / carrier.bandwidth;
        let signal_iq = carrier.generate(1_000_000);

        let slice: &[Complex<f32>] = &signal_iq;
        let (freqs, spectrum) = welch(slice, 1e6, 2048, Some(512), Some(2048), crate::spectrum::window::WindowType::Rectangular, Some(AveragingMethod::Median));

        if crate::test_utils::should_plot() {
            plot_spectrum(&freqs, &vector_ops::to_db(&spectrum), "SNR Comparison");
        }
    }

    pub fn new(carriers: Vec<Carrier>) -> Self {
        Channel {
            carriers,
            noise_floor_db: None,
            seed: None,
            impairments: Vec::new(),
        }
    }

    fn noise_floor_to_power(&self) -> f64 {
        match self.noise_floor_db {
            Some(db) => 10.0_f64.powf(db / 10.0),
            None => panic!("Noise floor must be set before generating channel. Use set_noise_floor_db() or set_noise_floor_linear()"),
        }
    }

    pub fn generate<T: Float>(&self, num_samples: usize) -> ComplexVec<T> {
        if self.carriers.is_empty() {
            panic!("Channel must have at least one carrier");
        }

        // Get noise power from the noise floor (required for SNR scaling)
        let noise_power = self.noise_floor_to_power();

        // Generate, scale, and combine carriers to achieve target SNRs
        let mut combined = ComplexVec::new();
        for carrier in &self.carriers {
            // Generate clean signal
            let clean_signal = carrier.generate_clean::<T>(num_samples);

            // Calculate oversample rate for power measurement (accounts for in-band power)
            let oversample_rate = 1.0 / carrier.bandwidth;

            // Scale to achieve target SNR
            // SNR = P_signal / P_noise  =>  P_signal = SNR_linear * P_noise
            let snr_linear = 10.0_f64.powf(carrier.snr_db / 10.0);
            let target_power = snr_linear * noise_power;

            // S

    pub fn generate_clean<T: Float>(&self, num_samples: usize) -> ComplexVec<T> {
        if self.carriers.is_empty() {
            panic!("Channel must have at least one carrier");
        }

        // Get noise power from the noise floor (required for SNR scaling)
        let noise_power = self.noise_floor_to_power();

        // Generate, scale, and combine clean signals from all carriers
        let mut combined = ComplexVec::new();
        for carrier in &self.carriers {
            // Generate clean signal
            let clean_signal = carrier.generate_clean::<T>(num_samples);

            // Calculate oversample rate for power measurement
            let oversample_rate = 1.0 / carrier.bandwidth;

            // Scale to achieve target SNR
            let snr_linear = 10.0_f64.powf(carrier.snr_db / 10.0);
            let target_power = snr_linear * noise_power;

            // Scale the carrier to target power
            let scaled_signal = clean_signal.scale_to_power(target_

    fn test_channel_creation() {
        let carrier1 = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.1,
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        let channel = Channel::new(vec![carrier1]);
        assert_eq!(channel.num_carriers(), 1);
    }

    fn test_channel_single_carrier() {
        let carrier = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.1,
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        let mut channel = Channel::new(vec![carrier]);
        channel.set_noise_floor_db(-100.0);

        let result = channel.generate::<f64>(1000);
        assert_eq!(result.len(), 1000);
    }

    fn test_channel_two_carriers() {
        let carrier1 = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.1,
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        let carrier2 = Carrier::new(
            ModType::_QPSK,
            0.1,
            -0.1,
            15.0,
            0.35,
            1e6,
            Some(43),
        );

        let mut channel = Channel::new(vec![carrier1, carrier2]);
        channel.set_noise_floor_db(-100.0);

        let result = channel.generate::<f64>(1000);
        assert_eq!(result.len(), 1000);
    }

    fn test_channel_generate_clean() {
        let carrier = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.1,
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        let mut channel = Channel::new(vec![carrier]);
        channel.set_noise_floor_db(-100.0);

        let clean = channel.generate_clean::<f64>(1000);
        assert_eq!(clean.len(), 1000);
    }

    fn test_channel_noise_floor_linear() {
        let carrier = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.1,
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        let mut channel = Channel::new(vec![carrier]);

        // Set noise floor as linear value
        channel.set_noise_floor_linear(0.01);

        let result = channel.generate::<f64>(1000);
        assert_eq!(result.len(), 1000);
    }

    fn test_channel_with_seed() {
        let carrier1 = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.1,
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        let carrier2 = Carrier::new(
            ModType::_QPSK,
            0.1,
            -0.1,
            10.0,
            0.35,
            1e6,
            Some(43),
        );

        // Generate with same seed - should get same result
        let mut channel1 = Channel::new(vec![carrier1.clone(), carrier2.clone()]);
        channel1.set_noise_floor_db(-100.0);
        channel1.set_seed(999);

        let mut channel2 = Channel::new(vec![carrier1, carrier2]);
        channel2.set_noise_floor_db(-100.0);
        channel2.set_seed(999);

        let result1 = channel1.generate::<f64>(1000);
        let result2 = channel2.generate::<f64>(1000);

        // Results should be identical with same seed
        for (s1, s2) in result1.iter().zip(result2.iter()) {
    

    fn test_channel_no_noise_floor() {
        let carrier = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.1,
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        let channel = Channel::new(vec![carrier]);

        // Should panic because noise floor was never set
        let _ = channel.generate::<f64>(1000);
    }

    fn test_channel_spectrum_with_snr_verification() {
        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;

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

        // Create two carriers with different target SNRs
        let carrier1 = Carrier::new(
            ModType::_QPSK,
            0.05,     // 50 kHz bandwidth (narrow)
            0.15,     // 150 kHz center freq
            10.0,     // Target 10 dB SNR
            0.35,     // RRC rolloff
            sample_rate,
            Some(42),
        );

        let carrier2 = Carrier::new(
            ModType::_16QAM,
            0.10,     // 100 kHz bandwidth (wider)
            -0.20,    // -200 kHz center freq
            5.0,      // Target 5 dB SNR
            0.35,
            sample_rate,
            Some(43),
        );

        // Def

    fn test_channel_with_single_impairment() {
        let carrier = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.1,
            10.0,
            0.35,
            1e6,
            Some(42),
        );
        let mut channel = Channel::new(vec![carrier]);
        channel.set_noise_floor_db(-100.0);
        channel.add_impairment(Impairment::DigitizerDroopAD9361);

        let result = channel.generate::<f64>(1000);
        assert_eq!(result.len(), 1000);
    }

    fn test_channel_with_multiple_impairments() {
        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;

        let sample_rate = 1e6;
        let num_samples = 1_000_000;

        let carrier = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.1,
            10.0,
            0.35,
            sample_rate,
            Some(1),
        );

        let carrier2 = Carrier::new(
            ModType::_QPSK,
            0.05,
            -0.25,
            3.0,
            0.20,
            sample_rate,
            Some(2),
        );

        let carrier3 = Carrier::new(
            ModType::_QPSK,
            0.02,
            0.2,
            5.0,
            0.15,
            sample_rate,
            Some(3),
        );

        let mut channel = Channel::new(vec![carrier, carrier2, carrier3]);
        channel.se

    fn test_channel_with_custom_digitizer_droop() {
        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;

        let sample_rate = 1e6;
        let num_samples = 25000;

        let carrier = Carrier::new(
            ModType::_QPSK,
            0.1,
            0.0,
            10.0,
            0.35,
            sample_rate,
            Some(42),
        );

        let mut channel = Channel::new(vec![carrier]);
        channel.set_noise_floor_db(-50.0);
        channel.set_seed(999);

        // Add custom digitizer droop (4th order, higher cutoff)
        channel.add_impairment(Impairment::DigitizerDroop {
            order: 9,
            cutoff: 0.45,
        });

        let result = channel.generate::<f64>(num_samples);
        assert_eq!(result.len(), num_samples);

        println!("\n=== Channel with Custom Digitizer Droop

    pub fn new(freq_hz: f64, sample_rate_hz: f64, block_size: usize) -> Self {
        CW {
            freq_hz: freq_hz,
            sample_rate_hz: sample_rate_hz,
            block_size: block_size,
            sample_num: 0,
        }
    }

    pub fn generate_block<T: Float>(&mut self) -> Vec<Complex<T>> {
        let mut cw_block = Vec::with_capacity(self.block_size);
        for _ in 0..self.block_size {
            let t: T = T::from(self.sample_num).unwrap() / T::from(self.sample_rate_hz).unwrap();
            let i = (T::from(2.0).unwrap() * T::from(PI).unwrap() * T::from(self.freq_hz).unwrap() * t).cos();
            let q = (T::from(2.0).unwrap() * T::from(PI).unwrap() * T::from(self.freq_hz).unwrap() * t).sin();
            cw_block.push(Complex::<T>::new(i, q));
            self.sample_num += 1;
        }
        cw_block
    }

    fn test_cw_block_gen() {
        let block_size: usize = 1024;
        let num_blocks = 100;
        let freq_hz = -100e3;
        let sample_rate_hz = 500e3;
        let mut generator = CW::new(freq_hz, sample_rate_hz, block_size);

        let mut samples: Vec<_> = (0..num_blocks)
            .flat_map(|_| generator.generate_block::<f64>())
            .collect();

        assert_eq!(samples.len(), block_size * num_blocks);

        if crate::test_utils::should_plot() {
            fft::fft::<f64>(&mut samples);
            let mut cw_fft = ComplexVec::from_vec(samples);
            let mut cw_fft_abs: Vec<f64> = cw_fft.abs();

            fft::fftshift::<f64>(&mut cw_fft_abs);
            let freqs: Vec<f64> = fft::fftfreqs::<f64>((-sample_rate_hz/2_f64) as f64, (sample_rate_hz/2_f64) as f64, cw_fft_abs.len());

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

    pub fn new(
        sample_rate_hz: f64,
        symbol_rate_hz: f64,
        carrier_freq_hz: f64,
        freq_deviation_hz: f64,
        block_size: usize,
        seed: Option<u64>,
    ) -> Self {
        let samples_per_symbol = (sample_rate_hz / symbol_rate_hz).round() as usize;

        let bit_gen = match seed {
            Some(s) => BitGenerator::new_from_seed(s),
            None => BitGenerator::new_from_entropy(),
        };

        FskCarrier {
            sample_rate_hz: T::from(sample_rate_hz).unwrap(),
            symbol_rate_hz: T::from(symbol_rate_hz).unwrap(),
            carrier_freq_hz: T::from(carrier_freq_hz).unwrap(),
            freq_deviation_hz: T::from(freq_deviation_hz).unwrap(),
            block_size,
            accumulated_phase: T::zero(),
            bit_gen,
            samples_per_symbol,
        }
    }

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

        let two_pi = T::from(2.0 * PI).unwrap();
        let half = T::from(0.5).unwrap();

        for _ in 0..num_symbols {
            // Get next bit
            let bit = self.bit_gen.next_bit();

            // Calculate instantaneous frequency for this symbol
            // bit 0 -> fc - Df/2
            // bit 1 -> fc + Df/2
            let freq = if bit {
                self.carrier_freq_hz + self.freq_deviation_hz * half
            } else {
                self.carrier_freq_hz - self.freq_deviation_hz * half
            };

            // Generate samples for this symbol
            for _ in 0..self.samples_per_symbol {
                // Calculate phase increment for this sample
                let phase_increment = two_pi * freq / self.sample_rate_hz;

                // Accumulate phas

    fn test_fsk_basic() {
        let sample_rate = 1e6;
        let symbol_rate = 1e5;
        let carrier_freq = 2.5e5;
        let deviation = 5e4;
        let block_size = 1000;

        let mut fsk: FskCarrier<f64> = FskCarrier::new(
            sample_rate,
            symbol_rate,
            carrier_freq,
            deviation,
            block_size,
            Some(42),
        );

        let block = fsk.generate_block();

        // Verify block size
        assert_eq!(block.len(), block_size);

        // Verify samples are unit magnitude (approximately)
        for i in 0..block.len() {
            let mag = block[i].norm();
            assert!((mag - 1.0).abs() < 1e-10, "Sample {} has magnitude {}, expected 1.0", i, mag);
        }
    }

    fn test_fsk_phase_continuity() {
        let sample_rate = 1e6;
        let symbol_rate = 1e5;
        let carrier_freq = 2.5e5;
        let deviation = 5e4;
        let block_size = 100;

        let mut fsk: FskCarrier<f64> = FskCarrier::new(
            sample_rate,
            symbol_rate,
            carrier_freq,
            deviation,
            block_size,
            Some(123),
        );

        // Generate multiple blocks and verify phase continuity
        let block1 = fsk.generate_block();
        let block2 = fsk.generate_block();

        // Blocks should not be identical (different random bits)
        let blocks_identical = block1.iter().zip(block2.iter()).all(|(a, b)| a == b);
        assert!(!blocks_identical, "Blocks should not be identical due to different random bits");

        // All samples should have unit magnitude
        for i in 0..block1.len() {
            assert!((block1[i].norm() - 1.0).abs() < 1e-10);
        }
        for i in 0..block2.len() {

    fn test_fsk_modulation_index() {
        let fsk: FskCarrier<f64> = FskCarrier::new(1e6, 1e5, 2.5e5, 5e4, 1000, Some(42));

        let h = fsk.modulation_index();

        // For these parameters: h = 50000 / 100000 = 0.5 (MSK)
        assert!((h - 0.5).abs() < 1e-10);
    }

    fn test_fsk_with_awgn_spectrum() {
        use crate::test_utils::should_plot;
        
        if !should_plot() {
            println!("Skipping FSK with AWGN spectrum plot (set PLOT=true to enable)");
            return;
        }

        println!("\n=== FSK with AWGN Spectrum Test ===");

        // FSK parameters
        let sample_rate = 1e6_f64;
        let symbol_rate = 1e5_f64;
        let carrier_freq = 2.5e5_f64;
        let deviation = 5e4_f64;  // Two frequencies: 200 kHz and 300 kHz
        let block_size = 4096;

        println!("Sample rate: {} Hz", sample_rate);
        println!("Symbol rate: {} Hz", symbol_rate);
        println!("Carrier freq: {} Hz", carrier_freq);
        println!("Deviation: {} Hz", deviation);
        println!("Frequencies: {} Hz and {} Hz",
                 carrier_freq - deviation/2.0,
                 carrier_freq + deviation/2.0);

        // Generate FSK signal
        let mut fsk: FskCarrier<f64> = FskCarrier::new(
            sample_

    pub fn apply(&self, signal: &mut [Complex<f64>]) {
        match self {
            Impairment::DigitizerDroop { order, cutoff } => {
                apply_butterworth_filter(signal, *order, *cutoff);
            }
            Impairment::DigitizerDroopAD9361 => {
                apply_digitizer_droop_ad9361(signal);
            }
            Impairment::DigitizerDroopTraditional => {
                apply_digitizer_droop_traditional(signal);
            }
            Impairment::FrequencyVariation {
                amplitude_db,
                cycles,
                phase_offset,
            } => {
                // Generate frequency-domain variation pattern
                let variation = frequency_dependent_amplitude_variation(
                    signal.len(),
                    *amplitude_db,
                    *cycles,
                    *phase_offset,
                );
                let variation_lin = to_linear(&variation);

                // Apply in frequency d

pub fn frequency_dependent_amplitude_variation<T: Float>(num_samples: usize, amplitude_db: T, cycles: T, phase_offset: T,
) -> Vec<T> {
    let pi = T::from(PI).unwrap();
    (0..num_samples)
        .into_iter()
        .map(|idx| {
            amplitude_db
                * (T::from(2.0).unwrap() * pi * T::from(idx).unwrap()
                    / T::from(num_samples).unwrap()
                    * cycles
                    + phase_offset)
                    .sin()
        })
            .collect()
    }

    fn test_freq_ampl_variation() {
        use crate::test_utils::should_plot;

        if !should_plot() {
            println!("Skipping frequency-dependent amplitude variation plot (set PLOT=true to enable)");
            return;
        }

        let sample_rate = 1e6_f64;
        let num_samples = (2.0_f64).powf(20.0) as usize;  // Use power of 2 for better FFT performance
        // Generate three different frequency variation patterns
        let freq_var_1 = frequency_dependent_amplitude_variation(num_samples, 1.0, 2.0, 0.0);
        let freq_var_2 = frequency_dependent_amplitude_variation(num_samples, 0.6, 5.0, 0.0);
        let freq_var_3 = frequency_dependent_amplitude_variation(num_samples, 0.4, 1.0, 1.5);
        // Combine all variations
        let total_variation_db = add(&add(&freq_var_1, &freq_var_2), &freq_var_3);
        let total_variation_lin = to_linear(&total_variation_db);
        // Generate white noise
        let mut awgn = AWGN::new_from_seed(sample_rate,

    fn test_random_freq_ampl_variation() {
        use std::env;
        use std::f64::consts::PI;

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

        let seed = 2; // or any u64
        let mut rng = StdRng::seed_from_u64(seed);
        let sample_rate = 1e6_f64;
        let num_samples = (2.0_f64).powf(20.0) as usize;  // Use power of 2 for better FFT performance

        let mut awgn = AWGN::new_from_seed(sample_rate, num_samples, 1.0, 0);
        let mut noise: ComplexVec<f32> = awgn.generate_block();
        fft::fft(&mut noise);

        // Apply frequency-dependent amplitude variation in frequency domain
        let num_variations = rand::thread_rng().gen_range(1..=5);
        let mut freq_var_total: Vec<f32> = vec![0.0; num_samples];
        for _ in 0..num_variations 

    fn test_digitizer_droop() {
        use std::env;

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

        println!("\n=== Digitizer Channel Droop Impairment ===\n");

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

        // Generate white noise
        println!("Generating AWGN...");
        let mut awgn = AWGN::new_from_seed(sample_rate, n_samples, 1.0, 0);
        let noise: Vec<Complex<f64>> = awgn.generate_block::<f64>().to_vec();

        // 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_digitizer_droop_comparison() {
        use std::env;

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

        println!("\n=== Digitizer Droop: Order Comparison ===\n");

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

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

        // Test different filter orders with same cutoff
        let cutoff = 0.45;
        let orders = vec![2, 3, 4, 6, 8];

        for order in orders {
            let mut filtered = noise.clone();
            apply_digitizer_droop(&mut filtered, order, cutoff);

            let (freqs, psd) = welch(
                &filtered,
              

    fn test_digitizer_droop_on_carrier() {
        use std::env;
        use crate::generate::carrier::Carrier;
        use crate::ModType;

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

        println!("\n=== Digitizer Droop Applied to QPSK Carrier ===\n");

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

        // Create a QPSK carrier
        let carrier = Carrier::new(
            ModType::_QPSK,
            0.1,        // 10% bandwidth
            0.0,        // Centered
            15.0,       // 15 dB SNR
            0.35,       // RRC rolloff
            sample_rate,
            Some(42),
        );

        // Generate clean carrier signal
        println!("Generating QPSK carrier (clean)...");
        let clean_signal: Vec<Complex<f64>> = car

    pub fn new(sample_rate_hz: T, symbol_rate_hz: T, mod_type: ModType, rolloff_factor: T, block_size: usize, filter_taps: usize, seed: Option<u64>) -> Self {
        let bit_gen = match seed {
            Some(s) => BitGenerator::new_from_seed(s),
            None => BitGenerator::new_from_entropy(),
        };

        // Calculate integer upsampling factor
        // For 1MHz sample rate / 800kHz symbol rate = 1.25, we round to nearest integer
        let sps_exact = sample_rate_hz.to_f64().unwrap() / symbol_rate_hz.to_f64().unwrap();
        let samples_per_symbol = sps_exact.round() as usize;
        if samples_per_symbol == 0 {
            panic!("samples_per_symbol must be at least 1");
        }

        // Build RRC filter once during initialization
        let rrc = RRCFilter::new(
            filter_taps,
            sample_rate_hz.to_f64().unwrap(),
            symbol_rate_hz.to_f64().unwrap(),
            rolloff_factor.to_f64().unwrap(),
        );
        let rrc_filter 

    pub fn generate_block(&mut self) -> ComplexVec<T> {
        // Calculate how many symbols we need to generate to get block_size samples
        // For sps=1: 1 symbol -> 1 sample (after upsampling and filtering)
        // We need to account for the overlap and filter delay
        let num_symbols = if self.samples_per_symbol == 1 {
            self.block_size
        } else {
            // For upsampling: symbols * sps ≈ samples
            (self.block_size + self.samples_per_symbol - 1) / self.samples_per_symbol
        };

        // Generate symbols for this block using the pre-built modulation
        let new_symbols = self.generate_symbols(num_symbols);

        // Combine overlap from previous block with new symbols
        let mut all_symbols = self.overlap_symbols.clone();
        all_symbols.extend(new_symbols.iter().cloned());

        // Upsample symbols by inserting zeros between them
        let upsampled = if self.samples_per_symbol > 1 {
            self.upsample_

    fn upsample_symbols(&self, symbols: &ComplexVec<T>) -> ComplexVec<T> {
        let mut upsampled = Vec::new();
        let zero = num_complex::Complex::new(T::zero(), T::zero());

        for i in 0..symbols.len() {
            upsampled.push(symbols[i]);
            // Insert (sps - 1) zeros after each symbol (except the last one)
            if i < symbols.len() - 1 {
                for _ in 0..(self.samples_per_symbol - 1) {
                    upsampled.push(zero);
                }
            }
        }

        ComplexVec::from_vec(upsampled)
    }

    fn generate_symbols(&mut self, count: usize) -> ComplexVec<T> {
        let mut symbols = Vec::new();

        match &self.modulation {
            modulation::Modulation::BPSK(bpsk) => {
                for _ in 0..count {
                    let bit = self.bit_gen.next_bit();
                    let bits = if bit { 1u8 } else { 0u8 };
                    if let Some(symbol) = bpsk.modulate(bits) {
                        symbols.push(symbol);
                    }
                }
            },
            modulation::Modulation::QPSK(qpsk) => {
                for _ in 0..count {
                    let bits = self.bit_gen.next_2_bits();
                    if let Some(symbol) = qpsk.modulate(bits) {
                        symbols.push(symbol);
                    }
                }
            },
            modulation::Modulation::PSK8(psk8) => {
                for _ in 0..count {
                    let bits = self.bit_gen.next_3_bits();
                    if let Some(

    fn extract_last_symbols(&self, symbols: &ComplexVec<T>, count: usize) -> ComplexVec<T> {
        let start_idx = if symbols.len() > count {
            symbols.len() - count
        } else {
            0
        };

        let mut last_symbols = Vec::new();
        for i in start_idx..symbols.len() {
            last_symbols.push(symbols[i]);
        }
        ComplexVec::from_vec(last_symbols)
    }

    fn test_qpsk_carrier_block_gen() {
        let block_size: usize = 1e6 as usize;
        let num_blocks = 1;
        let sample_rate_hz = 1e6_f64;
        let symbol_rate_hz = 800e3_f64;
        let filter_taps = 51;
        let rolloff_factor = 0.05_f64;

        let mut carrier = PskCarrier::new(
            sample_rate_hz,
            symbol_rate_hz,
            ModType::_QPSK,
            rolloff_factor,
            block_size,
            filter_taps,
            Some(42), // Seed for reproducibility
        );

        println!("\n=== QPSK Carrier Configuration ===");
        println!("Block size: {}", block_size);
        println!("Sample rate: {} Hz", sample_rate_hz);
        println!("Symbol rate: {} Hz", symbol_rate_hz);
        println!("Samples per symbol: {}", carrier.samples_per_symbol);
        println!("Filter taps: {}", filter_taps);
        println!("Rolloff factor: {}", rolloff_factor);

        // Generate blocks and merge into single ComplexVec
        let mut 

    fn test_qpsk_symbols_constellation() {
        use std::env;
        use plotly::{Plot, Scatter};
        use plotly::common::Mode;

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

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

        let mut carrier = PskCarrier::new(
            sample_rate_hz,
            symbol_rate_hz,
            ModType::_QPSK,
            rolloff_factor,
            1000, // Small block for visualization
            filter_taps,
            Some(42),
        );

        // Generate just the raw symbols (before pulse shaping) for visualization
        let num_symbols = 1000;
        let symbols = carrier.generate_symbols(num_symbols);

        // Extract I and Q v