diff --git a/python/tests/unit/aout/test_analyze_error_by_phase.py b/python/tests/unit/aout/test_analyze_error_by_phase.py
index 548fc6c..485510d 100644
--- a/python/tests/unit/aout/test_analyze_error_by_phase.py
+++ b/python/tests/unit/aout/test_analyze_error_by_phase.py
@@ -23,6 +23,7 @@ def test_analyze_error_by_phase_pure_cases():
     A = 0.49
     DC = 0.5
     phase_clean = 2 * np.pi * Fin * t
+    rng = np.random.default_rng(2026062210)
 
     print(f"\n[Config] Fs={Fs/1e6:.0f} MHz, Fin={Fin/1e6:.2f} MHz, N={N}, A={A}")
 
@@ -38,9 +39,9 @@ def test_analyze_error_by_phase_pure_cases():
 
     # Generate signals
     for case in test_cases:
-        am_noise = np.random.randn(N) * case['am'] if case['am'] > 0 else 0
-        pm_noise = np.random.randn(N) * case['pm'] / A if case['pm'] > 0 else 0
-        th_noise = np.random.randn(N) * case['thermal'] if case['thermal'] > 0 else 0
+        am_noise = rng.standard_normal(N) * case['am'] if case['am'] > 0 else 0
+        pm_noise = rng.standard_normal(N) * case['pm'] / A if case['pm'] > 0 else 0
+        th_noise = rng.standard_normal(N) * case['thermal'] if case['thermal'] > 0 else 0
         case['signal'] = (A + am_noise) * np.sin(phase_clean + pm_noise) + DC + th_noise
 
     # Test both baseline modes
@@ -93,6 +94,7 @@ def test_analyze_error_by_phase_mixed_cases():
     A = 0.49
     DC = 0.5
     phase_clean = 2 * np.pi * Fin * t
+    rng = np.random.default_rng(2026062211)
 
     # Define mixed test cases
     test_cases = [
@@ -106,9 +108,9 @@ def test_analyze_error_by_phase_mixed_cases():
 
     # Generate signals
     for case in test_cases:
-        am_noise = np.random.randn(N) * case['am'] if case['am'] > 0 else 0
-        pm_noise = np.random.randn(N) * case['pm'] / A if case['pm'] > 0 else 0
-        th_noise = np.random.randn(N) * case['thermal'] if case['thermal'] > 0 else 0
+        am_noise = rng.standard_normal(N) * case['am'] if case['am'] > 0 else 0
+        pm_noise = rng.standard_normal(N) * case['pm'] / A if case['pm'] > 0 else 0
+        th_noise = rng.standard_normal(N) * case['thermal'] if case['thermal'] > 0 else 0
         case['signal'] = (A + am_noise) * np.sin(phase_clean + pm_noise) + DC + th_noise
 
     fig, axes = plt.subplots(1, 3, figsize=(18, 8))
diff --git a/python/tests/unit/aout/test_analyze_error_by_value.py b/python/tests/unit/aout/test_analyze_error_by_value.py
index d025b53..f393a9e 100644
--- a/python/tests/unit/aout/test_analyze_error_by_value.py
+++ b/python/tests/unit/aout/test_analyze_error_by_value.py
@@ -22,15 +22,16 @@ def test_analyze_error_by_value_thermal_vs_nonlinearity():
     A = 0.49
     DC = 0.5
     base_noise = 50e-6
+    rng = np.random.default_rng(2026062212)
 
     print(f"\n[Config] Fs={Fs/1e6:.0f} MHz, Fin={Fin/1e6:.2f} MHz, N={N}")
 
     # Case 1: Ideal ADC with Thermal Noise
-    sig_noise = A * np.sin(2 * np.pi * Fin * t) + DC + np.random.randn(N) * base_noise
+    sig_noise = A * np.sin(2 * np.pi * Fin * t) + DC + rng.standard_normal(N) * base_noise
 
     # Case 2: ADC with 3rd Order Nonlinearity
     k3 = 0.01
-    sig_nonlin = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + np.random.randn(N) * base_noise
+    sig_nonlin = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + rng.standard_normal(N) * base_noise
 
     # Create figure with 3 subplots
     fig, axes = plt.subplots(1, 3, figsize=(18, 8))
@@ -65,8 +66,9 @@ def test_analyze_error_by_value_bin_count(n_bins):
     DC = 0.5
     base_noise = 50e-6
     k3 = 0.01
+    rng = np.random.default_rng(2026062213)
 
-    sig_nonlin = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + np.random.randn(N) * base_noise
+    sig_nonlin = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + rng.standard_normal(N) * base_noise
 
     # Create a simple plot
     fig, ax = plt.subplots(1, 1, figsize=(6, 5))
@@ -86,17 +88,18 @@ def test_analyze_error_by_value_different_nonlinearities():
     A = 0.49
     DC = 0.5
     base_noise = 50e-6
+    rng = np.random.default_rng(2026062214)
 
     # Case 1: Thermal only
-    sig_thermal = A * np.sin(2 * np.pi * Fin * t) + DC + np.random.randn(N) * base_noise
+    sig_thermal = A * np.sin(2 * np.pi * Fin * t) + DC + rng.standard_normal(N) * base_noise
 
     # Case 2: 2nd order nonlinearity
     k2 = 0.01
-    sig_k2 = A * np.sin(2 * np.pi * Fin * t) + DC + k2 * (A * np.sin(2 * np.pi * Fin * t))**2 + np.random.randn(N) * base_noise
+    sig_k2 = A * np.sin(2 * np.pi * Fin * t) + DC + k2 * (A * np.sin(2 * np.pi * Fin * t))**2 + rng.standard_normal(N) * base_noise
 
     # Case 3: 3rd order nonlinearity
     k3 = 0.01
-    sig_k3 = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + np.random.randn(N) * base_noise
+    sig_k3 = A * np.sin(2 * np.pi * Fin * t) + DC + k3 * (A * np.sin(2 * np.pi * Fin * t))**3 + rng.standard_normal(N) * base_noise
 
     # Create figure with 3 subplots
     fig, axes = plt.subplots(1, 3, figsize=(18, 8))
diff --git a/python/tests/unit/aout/test_decompose_harmonics.py b/python/tests/unit/aout/test_decompose_harmonics.py
index 946825b..7168eba 100644
--- a/python/tests/unit/aout/test_decompose_harmonics.py
+++ b/python/tests/unit/aout/test_decompose_harmonics.py
@@ -22,25 +22,26 @@ def test_decompose_harmonics_basic():
     A = 0.25
     DC = 0.5
     base_noise = 50e-6
+    rng = np.random.default_rng(2026062215)
 
     sig_ac = A * np.sin(2 * np.pi * Fin * t)
     sig_ideal = sig_ac + DC
     print(f"\n[Config] Fs={Fs/1e6:.0f} MHz, Fin={Fin/1e6:.2f} MHz, N={N}")
 
     # Case 1: Ideal ADC with Thermal Noise
-    sig_noise = sig_ideal + np.random.randn(N) * base_noise
+    sig_noise = sig_ideal + rng.standard_normal(N) * base_noise
 
     # Case 2: ADC with Nonlinearity
     k2 = 0.001
     k3 = 0.005
-    sig_nonlin = DC + sig_ac + k2 * sig_ac**2 + k3 * sig_ac**3 + np.random.randn(N) * base_noise
+    sig_nonlin = DC + sig_ac + k2 * sig_ac**2 + k3 * sig_ac**3 + rng.standard_normal(N) * base_noise
 
     # Case 3: ADC with Glitches
     glitch_prob = 0.01
     glitch_amplitude = 0.1
-    glitch_mask = np.random.rand(N) < glitch_prob
+    glitch_mask = rng.random(N) < glitch_prob
     glitch = glitch_mask * glitch_amplitude
-    sig_glitch = sig_ideal + glitch + np.random.randn(N) * base_noise
+    sig_glitch = sig_ideal + glitch + rng.standard_normal(N) * base_noise
 
     # Create figure
     fig, axes = plt.subplots(1, 3, figsize=(18, 8))
@@ -74,9 +75,10 @@ def test_decompose_harmonics_nonlinearity_levels(k3_value):
     A = 0.25
     DC = 0.5
     base_noise = 50e-6
+    rng = np.random.default_rng(2026062216)
 
     sig_ac = A * np.sin(2 * np.pi * Fin * t)
-    sig_nonlin = DC + sig_ac + k3_value * sig_ac**3 + np.random.randn(N) * base_noise
+    sig_nonlin = DC + sig_ac + k3_value * sig_ac**3 + rng.standard_normal(N) * base_noise
 
     # Just run the analysis without creating plot
     fig, ax = plt.subplots(1, 1, figsize=(6, 5))
diff --git a/python/tests/unit/aout/test_decompose_harmonics_polar.py b/python/tests/unit/aout/test_decompose_harmonics_polar.py
index 9fb5f3b..164e954 100644
--- a/python/tests/unit/aout/test_decompose_harmonics_polar.py
+++ b/python/tests/unit/aout/test_decompose_harmonics_polar.py
@@ -24,25 +24,26 @@ def test_decompose_harmonics_polar_basic():
     DC = 0.5
     base_noise = 50e-6
     adc_range = [0, 1]
+    rng = np.random.default_rng(2026062217)
 
     sig_ac = A * np.sin(2 * np.pi * Fin * t)
     sig_ideal = sig_ac + DC
     print(f"\n[Config] Fs={Fs/1e6:.0f} MHz, Fin={Fin/1e6:.6f} MHz, Bin={Fin_bin}, N={N}")
 
     # Case 1: Ideal ADC with Thermal Noise
-    sig_noise = sig_ideal + np.random.randn(N) * base_noise
+    sig_noise = sig_ideal + rng.standard_normal(N) * base_noise
 
     # Case 2: ADC with Nonlinearity
     k2 = 0.001
     k3 = 0.005
-    sig_nonlin = DC + sig_ac + k2 * sig_ac**2 + k3 * sig_ac**3 + np.random.randn(N) * base_noise
+    sig_nonlin = DC + sig_ac + k2 * sig_ac**2 + k3 * sig_ac**3 + rng.standard_normal(N) * base_noise
 
     # Case 3: ADC with Glitches
     glitch_prob = 0.01
     glitch_amplitude = 0.1
-    glitch_mask = np.random.rand(N) < glitch_prob
+    glitch_mask = rng.random(N) < glitch_prob
     glitch = glitch_mask * glitch_amplitude
-    sig_glitch = sig_ideal + glitch + np.random.randn(N) * base_noise
+    sig_glitch = sig_ideal + glitch + rng.standard_normal(N) * base_noise
 
     # Create 2x3 subplot grid
     fig = plt.figure(figsize=(18, 12))
diff --git a/python/tests/unit/aout/test_fit_sine_4param.py b/python/tests/unit/aout/test_fit_sine_4param.py
index 34d59e5..73b67da 100644
--- a/python/tests/unit/aout/test_fit_sine_4param.py
+++ b/python/tests/unit/aout/test_fit_sine_4param.py
@@ -21,11 +21,12 @@ def test_fit_sine_4param_basic():
     A = 0.499
     DC = 0.5
     noise_rms = 20e-3
+    rng = np.random.default_rng(2026062218)
 
     print(f"\n[Sinewave] [Fs={Fs/1e6:.1f} MHz] [Fin={Fin/1e6:.6f} MHz] [Amplitude={A:.3f} V] [DC={DC:.3f} V] [Noise RMS={noise_rms*1e3:.2f} mV]")
 
     sig_ideal = A * np.cos(2 * np.pi * Fin * t) + DC
-    sig_noisy = sig_ideal + np.random.randn(N) * noise_rms
+    sig_noisy = sig_ideal + rng.standard_normal(N) * noise_rms
 
     # Fit 4-parameter sine wave
     result = fit_sine_4param(sig_noisy)
@@ -110,9 +111,10 @@ def test_fit_sine_4param_noise_levels(noise_level):
     t = np.arange(N) / Fs
     A = 0.499
     DC = 0.5
+    rng = np.random.default_rng(2026062219)
 
     sig_ideal = A * np.cos(2 * np.pi * Fin * t) + DC
-    sig_noisy = sig_ideal + np.random.randn(N) * noise_level
+    sig_noisy = sig_ideal + rng.standard_normal(N) * noise_level
 
     result = fit_sine_4param(sig_noisy)
     sig_fit = result['fitted_signal']
diff --git a/python/tests/unit/aout/test_inl_from_sine_sweep_length.py b/python/tests/unit/aout/test_inl_from_sine_sweep_length.py
index f4a0c01..40b0773 100644
--- a/python/tests/unit/aout/test_inl_from_sine_sweep_length.py
+++ b/python/tests/unit/aout/test_inl_from_sine_sweep_length.py
@@ -25,6 +25,7 @@ def test_inl_from_sine_sweep_length():
     DC = 0.5
     base_noise = 50e-6
     hd2_dB, hd3_dB = -80, -66
+    rng = np.random.default_rng(2026062220)
 
     # Compute HD coefficients
     hd2_amp = 10**(hd2_dB/20)
@@ -43,7 +44,7 @@ def test_inl_from_sine_sweep_length():
         fin, J = find_coherent_frequency(fs, fin_target, N)
         t = np.arange(N) / fs
         sinewave = A * np.sin(2 * np.pi * fin * t)
-        signal_distorted = sinewave + k2 * sinewave**2 + k3 * sinewave**3 + DC + np.random.randn(N) * base_noise
+        signal_distorted = sinewave + k2 * sinewave**2 + k3 * sinewave**3 + DC + rng.standard_normal(N) * base_noise
 
         result = analyze_spectrum(signal_distorted, fs=fs, create_plot=False)
 
diff --git a/python/tests/unit/aout/test_jitter_calculation.py b/python/tests/unit/aout/test_jitter_calculation.py
index 9cc8ca7..4a99ae2 100644
--- a/python/tests/unit/aout/test_jitter_calculation.py
+++ b/python/tests/unit/aout/test_jitter_calculation.py
@@ -21,6 +21,7 @@ def test_jitter_recovery_at_1ghz():
     A = 0.49
     DC = 0.0
     base_noise = 50e-6
+    rng = np.random.default_rng(2026062221)
 
     # Test 3 frequencies
     fin_targets = [100e6, 1000e6, 2000e6]  # 100 MHz, 1 GHz, 2 GHz
@@ -52,9 +53,9 @@ def test_jitter_recovery_at_1ghz():
 
                 # Phase jitter model
                 phase_noise_rms = 2 * np.pi * Fin * jitter_rms
-                phase_jitter = np.random.randn(N) * phase_noise_rms
+                phase_jitter = rng.standard_normal(N) * phase_noise_rms
 
-                signal = A * np.sin(2*np.pi*Fin*t + phase_jitter) + DC + np.random.randn(N) * noise_level
+                signal = A * np.sin(2*np.pi*Fin*t + phase_jitter) + DC + rng.standard_normal(N) * noise_level
 
                 # Measure jitter using analyze_error_by_phase
                 results = analyze_error_by_phase(signal, norm_freq=Fin/Fs, n_bins=100,
@@ -164,6 +165,7 @@ def test_jitter_recovery_frequency_sweep(fin_target, expected_correlation):
     Fs = 7e9
     A = 0.49
     DC = 0.0
+    rng = np.random.default_rng(2026062222)
 
     # Find coherent frequency
     Fin, Fin_bin = find_coherent_frequency(fs=Fs, fin_target=fin_target, n_fft=N)
@@ -179,7 +181,7 @@ def test_jitter_recovery_frequency_sweep(fin_target, expected_correlation):
         # Generate signal
         t = np.arange(N) / Fs
         phase_noise_rms = 2 * np.pi * Fin * jitter_rms
-        phase_jitter = np.random.randn(N) * phase_noise_rms
+        phase_jitter = rng.standard_normal(N) * phase_noise_rms
         signal = A * np.sin(2*np.pi*Fin*t + phase_jitter) + DC
 
         # Measure jitter
diff --git a/python/tests/unit/aout/test_verify_spec_plot_phase.py b/python/tests/unit/aout/test_verify_spec_plot_phase.py
index ae22a82..a6448ff 100644
--- a/python/tests/unit/aout/test_verify_spec_plot_phase.py
+++ b/python/tests/unit/aout/test_verify_spec_plot_phase.py
@@ -21,8 +21,7 @@ def generate_test_signal():
         signal: Test signal with 4 harmonics + noise
         params: Dictionary of expected values
     """
-    # Set random seed for reproducibility
-    np.random.seed(42)
+    rng = np.random.default_rng(42)
 
     N = 8192           # Number of samples
     Fs = 1e9           # Sampling frequency
@@ -49,7 +48,7 @@ def generate_test_signal():
               A_HD4 * np.sin(2*np.pi*4*Fin*t + phi_HD4))
 
     # Add DC offset and small noise
-    signal = signal + 0.5 + np.random.randn(N) * 1e-5
+    signal = signal + 0.5 + rng.standard_normal(N) * 1e-5
 
     # Calculate maxSignal for normalization (needed for expected values)
     maxSignal = np.max(signal) - np.min(signal)
diff --git a/python/tests/unit/calibration/test_verify_calibration_full.py b/python/tests/unit/calibration/test_verify_calibration_full.py
index ddbb911..fb1cd87 100644
--- a/python/tests/unit/calibration/test_verify_calibration_full.py
+++ b/python/tests/unit/calibration/test_verify_calibration_full.py
@@ -62,7 +62,7 @@ def test_calibration_single_dataset_shuffled():
             shift_amounts_order = np.arange(bit_width - 1, -1, -1)
 
 
-            # shuffled_indices = np.random.permutation(bit_width)
+            # shuffled_indices = np.random.default_rng(2026062200).permutation(bit_width)
             shuffled_indices = np.arange(bit_width)
             shuffled_weights = true_weights[shuffled_indices]
             current_shifts = shift_amounts_order[shuffled_indices]
@@ -170,7 +170,7 @@ def test_calibration_single_dataset_shuffled_search_freq():
             shift_amounts_order = np.arange(bit_width - 1, -1, -1)
 
 
-            # shuffled_indices = np.random.permutation(bit_width)
+            # shuffled_indices = np.random.default_rng(2026062200).permutation(bit_width)
             shuffled_indices = np.arange(bit_width)
             shuffled_weights = true_weights[shuffled_indices]
             current_shifts = shift_amounts_order[shuffled_indices]
diff --git a/python/tests/unit/calibration/test_verify_lstsq_solver.py b/python/tests/unit/calibration/test_verify_lstsq_solver.py
index 0057bdd..329af11 100644
--- a/python/tests/unit/calibration/test_verify_lstsq_solver.py
+++ b/python/tests/unit/calibration/test_verify_lstsq_solver.py
@@ -142,12 +142,13 @@ def test_weight_recovery_absolute():
     4. Force the Solver's 'b' vector to match our Target Signal.
     """
     n_samples = 1000
+    rng = np.random.default_rng(2026062205)
     # 1. Define Ground Truth
     true_weights = np.array([0.25, 0.5, 1.0])
     bit_width = len(true_weights)
     
     # 2. Generate random bit matrix (n_samples, 3)
-    bits = np.random.randint(0, 2, (n_samples, bit_width)).astype(float)
+    bits = rng.integers(0, 2, (n_samples, bit_width)).astype(float)
     
     # 3. Construct the Target Signal directly from bits and weights
     # This signal is exactly what the solver should aim to reconstruct
@@ -238,15 +239,16 @@ def test_solve_weights_shared_recovery():
     Test 6: Verify that weights are shared across datasets but DCs are independent.
     """
     n_samples = 500
+    rng = np.random.default_rng(2026062206)
     true_weights = np.array([10.0, 20.0]) # Shared
     dc1, dc2 = 5.0, -3.0 # Different DCs
     
     # Dataset 1: bits @ weights + dc1
-    bits1 = np.random.randint(0, 2, (n_samples, 2)).astype(float)
+    bits1 = rng.integers(0, 2, (n_samples, 2)).astype(float)
     sig1 = bits1 @ true_weights + dc1
     
     # Dataset 2: bits @ weights + dc2
-    bits2 = np.random.randint(0, 2, (n_samples, 2)).astype(float)
+    bits2 = rng.integers(0, 2, (n_samples, 2)).astype(float)
     sig2 = bits2 @ true_weights + dc2
     
     # Mocking the basis to be the signals themselves to force direct recovery
diff --git a/python/tests/unit/calibration/test_verify_patch_rank_deficiency.py b/python/tests/unit/calibration/test_verify_patch_rank_deficiency.py
index a5ba574..29bb444 100644
--- a/python/tests/unit/calibration/test_verify_patch_rank_deficiency.py
+++ b/python/tests/unit/calibration/test_verify_patch_rank_deficiency.py
@@ -10,7 +10,8 @@ def test_patch_rank_deficiency_logic():
     """
     nominal_weights = [128, 64, 32, 16, 8, 4, 2, 1]
 
-    bits_input = np.random.randint(0, 2, (2048, len(nominal_weights)))
+    rng = np.random.default_rng(2026062207)
+    bits_input = rng.integers(0, 2, (2048, len(nominal_weights)))
     # Create dependencies:
     bits_input[:, 1] = 0
     bits_input[:, 2] = 1
@@ -54,7 +55,8 @@ def test_patch_rank_deficiency_logic_reverse():
     """
     nominal_weights = [128, 64, 32, 16, 8, 4, 2, 1][::-1]
 
-    bits_input = np.random.randint(0, 2, (2048, len(nominal_weights)))
+    rng = np.random.default_rng(2026062208)
+    bits_input = rng.integers(0, 2, (2048, len(nominal_weights)))
     # Create dependencies:
     bits_input[:, 1] = 0
     bits_input[:, 2] = 1
@@ -99,7 +101,8 @@ def test_patch_rank_deficiency_logic_recover():
     """
     nominal_weights = [128, 64, 32, 16, 8, 4, 2, 1][::-1]
 
-    bits_input = np.random.randint(0, 2, (2048, len(nominal_weights)))
+    rng = np.random.default_rng(2026062209)
+    bits_input = rng.integers(0, 2, (2048, len(nominal_weights)))
     # Create dependencies:
     bits_input[:, 1] = 0
     bits_input[:, 2] = 1
@@ -124,4 +127,4 @@ def test_patch_rank_deficiency_logic_recover():
 
     # --- Assertions ---
     assert weights_recovered[1] == 0.0
-    assert weights_recovered[2] == 0.0
\ No newline at end of file
+    assert weights_recovered[2] == 0.0
diff --git a/python/tests/unit/calibration/test_verify_prepare_input.py b/python/tests/unit/calibration/test_verify_prepare_input.py
index f2db145..fd52eb8 100644
--- a/python/tests/unit/calibration/test_verify_prepare_input.py
+++ b/python/tests/unit/calibration/test_verify_prepare_input.py
@@ -3,8 +3,10 @@
 
 def test_single_dataset():
     """Verify standard flow with N (samples) > M (bits)."""
+    rng = np.random.default_rng(2026062201)
+
     # --- Case 1: Proper Orientation ---
-    data = np.random.randint(0, 2, (1024, 8)).astype(float)
+    data = rng.integers(0, 2, (1024, 8)).astype(float)
     print()
     res = _prepare_input(data, verbose=2)
     
@@ -14,7 +16,7 @@ def test_single_dataset():
     np.testing.assert_array_equal(res["bits_segments"][0], data)
 
     # --- Case 2: Transpose Triggered (N < M) ---
-    data = np.random.randint(0, 2, (8, 1024)).astype(float)
+    data = rng.integers(0, 2, (8, 1024)).astype(float)
     print()
     res = _prepare_input(data, verbose=2)
     
@@ -24,7 +26,7 @@ def test_single_dataset():
     np.testing.assert_array_equal(res["bits_segments"][0], data.T)
 
     # --- Case 3: Square Matrix Boundary ---
-    data = np.random.randint(0, 2, (16, 16)).astype(float)
+    data = rng.integers(0, 2, (16, 16)).astype(float)
     print()
     res = _prepare_input(data, verbose=2)
     
@@ -36,11 +38,13 @@ def test_single_dataset():
 
 def test_multi_dataset():
     """Verify vertical stacking and segment tracking for multiple datasets."""
+    rng = np.random.default_rng(2026062202)
+
     # Create 4 datasets with varying sample counts but same bit width
-    d1 = np.random.randint(0, 2, (1024, 8)).astype(float)
-    d2 = np.random.randint(0, 2, (512, 8)).astype(float)
-    d3 = np.random.randint(0, 2, (8, 8192)).astype(float)
-    d4 = np.random.randint(0, 2, (8, 16384)).astype(float)
+    d1 = rng.integers(0, 2, (1024, 8)).astype(float)
+    d2 = rng.integers(0, 2, (512, 8)).astype(float)
+    d3 = rng.integers(0, 2, (8, 8192)).astype(float)
+    d4 = rng.integers(0, 2, (8, 16384)).astype(float)
     data = [d1, d2, d3, d4]
     print()
     res = _prepare_input(data, verbose=2)
diff --git a/python/tests/unit/calibration/test_verify_scale_columns_for_conditioning.py b/python/tests/unit/calibration/test_verify_scale_columns_for_conditioning.py
index b0f353f..e8f5f8a 100644
--- a/python/tests/unit/calibration/test_verify_scale_columns_for_conditioning.py
+++ b/python/tests/unit/calibration/test_verify_scale_columns_for_conditioning.py
@@ -8,7 +8,8 @@ def test_scaling_and_near_zero():
     # Col 1: Range [0, 1e6] (Large magnitude)
     # Col 2: Range [0, 0] (Dead bit / Near-zero)
     
-    bits = np.random.rand(1024, 12)    
+    rng = np.random.default_rng(2026062203)
+    bits = rng.random((1024, 12))
     bits[:, 3] *= 1e6  # Amplify col 3 to test large scaling
     bits[:, 5] *= 0    # Wipe col 5 to test near-zero protection
     bits[:, 11] *= 0    # Wipe col 11 to test near-zero protection
@@ -29,8 +30,9 @@ def test_scaling_and_near_zero():
 
 def test_scaling_and_recover():
     n_bits = 16
+    rng = np.random.default_rng(2026062204)
     nominal_weights = 2 ** np.arange(n_bits-1, -1, -1) / 2**(n_bits-1)
-    bits = np.random.randint(0, 2, (1024, n_bits)) * nominal_weights
+    bits = rng.integers(0, 2, (1024, n_bits)) * nominal_weights
 
     print()
     bits_effective, bit_scales = _scale_columns_for_conditioning(bits, verbose=2)
diff --git a/python/tests/unit/dout/test_plot_residual_scatter.py b/python/tests/unit/dout/test_plot_residual_scatter.py
index 128affc..33dc80a 100644
--- a/python/tests/unit/dout/test_plot_residual_scatter.py
+++ b/python/tests/unit/dout/test_plot_residual_scatter.py
@@ -8,9 +8,8 @@
 from adctoolbox import plot_residual_scatter
 
 
-def _make_adc_data(n=1024, m=6, seed=42):
+def _make_adc_data(n=1024, m=6):
     """Generate simple ADC signal and bit matrix for testing."""
-    np.random.seed(seed)
     t = np.arange(n)
     sig = (np.sin(2 * np.pi * 3 * t / n) / 2 + 0.5) * (2**m - 1)
     code = np.clip(np.round(sig).astype(int), 0, 2**m - 1)
@@ -18,6 +17,21 @@ def _make_adc_data(n=1024, m=6, seed=42):
     return sig, bits
 
 
+def _visible_axes(fig):
+    return [ax for ax in fig.axes if ax.get_visible()]
+
+
+def _hidden_axes(fig):
+    return [ax for ax in fig.axes if not ax.get_visible()]
+
+
+def _assert_scatter_axis(ax, *, n_points, xlabel, ylabel):
+    assert ax.get_xlabel() == xlabel
+    assert ax.get_ylabel() == ylabel
+    assert len(ax.collections) == 1
+    assert ax.collections[0].get_offsets().shape == (n_points, 2)
+
+
 def test_basic_output():
     """Basic call returns dict with expected keys."""
     sig, bits = _make_adc_data()
@@ -82,9 +96,43 @@ def test_residual_stage_zero():
 
 
 def test_plot_creation():
-    """Plot creates without error."""
+    """Plot creates one scatter axis per requested pair."""
     sig, bits = _make_adc_data()
-    result = plot_residual_scatter(sig, bits, pairs=[(0, 6), (3, 6)],
-                                   create_plot=True)
-    assert result is not None
-    plt.close('all')
+    pairs = [(0, 6), (3, 6)]
+    result = plot_residual_scatter(sig, bits, pairs=pairs, create_plot=True)
+    fig = plt.gcf()
+    visible_axes = _visible_axes(fig)
+
+    assert result['pairs'] == pairs
+    assert len(fig.axes) == 2
+    assert len(visible_axes) == len(pairs)
+    _assert_scatter_axis(
+        visible_axes[0],
+        n_points=len(sig),
+        xlabel='Signal',
+        ylabel='Res. of bit #6',
+    )
+    _assert_scatter_axis(
+        visible_axes[1],
+        n_points=len(sig),
+        xlabel='Res. of bit #3',
+        ylabel='Res. of bit #6',
+    )
+    plt.close(fig)
+
+
+def test_plot_creation_hides_unused_subplots():
+    """Unused axes in the subplot grid should be hidden."""
+    sig, bits = _make_adc_data(n=128)
+    pairs = [(0, 6), (1, 6), (2, 6), (3, 6), (4, 6)]
+    result = plot_residual_scatter(sig, bits, pairs=pairs, create_plot=True)
+    fig = plt.gcf()
+
+    assert result['pairs'] == pairs
+    assert len(fig.axes) == 8
+    assert len(_visible_axes(fig)) == len(pairs)
+    assert len(_hidden_axes(fig)) == 3
+    for ax in _visible_axes(fig):
+        assert len(ax.collections) == 1
+        assert ax.collections[0].get_offsets().shape == (len(sig), 2)
+    plt.close(fig)
diff --git a/python/tests/unit/fundamentals/test_verify_fit_sine_4param.py b/python/tests/unit/fundamentals/test_verify_fit_sine_4param.py
index 8b929a1..33b7841 100644
--- a/python/tests/unit/fundamentals/test_verify_fit_sine_4param.py
+++ b/python/tests/unit/fundamentals/test_verify_fit_sine_4param.py
@@ -54,7 +54,7 @@ def test_verify_fit_sine_4param_noisy_signal():
     2. Fit using fit_sine_4param
     3. Assert: Fitted parameters close to true values (within noise margin)
     """
-    np.random.seed(42)
+    rng = np.random.default_rng(42)
     N = 1000
     t = np.arange(N)
 
@@ -66,7 +66,7 @@ def test_verify_fit_sine_4param_noisy_signal():
 
     # Generate noisy signal
     sig_ideal = A_true * np.sin(2*np.pi*freq_true*t) + dc_true
-    noise = np.random.normal(0, noise_std, N)
+    noise = rng.normal(0, noise_std, N)
     sig_noisy = sig_ideal + noise
 
     # Fit the signal
diff --git a/python/tests/unit/oversampling/test_matlab_compat_api.py b/python/tests/unit/oversampling/test_matlab_compat_api.py
index 77cecca..e566ca8 100644
--- a/python/tests/unit/oversampling/test_matlab_compat_api.py
+++ b/python/tests/unit/oversampling/test_matlab_compat_api.py
@@ -55,8 +55,9 @@ def test_ifilter_rejects_invalid_frequency_bands():
 
 def test_perfosr_returns_matlab_order_and_matches_existing_core():
     n = 1024
+    rng = np.random.default_rng(0)
     t = np.arange(n)
-    data = 0.5 * np.sin(2 * np.pi * 0.05 * t) + 0.002 * np.random.RandomState(0).normal(size=n)
+    data = 0.5 * np.sin(2 * np.pi * 0.05 * t) + 0.002 * rng.normal(size=n)
     osr = np.array([2, 4, 8, 16])
 
     osr_out, sndr, sfdr, enob = perfosr(data, osr=osr, disp=False)
diff --git a/python/tests/unit/spectrum/test_align_spectrum_phase.py b/python/tests/unit/spectrum/test_align_spectrum_phase.py
index e096758..92716e8 100644
--- a/python/tests/unit/spectrum/test_align_spectrum_phase.py
+++ b/python/tests/unit/spectrum/test_align_spectrum_phase.py
@@ -63,9 +63,10 @@ def test_magnitude_preservation(self):
         n_fft = 256
         bin_idx = 25
         bin_r = 25.0
+        rng = np.random.default_rng(2026062223)
 
         # Create FFT data
-        fft_data = np.random.randn(n_fft) + 1j * np.random.randn(n_fft)
+        fft_data = rng.standard_normal(n_fft) + 1j * rng.standard_normal(n_fft)
         original_mag = np.abs(fft_data)
 
         fft_aligned = _align_spectrum_phase(fft_data, bin_idx, bin_r, n_fft)
@@ -178,8 +179,9 @@ def test_output_is_copy(self):
         n_fft = 128
         bin_idx = 10
         bin_r = 10.0
+        rng = np.random.default_rng(2026062224)
 
-        fft_data = np.random.randn(n_fft) + 1j * np.random.randn(n_fft)
+        fft_data = rng.standard_normal(n_fft) + 1j * rng.standard_normal(n_fft)
         fft_data_original = fft_data.copy()
 
         _align_spectrum_phase(fft_data, bin_idx, bin_r, n_fft)
@@ -229,8 +231,9 @@ def test_large_fft_size(self):
         n_fft = 16384
         bin_idx = 1000
         bin_r = 1000.0
+        rng = np.random.default_rng(2026062225)
 
-        fft_data = np.random.randn(n_fft) + 1j * np.random.randn(n_fft)
+        fft_data = rng.standard_normal(n_fft) + 1j * rng.standard_normal(n_fft)
         fft_data[bin_idx] = 100.0 * np.exp(1j * np.pi / 3)
 
         # Should complete without error
diff --git a/python/tests/unit/spectrum/test_assumed_signal_power.py b/python/tests/unit/spectrum/test_assumed_signal_power.py
index 59f36e8..43e32b5 100644
--- a/python/tests/unit/spectrum/test_assumed_signal_power.py
+++ b/python/tests/unit/spectrum/test_assumed_signal_power.py
@@ -19,9 +19,10 @@ def test_assumed_signal_power(assumed_sig_pwr_dbfs):
     Fs = 100e6
     Fin = 123 / N_fft * Fs  # Coherent frequency
     t = np.arange(N_fft) / Fs
+    rng = np.random.default_rng(2026062231)
 
     # Generate pure tone
-    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + np.random.randn(N_fft) * noise_rms
+    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + rng.standard_normal(N_fft) * noise_rms
 
     # Compute spectrum WITH assumed signal power
     result = compute_spectrum(signal, fs=Fs, win_type='hann', side_bin=1, assumed_sig_pwr_dbfs=assumed_sig_pwr_dbfs)
diff --git a/python/tests/unit/spectrum/test_compute_spectrum_peak_power.py b/python/tests/unit/spectrum/test_compute_spectrum_peak_power.py
index da216a2..11779e3 100644
--- a/python/tests/unit/spectrum/test_compute_spectrum_peak_power.py
+++ b/python/tests/unit/spectrum/test_compute_spectrum_peak_power.py
@@ -56,9 +56,10 @@ def test_peak_power_with_windows(win_type, side_bin):
     Fs = 100e6
     Fin = 123 / N_fft * Fs  # Coherent frequency
     t = np.arange(N_fft) / Fs
+    rng = np.random.default_rng(2026062232)
 
     # Generate pure tone
-    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + np.random.randn(N_fft) * 1e-6
+    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + rng.standard_normal(N_fft) * 1e-6
 
     # Compute spectrum
     result = compute_spectrum(signal, fs=Fs, win_type=win_type, side_bin=side_bin)
@@ -101,9 +102,10 @@ def test_peak_power_with_windows_max_scale_range(win_type, side_bin):
     Fs = 100e6
     Fin = 123 / N_fft * Fs  # Coherent frequency
     t = np.arange(N_fft) / Fs
+    rng = np.random.default_rng(2026062233)
 
     # Generate pure tone
-    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + np.random.randn(N_fft) * 1e-6
+    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + rng.standard_normal(N_fft) * 1e-6
 
     # Compute spectrum
     max_scale_range = [-1, 1]
diff --git a/python/tests/unit/spectrum/test_compute_spectrum_sndr.py b/python/tests/unit/spectrum/test_compute_spectrum_sndr.py
index 5344f2d..3f1df2a 100644
--- a/python/tests/unit/spectrum/test_compute_spectrum_sndr.py
+++ b/python/tests/unit/spectrum/test_compute_spectrum_sndr.py
@@ -150,8 +150,9 @@ def test_compute_spectrum_sndr(signal_amplitude, noise_rms):
     Fs = 100e6     # 100 MHz sampling rate
     Fin = 123 / N_fft * Fs  # Coherent frequency
     t = np.arange(N_fft) / Fs
+    rng = np.random.default_rng(2026062226)
 
-    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + np.random.randn(N_fft) * noise_rms
+    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + rng.standard_normal(N_fft) * noise_rms
 
     # Define full scale range as internal parameter
     max_scale_range = [-1, 1]
@@ -188,8 +189,9 @@ def test_compute_spectrum_sndr_no_max_scale_range(signal_amplitude, noise_rms):
     Fs = 100e6     # 100 MHz sampling rate
     Fin = 123 / N_fft * Fs  # Coherent frequency
     t = np.arange(N_fft) / Fs
+    rng = np.random.default_rng(2026062227)
 
-    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + np.random.randn(N_fft) * noise_rms
+    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + rng.standard_normal(N_fft) * noise_rms
 
     # Use hann window
     result = compute_spectrum(signal, fs=Fs, win_type='hann', verbose=0)
@@ -239,8 +241,9 @@ def test_compute_spectrum_window_sweep(win_type, side_bin):
     Fs = 100e6     # 100 MHz sampling rate
     Fin = 123 / N_fft * Fs  # Coherent frequency
     t = np.arange(N_fft) / Fs
+    rng = np.random.default_rng(2026062228)
 
-    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + np.random.randn(N_fft) * noise_rms
+    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + rng.standard_normal(N_fft) * noise_rms
 
     # Define full scale range as internal parameter
     max_scale_range = [-2, 2]
@@ -294,8 +297,9 @@ def test_compute_spectrum_window_sweep_no_max_scale_range(win_type, side_bin):
     Fs = 100e6     # 100 MHz sampling rate
     Fin = 123 / N_fft * Fs  # Coherent frequency
     t = np.arange(N_fft) / Fs
+    rng = np.random.default_rng(2026062229)
 
-    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + np.random.randn(N_fft) * noise_rms
+    signal = signal_amplitude * np.sin(2 * np.pi * Fin * t) + rng.standard_normal(N_fft) * noise_rms
 
     print(f"\n[Window Test] win_type={win_type}, side_bin={side_bin}, A={signal_amplitude:.2f}, noise_rms={noise_rms*1e6:.2f}uV")
 
diff --git a/python/tests/unit/spectrum/test_compute_spectrum_sndr_averaging.py b/python/tests/unit/spectrum/test_compute_spectrum_sndr_averaging.py
index b7f2dc1..7117499 100644
--- a/python/tests/unit/spectrum/test_compute_spectrum_sndr_averaging.py
+++ b/python/tests/unit/spectrum/test_compute_spectrum_sndr_averaging.py
@@ -23,8 +23,9 @@ def test_power_averaging_sndr(M):
     t = np.arange(N_fft) / Fs
     signal_clean = signal_amplitude * np.sin(2 * np.pi * Fin * t)
     signal_runs = np.zeros((M, N_fft))
+    rng = np.random.default_rng(2026062234)
     for m in range(M):
-        signal_runs[m, :] = signal_clean + np.random.randn(N_fft) * noise_rms
+        signal_runs[m, :] = signal_clean + rng.standard_normal(N_fft) * noise_rms
 
     # Power averaging
     result = compute_spectrum(signal_runs, fs=Fs, coherent_averaging=False)
@@ -68,8 +69,9 @@ def test_coherent_averaging_sndr(M):
     t = np.arange(N_fft) / Fs
     signal_clean = signal_amplitude * np.sin(2 * np.pi * Fin * t)
     signal_runs = np.zeros((M, N_fft))
+    rng = np.random.default_rng(2026062235)
     for m in range(M):
-        signal_runs[m, :] = signal_clean + np.random.randn(N_fft) * noise_rms
+        signal_runs[m, :] = signal_clean + rng.standard_normal(N_fft) * noise_rms
 
     # Coherent averaging
     result = compute_spectrum(signal_runs, fs=Fs, coherent_averaging=True)
diff --git a/python/tests/unit/spectrum/test_extract_harmonic_powers.py b/python/tests/unit/spectrum/test_extract_harmonic_powers.py
index fd8241f..434d87a 100644
--- a/python/tests/unit/spectrum/test_extract_harmonic_powers.py
+++ b/python/tests/unit/spectrum/test_extract_harmonic_powers.py
@@ -28,7 +28,8 @@ def test_harmonic_detection(hd2_target, hd3_target):
     # Generate signal with both HD2 and HD3 distortion
     t = np.arange(N_fft) / Fs
     sig_ideal = A * np.sin(2 * np.pi * Fin * t)
-    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + np.random.randn(N_fft) * noise_rms
+    rng = np.random.default_rng(2026062236)
+    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + rng.standard_normal(N_fft) * noise_rms
 
     # Run spectrum analysis
     result = compute_spectrum(signal, fs=Fs, max_harmonic=6, side_bin=1)
diff --git a/python/tests/unit/spectrum/test_locate_bins.py b/python/tests/unit/spectrum/test_locate_bins.py
index eb127fb..f84c7c5 100644
--- a/python/tests/unit/spectrum/test_locate_bins.py
+++ b/python/tests/unit/spectrum/test_locate_bins.py
@@ -19,7 +19,8 @@ def test_locate_bins(bin_target, side_bin, max_harmonic):
 
     f_sig = bin_target * fs / n_fft
     t = np.arange(n_fft) / fs
-    signal = np.sin(2 * np.pi * f_sig * t) + np.random.randn(n_fft) * noise_rms
+    rng = np.random.default_rng(2026062237)
+    signal = np.sin(2 * np.pi * f_sig * t) + rng.standard_normal(n_fft) * noise_rms
 
     result = compute_spectrum(signal, fs=fs, osr=osr, side_bin=side_bin, max_harmonic=max_harmonic)
 
@@ -78,7 +79,8 @@ def test_locate_bins_noncoherent(fin, max_harmonic):
 
     # Generate signal at non-coherent frequency
     t = np.arange(n_fft) / fs
-    signal = np.sin(2 * np.pi * fin * t) + np.random.randn(n_fft) * noise_rms
+    rng = np.random.default_rng(2026062238)
+    signal = np.sin(2 * np.pi * fin * t) + rng.standard_normal(n_fft) * noise_rms
 
     result = compute_spectrum(signal, fs=fs, osr=osr, side_bin=side_bin, max_harmonic=max_harmonic)
 
diff --git a/python/tests/unit/spectrum/test_nf_methods.py b/python/tests/unit/spectrum/test_nf_methods.py
index 3e3b077..8402b59 100644
--- a/python/tests/unit/spectrum/test_nf_methods.py
+++ b/python/tests/unit/spectrum/test_nf_methods.py
@@ -20,7 +20,8 @@ def test_nf_methods_comparison():
     Fin_coherent = bin_target * Fs / N_fft
 
     t = np.arange(N_fft) / Fs
-    signal = A * np.sin(2*np.pi*Fin_coherent*t) + np.random.randn(N_fft) * noise_rms
+    rng = np.random.default_rng(2026062239)
+    signal = A * np.sin(2*np.pi*Fin_coherent*t) + rng.standard_normal(N_fft) * noise_rms
 
     # MATLAB numbering: 0=auto, 1=median, 2=trimmed, 3=exclude
     results_auto = compute_spectrum(signal, fs=Fs, nf_method=0)
diff --git a/python/tests/unit/spectrum/test_noise_shaping_spectrum.py b/python/tests/unit/spectrum/test_noise_shaping_spectrum.py
index 4d16777..d588505 100644
--- a/python/tests/unit/spectrum/test_noise_shaping_spectrum.py
+++ b/python/tests/unit/spectrum/test_noise_shaping_spectrum.py
@@ -26,6 +26,7 @@ def test_first_order_noise_shaping_spectrum():
     OSR = 64  # Oversampling ratio
     Fin_target = Fs / (2 * OSR) / 14  # Signal in lower quarter of Nyquist band
     A = 0.49
+    rng = np.random.default_rng(2026062230)
 
     # Find coherent frequency
     Fin, Fin_bin = find_coherent_frequency(fs=Fs, fin_target=Fin_target, n_fft=N)
@@ -35,10 +36,10 @@ def test_first_order_noise_shaping_spectrum():
     signal_ideal = A * np.sin(2 * np.pi * Fin * t)
 
     # Thermal noise (white, cannot be shaped) - keep small
-    thermal_noise = np.random.randn(N) * 100e-6
+    thermal_noise = rng.standard_normal(N) * 100e-6
 
     # Quantization noise (white noise as approximation, can be shaped) - make dominant
-    quant_noise_white = np.random.randn(N) * 2000e-6
+    quant_noise_white = rng.standard_normal(N) * 2000e-6
 
     # Apply 1st order noise shaping to quantization noise: NTF(z) = 1 - z^-1
     # noise_shaped[n] = noise[n] - noise[n-1]
@@ -68,8 +69,8 @@ def test_first_order_noise_shaping_spectrum():
 
     for run in range(N_runs):
         # Generate new noise for each run
-        thermal_run = np.random.randn(N) * 100e-6
-        quant_run = np.random.randn(N) * 2000e-6
+        thermal_run = rng.standard_normal(N) * 100e-6
+        quant_run = rng.standard_normal(N) * 2000e-6
 
         # Apply noise shaping to quantization noise
         quant_shaped_1st_run = np.zeros(N)
diff --git a/python/tests/unit/spectrum/test_plot_spectrum.py b/python/tests/unit/spectrum/test_plot_spectrum.py
index 9d383b0..3032594 100644
--- a/python/tests/unit/spectrum/test_plot_spectrum.py
+++ b/python/tests/unit/spectrum/test_plot_spectrum.py
@@ -17,6 +17,41 @@
 output_dir.mkdir(exist_ok=True)
 
 
+def _text_values(ax):
+    return [text.get_text() for text in ax.texts]
+
+
+def _assert_spectrum_axes(ax, result, *, show_label=True, expected_title='Power Spectrum'):
+    """Verify that plot_spectrum rendered the expected plot structure."""
+    spectrum_line = ax.lines[0]
+    assert len(spectrum_line.get_xdata()) == len(result['plot_data']['freq'])
+    assert len(spectrum_line.get_ydata()) == len(result['plot_data']['power_spectrum_db_plot'])
+
+    ymin, ymax = ax.get_ylim()
+    xmin, xmax = ax.get_xlim()
+    assert xmin > 0
+    assert xmax == pytest.approx(result['fs'] / 2)
+    assert ymin < 0
+    assert ymax == 0
+
+    if expected_title is not None:
+        assert ax.get_title() == expected_title
+
+    if show_label:
+        assert ax.get_xlabel() == 'Freq (Hz)'
+        assert ax.get_ylabel() == 'dBFS'
+        texts = _text_values(ax)
+        for prefix in ('Fin/fs =', 'ENoB =', 'SNDR =', 'Sig ='):
+            assert any(text.startswith(prefix) for text in texts)
+        assert 'MaxSpur' in texts
+        assert any(line.get_marker() == 'o' for line in ax.lines)
+        assert any(line.get_marker() == 'd' for line in ax.lines)
+    else:
+        assert ax.get_xlabel() == ''
+        assert ax.get_ylabel() == ''
+        assert _text_values(ax) == []
+
+
 @pytest.mark.parametrize(
     "nf_line_level,expected",
     [
@@ -132,7 +167,8 @@ def test_plot_spectrum_distorted_sine(hd2_target, hd3_target):
     # Generate signal with both HD2 and HD3 distortion
     t = np.arange(N_fft) / Fs
     sig_ideal = A * np.sin(2 * np.pi * Fin * t)
-    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + np.random.randn(N_fft) * noise_rms
+    rng = np.random.default_rng(2026062246)
+    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + rng.standard_normal(N_fft) * noise_rms
 
     # Compute spectrum
     result = compute_spectrum(signal, fs=Fs, max_harmonic=6, side_bin=1)
@@ -140,12 +176,14 @@ def test_plot_spectrum_distorted_sine(hd2_target, hd3_target):
     # Create figure and plot
     fig, ax = plt.subplots(figsize=(10, 6))
     plot_spectrum(result, show_title=True, show_label=True, plot_harmonics_up_to=3, ax=ax)
+    _assert_spectrum_axes(ax, result)
+    assert len(ax.lines) >= 5
 
     # Save figure
     fig_path = output_dir / f'test_plot_distorted_hd2_{hd2_target}_hd3_{hd3_target}.png'
     plt.tight_layout()
     plt.savefig(fig_path, dpi=150, bbox_inches='tight')
-    plt.close()
+    plt.close(fig)
 
     # Verify file was created
     assert fig_path.exists(), f"Figure file not created: {fig_path}"
@@ -165,7 +203,8 @@ def test_plot_spectrum_clean_sine():
 
     # Generate clean signal
     t = np.arange(N_fft) / Fs
-    signal = A * np.sin(2 * np.pi * Fin * t) + np.random.randn(N_fft) * noise_rms
+    rng = np.random.default_rng(2026062247)
+    signal = A * np.sin(2 * np.pi * Fin * t) + rng.standard_normal(N_fft) * noise_rms
 
     # Compute spectrum
     result = compute_spectrum(signal, fs=Fs, max_harmonic=5, side_bin=1)
@@ -173,12 +212,13 @@ def test_plot_spectrum_clean_sine():
     # Create figure and plot
     fig, ax = plt.subplots(figsize=(10, 6))
     plot_spectrum(result, show_title=True, show_label=True, plot_harmonics_up_to=3, ax=ax)
+    _assert_spectrum_axes(ax, result)
 
     # Save figure
     fig_path = output_dir / 'test_plot_clean_sine.png'
     plt.tight_layout()
     plt.savefig(fig_path, dpi=150, bbox_inches='tight')
-    plt.close()
+    plt.close(fig)
 
     # Verify file was created
     assert fig_path.exists(), f"Figure file not created: {fig_path}"
@@ -214,7 +254,8 @@ def test_plot_spectrum_comparison():
         # Generate signal
         t = np.arange(N_fft) / Fs
         sig_ideal = A * np.sin(2 * np.pi * Fin * t)
-        signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + np.random.randn(N_fft) * noise_rms
+        rng = np.random.default_rng(2026062248 + idx)
+        signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + rng.standard_normal(N_fft) * noise_rms
 
         # Compute spectrum
         result = compute_spectrum(signal, fs=Fs, max_harmonic=6, side_bin=1)
@@ -224,14 +265,17 @@ def test_plot_spectrum_comparison():
         plot_spectrum(result, show_title=False, show_label=True, plot_harmonics_up_to=3, ax=axes[idx])
         axes[idx].set_title(f'{config["label"]}: HD2={config["hd2"]} dBc, HD3={config["hd3"]} dBc',
                            fontsize=12, fontweight='bold')
+        _assert_spectrum_axes(axes[idx], result, expected_title=None)
+        assert config['label'] in axes[idx].get_title()
 
     # Save figure
     fig_path = output_dir / 'test_plot_comparison.png'
     plt.tight_layout()
     plt.savefig(fig_path, dpi=150, bbox_inches='tight')
-    plt.close()
+    plt.close(fig)
 
     # Verify file was created
+    assert len(axes) == len(distortion_configs)
     assert fig_path.exists(), f"Figure file not created: {fig_path}"
     assert fig_path.stat().st_size > 0, f"Figure file is empty: {fig_path}"
 
@@ -256,12 +300,14 @@ def test_plot_spectrum_no_labels():
     # Create figure and plot without labels
     fig, ax = plt.subplots(figsize=(10, 6))
     plot_spectrum(result, show_title=False, show_label=False, plot_harmonics_up_to=0, ax=ax)
+    _assert_spectrum_axes(ax, result, show_label=False, expected_title='')
+    assert len(ax.lines) == 1
 
     # Save figure
     fig_path = output_dir / 'test_plot_no_labels.png'
     plt.tight_layout()
     plt.savefig(fig_path, dpi=150, bbox_inches='tight')
-    plt.close()
+    plt.close(fig)
 
     # Verify file was created
     assert fig_path.exists(), f"Figure file not created: {fig_path}"
@@ -295,12 +341,16 @@ def test_plot_spectrum_high_harmonics():
     # Create figure and plot with many harmonics
     fig, ax = plt.subplots(figsize=(12, 7))
     plot_spectrum(result, show_title=True, show_label=True, plot_harmonics_up_to=7, ax=ax)
+    _assert_spectrum_axes(ax, result)
+    texts = set(_text_values(ax))
+    assert {'2', '3', '4'}.issubset(texts)
+    assert sum(line.get_marker() == 's' for line in ax.lines) >= 3
 
     # Save figure
     fig_path = output_dir / 'test_plot_high_harmonics.png'
     plt.tight_layout()
     plt.savefig(fig_path, dpi=150, bbox_inches='tight')
-    plt.close()
+    plt.close(fig)
 
     # Verify file was created
     assert fig_path.exists(), f"Figure file not created: {fig_path}"
diff --git a/python/tests/unit/spectrum/test_prepare_fft_input.py b/python/tests/unit/spectrum/test_prepare_fft_input.py
index 2f414fc..0609582 100644
--- a/python/tests/unit/spectrum/test_prepare_fft_input.py
+++ b/python/tests/unit/spectrum/test_prepare_fft_input.py
@@ -27,7 +27,8 @@ def test_prepare_fft_input_transpose_handling(input_shape, expected_output_shape
     Verifies correct shape transformation for various input formats.
     """
 
-    test_data = np.random.randn(*input_shape)
+    rng = np.random.default_rng(2026062240)
+    test_data = rng.standard_normal(input_shape)
     processed = _prepare_fft_input(test_data)
 
     assert processed.shape == expected_output_shape
@@ -53,7 +54,8 @@ def test_prepare_fft_input_transpose_handling(input_shape, expected_output_shape
 ])
 def test_prepare_fft_input_max_scale_range(signal_range, max_scale_range, expected_range):
     """Test max_scale_range normalization."""
-    test_signal = np.random.uniform(*signal_range, 10000)
+    rng = np.random.default_rng(2026062241)
+    test_signal = rng.uniform(*signal_range, 10000)
     processed = _prepare_fft_input(test_signal, max_scale_range=max_scale_range)
 
     assert np.abs(np.mean(processed)) < 1e-10
diff --git a/python/tests/unit/spectrum/test_spectrum_averaging.py b/python/tests/unit/spectrum/test_spectrum_averaging.py
index 844d289..0e197eb 100644
--- a/python/tests/unit/spectrum/test_spectrum_averaging.py
+++ b/python/tests/unit/spectrum/test_spectrum_averaging.py
@@ -15,7 +15,8 @@
 ])
 def test_power_average(M, N):
     """Test power-averaged spectrum computation."""
-    data = np.random.randn(M, N)
+    rng = np.random.default_rng(2026062242)
+    data = rng.standard_normal((M, N))
 
     spectrum_power, _ = _power_average(data)
 
@@ -64,7 +65,8 @@ def test_coherent_average(M, N, osr):
     bin_target = N // 15
     A = 1.0  # Sine wave amplitude
     # Vectorized sine wave generation with random phases
-    phases = 2 * np.pi * np.random.rand(M, 1)
+    rng = np.random.default_rng(2026062243)
+    phases = 2 * np.pi * rng.random((M, 1))
     t = np.arange(N)
     data = A * np.sin(2 * np.pi * bin_target * t / N + phases)
 
diff --git a/python/tests/unit/spectrum/test_spectrum_values.py b/python/tests/unit/spectrum/test_spectrum_values.py
index b7334c6..0776761 100644
--- a/python/tests/unit/spectrum/test_spectrum_values.py
+++ b/python/tests/unit/spectrum/test_spectrum_values.py
@@ -28,7 +28,8 @@ def test_spectrum_values_at_bins(hd2_target, hd3_target):
     # Generate signal with both HD2 and HD3 distortion
     t = np.arange(N_fft) / Fs
     sig_ideal = A * np.sin(2 * np.pi * Fin * t)
-    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + np.random.randn(N_fft) * noise_rms
+    rng = np.random.default_rng(2026062244)
+    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + rng.standard_normal(N_fft) * noise_rms
 
     # Run spectrum analysis
     result = compute_spectrum(signal, fs=Fs, max_harmonic=6, side_bin=1)
@@ -93,7 +94,8 @@ def test_spectrum_values_at_bins_max_scale_range(hd2_target, hd3_target):
     # Generate signal with both HD2 and HD3 distortion
     t = np.arange(N_fft) / Fs
     sig_ideal = A * np.sin(2 * np.pi * Fin * t)
-    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + np.random.randn(N_fft) * noise_rms
+    rng = np.random.default_rng(2026062245)
+    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + rng.standard_normal(N_fft) * noise_rms
 
     # Run spectrum analysis
     result = compute_spectrum(signal, fs=Fs, max_harmonic=6, side_bin=1, max_scale_range=[-1, 1])
diff --git a/python/tests/unit/spectrum/test_sweep_performance_vs_osr.py b/python/tests/unit/spectrum/test_sweep_performance_vs_osr.py
index 9f78c43..a794005 100644
--- a/python/tests/unit/spectrum/test_sweep_performance_vs_osr.py
+++ b/python/tests/unit/spectrum/test_sweep_performance_vs_osr.py
@@ -9,6 +9,27 @@
 from adctoolbox.spectrum._bin_ranges import rfft_inband_bin_count
 
 
+def _legend_labels(ax):
+    legend = ax.get_legend()
+    assert legend is not None
+    return [text.get_text() for text in legend.get_texts()]
+
+
+def _assert_performance_axis(ax):
+    assert ax.get_xlabel() == 'OSR'
+    assert ax.get_ylabel() == 'SNDR / SFDR (dB)'
+    assert ax.get_title() == 'Performance vs OSR'
+    assert ax.get_xscale() == 'log'
+    assert len(ax.lines) >= 2
+    assert {'SNDR (ENOB)', 'SFDR'}.issubset(_legend_labels(ax))
+
+
+def _assert_enob_axis(ax):
+    assert ax.get_ylabel() == 'ENOB (bits)'
+    ymin, ymax = ax.get_ylim()
+    assert ymin < ymax
+
+
 def test_clean_sine_high_sndr():
     """Clean sine: SNDR should be very high at all OSR."""
     N = 1024
@@ -33,12 +54,12 @@ def test_clean_sine_high_sndr():
 
 def test_noisy_sine_sndr_increases_with_osr():
     """Sine + white noise: SNDR should increase with OSR (~3 dB per doubling)."""
-    np.random.seed(42)
+    rng = np.random.default_rng(42)
     N = 4096
     t = np.arange(N)
     freq = 0.05
     sig = 0.5 * np.sin(2 * np.pi * freq * t)
-    noise = np.random.normal(0, 0.01, N)
+    noise = rng.normal(0, 0.01, N)
     data = sig + noise
 
     osr_values = np.array([2, 4, 8, 16, 32, 64])
@@ -92,22 +113,46 @@ def test_plot_with_ax():
     N = 512
     t = np.arange(N)
     sig = 0.5 * np.sin(2 * np.pi * 0.1 * t)
+    osr = np.array([2, 4, 8])
 
     fig, ax = plt.subplots()
-    result = sweep_performance_vs_osr(sig, osr=np.array([2, 4, 8]), ax=ax)
-    assert result is not None
-    plt.close('all')
+    result = sweep_performance_vs_osr(sig, osr=osr, ax=ax)
+
+    np.testing.assert_array_equal(result['osr'], osr)
+    assert len(result['sndr']) == len(osr)
+    assert fig.axes[0] is ax
+    assert len(fig.axes) == 2
+    _assert_performance_axis(ax)
+    _assert_enob_axis(fig.axes[1])
+    plt.close(fig)
 
 
 def test_plot_auto_subplots():
     """No ax provided should create 2-subplot figure."""
+    plt.close('all')
     N = 512
+    rng = np.random.default_rng(0)
     t = np.arange(N)
-    sig = 0.5 * np.sin(2 * np.pi * 0.1 * t) + np.random.RandomState(0).normal(0, 0.01, N)
+    sig = 0.5 * np.sin(2 * np.pi * 0.1 * t) + rng.normal(0, 0.01, N)
 
     result = sweep_performance_vs_osr(sig, osr=np.array([2, 4, 8, 16]))
-    assert result is not None
-    plt.close('all')
+    fig = plt.gcf()
+    main_axes = [ax for ax in fig.axes if ax.get_title() == 'Performance vs OSR']
+    slope_axes = [ax for ax in fig.axes if ax.get_ylabel() == 'SNDR Slope (dB/decade)']
+    enob_axes = [ax for ax in fig.axes if ax.get_ylabel() == 'ENOB (bits)']
+
+    assert len(result['osr']) == 4
+    assert len(fig.axes) == 3
+    assert len(main_axes) == 1
+    assert len(slope_axes) == 1
+    assert len(enob_axes) == 1
+    _assert_performance_axis(main_axes[0])
+    _assert_enob_axis(enob_axes[0])
+    assert slope_axes[0].get_xlabel() == 'OSR'
+    assert slope_axes[0].get_xscale() == 'log'
+    assert len(slope_axes[0].lines) >= 2
+    assert any(text.get_text() == 'White Noise Limit' for text in slope_axes[0].texts)
+    plt.close(fig)
 
 
 def test_enob_formula():
diff --git a/python/tests/unit/spectrum/test_verify_compute_spectrum.py b/python/tests/unit/spectrum/test_verify_compute_spectrum.py
index ff3377d..5c6e2a6 100644
--- a/python/tests/unit/spectrum/test_verify_compute_spectrum.py
+++ b/python/tests/unit/spectrum/test_verify_compute_spectrum.py
@@ -91,7 +91,8 @@ def _run_noise_accuracy_test(A, noise_rms, Fs, Fin_target, max_scale_range=None,
 
     # Generate signal with known noise
     t = np.arange(N) / Fs
-    signal = A * np.sin(2*np.pi*Fin*t) + np.random.randn(N) * noise_rms
+    rng = np.random.default_rng(2026062249)
+    signal = A * np.sin(2*np.pi*Fin*t) + rng.standard_normal(N) * noise_rms
 
     # Compute spectrum (use boxcar window for accurate NSD comparison)
     # Boxcar has ENBW=1.0, so NSD calculation matches theoretical formula
@@ -291,7 +292,8 @@ def test_verify_compute_spectrum_2d_input():
     2. Assert: Averages correctly over runs
     """
     M, N = 3, 512
-    data_2d = 0.4 * np.sin(2*np.pi*0.1*np.arange(N)[np.newaxis, :] + np.random.randn(M, 1) * 0.01)
+    rng = np.random.default_rng(2026062250)
+    data_2d = 0.4 * np.sin(2*np.pi*0.1*np.arange(N)[np.newaxis, :] + rng.standard_normal((M, 1)) * 0.01)
 
     result = compute_spectrum(data_2d, fs=1.0, verbose=0)
 
@@ -342,7 +344,8 @@ def test_verify_compute_spectrum_coherent_vs_power():
     """
     M, N = 2, 512
     t = np.arange(N) / 1000.0
-    data = 0.4 * np.sin(2*np.pi*100*t)[np.newaxis, :] + np.random.randn(M, N) * 0.01
+    rng = np.random.default_rng(2026062251)
+    data = 0.4 * np.sin(2*np.pi*100*t)[np.newaxis, :] + rng.standard_normal((M, N)) * 0.01
     data = np.vstack([data, data])  # 4 runs
 
     result_power = compute_spectrum(data, fs=1000.0, coherent_averaging=False, verbose=0)
diff --git a/python/tests/unit/spectrum/test_window_sweep.py b/python/tests/unit/spectrum/test_window_sweep.py
index 332d277..d480d46 100644
--- a/python/tests/unit/spectrum/test_window_sweep.py
+++ b/python/tests/unit/spectrum/test_window_sweep.py
@@ -56,7 +56,8 @@ def test_window_sweep_coherent():
     # Generate signal with harmonics and noise
     t = np.arange(N_fft) / Fs
     sig_ideal = A * np.sin(2*np.pi*Fin*t)
-    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + np.random.randn(N_fft) * noise_rms
+    rng = np.random.default_rng(2026062252)
+    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + rng.standard_normal(N_fft) * noise_rms
 
     print('\n' + '='*80)
     print('COHERENT SIGNAL - WINDOW SWEEP (Default side_bin)')
@@ -147,7 +148,8 @@ def test_window_sweep_noncoherent():
     # Generate signal with harmonics and noise
     t = np.arange(N_fft) / Fs
     sig_ideal = A * np.sin(2*np.pi*Fin*t)
-    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + np.random.randn(N_fft) * noise_rms
+    rng = np.random.default_rng(2026062253)
+    signal = sig_ideal + k2 * sig_ideal**2 + k3 * sig_ideal**3 + rng.standard_normal(N_fft) * noise_rms
 
     print('\n' + '='*80)
     print('NON-COHERENT SIGNAL - WINDOW SWEEP (Default side_bin)')
@@ -250,7 +252,8 @@ def test_noise_accuracy_with_windows(win_type):
     # Generate signal with known noise
     t = np.arange(N_fft) / Fs
     sig_ideal = A * np.sin(2*np.pi*Fin*t)
-    signal = sig_ideal + np.random.randn(N_fft) * noise_rms
+    rng = np.random.default_rng(2026062254)
+    signal = sig_ideal + rng.standard_normal(N_fft) * noise_rms
 
     # Calculate theoretical values
     sig_pwr_theory = 0.0  # dBFS (auto-detect uses peak as reference)