diff --git a/python/packages/isce3/signal/compute_evd_cpi.py b/python/packages/isce3/signal/compute_evd_cpi.py index fea0191ff..f525ec376 100644 --- a/python/packages/isce3/signal/compute_evd_cpi.py +++ b/python/packages/isce3/signal/compute_evd_cpi.py @@ -166,6 +166,9 @@ def compute_evd_tb( eig_val_sort_array = np.zeros([num_cpi, cpi_len], dtype="f4") eig_vec_sort_array = np.zeros((num_cpi, cpi_len, cpi_len), dtype="complex64") + diag_power_array = np.zeros((num_cpi, cpi_len), dtype=np.float32) + diag_valid_array = np.zeros((num_cpi, cpi_len), dtype=bool) + tb_is_valid = True # Compute Eigenvalue and Eigenvector pairs for each CPI @@ -174,7 +177,7 @@ def compute_evd_tb( ): data_cpi = raw_data[cpi_slow_time] mask_valid_cpi = None if mask_valid is None else mask_valid[cpi_slow_time] - eig_val_sort, eig_vec_sort = compute_evd( + eig_val_sort, eig_vec_sort, diag_power_cpi, diag_valid_idx = compute_evd( data_cpi, mask_valid_cpi=mask_valid_cpi, off_diag_overlap_ratio=off_diag_overlap_ratio, @@ -182,17 +185,27 @@ def compute_evd_tb( ) # Verify if the eigenvalue of CPI at index defind by min_ev_valid_idx is meaningful - eig_val_sort_abs = np.maximum(np.abs(eig_val_sort), 1e-30) + eig_val_sort_abs = np.abs(eig_val_sort) + + # If any eigenvalue is NaN/inf OR leading eigenvalue is non-positive → invalid TB + if (not np.all(np.isfinite(eig_val_sort_abs))) or (eig_val_sort_abs[0] <= 0): + tb_is_valid = False + break + + eig_val_sort_abs = np.maximum(eig_val_sort_abs, 1e-30) noise_ev_norm_db = 10 * np.log10(eig_val_sort_abs[min_ev_valid_idx] / eig_val_sort_abs[0]) - if -noise_ev_norm_db > rx_dynamic_range_db: + if (not np.isfinite(noise_ev_norm_db)) or -noise_ev_norm_db > rx_dynamic_range_db: tb_is_valid = False break eig_val_sort_array[idx_cpi] = eig_val_sort eig_vec_sort_array[idx_cpi] = eig_vec_sort - return eig_val_sort_array, eig_vec_sort_array, tb_is_valid + diag_power_array[idx_cpi] = diag_power_cpi + diag_valid_array[idx_cpi] = diag_valid_idx + + return eig_val_sort_array, eig_vec_sort_array, diag_power_array, diag_valid_array, tb_is_valid def compute_evd( raw_data: np.ndarray, @@ -234,6 +247,7 @@ def compute_evd( # Application in Narrow-Band Interference Suppression for SAR”, IEEE Geoscience # and Remote Sensing Letters, vol. 4, no. 1, pp. 76,2007. + if mask_valid_cpi is not None: mask_valid_cpi = mask_valid_cpi.astype(bool, copy=False) @@ -242,7 +256,7 @@ def compute_evd( f"Valid CPI mask shape {mask_valid_cpi.shape} != CPI data shape {raw_data.shape}" ) - cov_cpi = compute_gap_exclusion_cov( + cov_cpi, diag_valid_idx = compute_gap_exclusion_cov( raw_data, mask_valid_cpi=mask_valid_cpi, off_diag_overlap_ratio=off_diag_overlap_ratio, @@ -251,10 +265,13 @@ def compute_evd( else: num_rng_samples = raw_data.shape[1] cov_cpi = (raw_data @ raw_data.conj().T) / num_rng_samples + diag_valid_idx = np.ones(raw_data.shape[0], dtype=bool) + + diag_power_cpi = np.real(np.diag(cov_cpi)) eig_val_sort, eig_vec_sort = eigen_decomp_sort(cov_cpi) - return eig_val_sort, eig_vec_sort + return eig_val_sort, eig_vec_sort, diag_power_cpi, diag_valid_idx def compute_gap_exclusion_cov( @@ -329,7 +346,7 @@ def compute_gap_exclusion_cov( if min_valid_diag < rng_samples_min: warnings.warn(f""" Minimum number of samples required per pulse to estimate sample covariance matrix - is {rng_samples_min}. The number of valid diagonal samples is {min_valid_off_diag}. + is {rng_samples_min}. The number of valid diagonal samples is {min_valid_diag}. """) # Zero-out invalid samples @@ -379,4 +396,4 @@ def compute_gap_exclusion_cov( # Ensure Hermitian numerically cov = (0.5 * (cov + cov.conj().T)).astype(np.complex64) - return cov + return cov, diag_valid_idx diff --git a/python/packages/isce3/signal/rfi_detection_evd.py b/python/packages/isce3/signal/rfi_detection_evd.py index 82cfd23cd..4432384f1 100644 --- a/python/packages/isce3/signal/rfi_detection_evd.py +++ b/python/packages/isce3/signal/rfi_detection_evd.py @@ -1,6 +1,6 @@ """ -Performs RFI detection of input data using Slow-Time Eigenvalue Slope -Thresholding algorithm (ST-EST). +RFI Detection using DUAL-CHECK: Condition Number + Power Stationarity +Enhanced version - uses BOTH condition number variability AND power spread check """ import numpy as np from isce3.signal.compute_evd_cpi import compute_evd_tb @@ -10,7 +10,7 @@ @dataclass class ThresholdParams: - """This dataclass computes the interpolated value of the number + """This dataclass computes the interpolated value of the number of sigmas (standard deviation) of the first difference of minimum Eigenvalues across threshold block. @@ -28,7 +28,6 @@ class ThresholdParams: final threshold. The values of x outside of the defined range of y are extrapolated. Defaults are [5.0, 2.0] - """ x: List[float] = field(default_factory=lambda: [2.0, 20.0]) y: List[float] = field(default_factory=lambda: [5.0, 2.0]) @@ -39,13 +38,13 @@ def __post_init__(self) -> None: if len(self.x) < 2: raise ValueError("At least two points are required") - def rfi_detect( raw_data, cpi_len, max_deg_freedom, min_ev_valid_idx, *, + num_rfi_buffer=3, num_max_trim=0, num_min_trim=0, max_num_rfi_ev=2, @@ -53,12 +52,34 @@ def rfi_detect( diag_valid_ratio=0.20, rx_dynamic_range_db=50.0, mask_valid=None, + threshold_method='max_ev', threshold_params: ThresholdParams = ThresholdParams(), + rfi_fig_merit_thresh=1.0, + bright_target_check=True, + bt_condition_num_thresh_db=2.0, + bt_pwr_spread_thresh_db=6.0, + pwr_ref_percentile=50, + pwr_upper_percentile=99.5, + sig_ev_margin_upper_db=3.0, + sig_ev_margin_lower_db=0.0, + pcr_range=[0.1, 0.4], ): - """This wrapper performs Eigenvalue Decomposition of input raw data as well as + """This wrapper performs Eigenvalue Decomposition of input raw data as well as RFI Eigenvalue slope threshold estimation and detection. + DUAL-CHECK VERSION: Uses condition number check AND power spread check + + Workflow: + 1. Perform EVD on all CPIs in the threshold block + 2. Validate TB (check for sufficient usable eigenvalues) + 3. DUAL-CHECK to filter out bright targets: + a. Condition number variability check + b. Power spread check (upper tail spread of diagonal power) + SKIP if BOTH checks pass (stable eigenvalue structure + low power spread) + 4. If TB fails either check, apply selected threshold method (ev_slope/max_ev) + 5. Detect RFI eigenvalues per CPI based on computed threshold + Parameters ------------ raw_data: array-like complex [num_pulses x num_rng_samples] @@ -72,14 +93,18 @@ def rfi_detect( Eigenvalue index used by threshold estimation to estimate the slow-time minimum Eigenvalue slope. This parameter is also used to validate that the threshold block has enough usable eigenvalues for robust sample covaraince estimation of a CPI. + num_rfi_buffer: int, default=3 + Number of buffer eigenvalue indices to skip after last possible RFI EV before + starting clean segment interpolation for 'max_ev' method. The clean segment + starts at index (max_deg_freedom + num_rfi_buffer - 1). num_max_trim: int, default=0 Number of large-value outliers to be trimmed in slow-time minimum Eigenvalues. num_min_trim: int, default=0 Number of small-value outliers to be trimmed in slow-time minimum Eigenvalues max_num_rfi_ev: int, default=2 - A detection error (miss) happens when a maximum power RFI emitter contaminates - multiple consecutive CPIs, resulting in a flat maximum Eigenvalue slope in slow - time. Hence the standard (STD) deviation of multiple dominant EVs across slow time + A detection error (miss) happens when a maximum power RFI emitter contaminates + multiple consecutive CPIs, resulting in a flat maximum Eigenvalue slope in slow + time. Hence the standard (STD) deviation of multiple dominant EVs across slow time defined by this parameter are compared. The one with the maximum STD is used for RFI Eigenvalue first difference computation. off_diag_overlap_ratio : float, default=0.25 @@ -95,11 +120,60 @@ def rfi_detect( Valid-sample mask with same shape as raw_data If provided, it has the shape of [num_pulses x num_rng_samples]. CPI sample covariance matrix will be normalized differently by excluding the invalid data gaps. + threshold_method : str, default='max_ev' + Detection method: 'ev_slope' or 'max_ev' threshold_params: ThresholdParams dataclass object, default=ThresholdParams() RFI detection threshold interpolation parameters. The x field defines STD ratio between maximum and minimum Eigenvalue slopes (MMES) of the slow-time threshold interval. The y field defines the number of sigma (STD) from the mean of MMES. + rfi_fig_merit_thresh : float or None, default=1.0 + Figure of merit threshold for 'max_ev' method. Default of 1.0 ensures + aggressive RFI detection - checks as many TBs as possible. + bright_target_check : bool, default=True + Enable dual bright target rejection check using BOTH condition number variability + AND power spread. When True, TB is skipped only if BOTH checks pass: + - Condition number std <= bt_condition_num_thresh_db (stable eigenvalue structure) + - Upper tail spread <= bt_pwr_spread_thresh_db (low power spread) + When False, bright target checks are disabled (all TBs are processed). + bt_condition_num_thresh_db : float, default=2.0 + Condition number std threshold in dB. TBs with std(cond#) <= this value + pass the condition# check (stable eigenvalue structure). + Only used when bright_target_check=True. + bt_pwr_spread_thresh_db : float, default=6.0 + Power spread threshold in dB. TBs with upper tail spread <= this value + pass the power check (low power spread across slow time). + Upper tail spread = pwr_upper_percentile - pwr_ref_percentile of diagonal power. + Only used when bright_target_check=True. + pwr_ref_percentile : float, default=50 + Reference (baseline) percentile for power spread computation. + Anchors the lower end of the spread measurement. Only used when bright_target_check=True. + pwr_upper_percentile : float, default=99.5 + Upper tail percentile for power spread computation (recommended: 98.0 or 99.5). + 99.5 = top 0.5% of power samples, 98.0 = top 2% of power samples. + Only used when bright_target_check=True. + sig_ev_margin_upper_db : float, default=3.0 + Conservative safety margin in dB applied when PCR is low (weak RFI). + Upper bound of the adaptive margin range. Only used when threshold_method='max_ev'. + sig_ev_margin_lower_db : float, default=0.0 + Aggressive safety margin in dB applied when PCR is high (strong RFI). + Lower bound of the adaptive margin range. Only used when threshold_method='max_ev'. + pcr_range : list of 2 floats, default=[0.1, 0.4] + (lower, upper) bounds of the PCR interpolation range. + CPIs with PCR <= lower use sig_ev_margin_upper_db (conservative). + CPIs with PCR >= upper use sig_ev_margin_lower_db (aggressive). + Only used when threshold_method='max_ev'. + + DUAL-CHECK LOGIC (when bright_target_check=True): + -------------------------------------------------- + TB is skipped only if BOTH checks pass (AND logic): + IF (cond# std <= bt_condition_num_thresh_db) AND (upper tail spread <= bt_pwr_spread_thresh_db): + SKIP TB (Clean or Bright Target) + ELSE: + PROCEED with RFI detection + + When bright_target_check=False: + All TBs proceed to RFI detection (no bright target filtering) Returns -------- @@ -117,10 +191,11 @@ def rfi_detect( "Total number of pulses must be greater or equal to number of pulses per single CPI." ) - # Need to validate sample covariance rank ( - eig_val_sort_array, + eig_val_sort_array, eig_vec_sort_array, + diag_power_array, + diag_valid_array, tb_is_valid, ) = compute_evd_tb( raw_data, @@ -131,34 +206,270 @@ def rfi_detect( min_ev_valid_idx=min_ev_valid_idx, rx_dynamic_range_db=rx_dynamic_range_db, ) + # If any CPI within a threshold block is determined to be invalid # Then skip threshold computation for this block by setting rfi_cpi_flag_array # to all zeros if not tb_is_valid: num_cpi = eig_val_sort_array.shape[0] rfi_cpi_flag_array = np.zeros((num_cpi, cpi_len), dtype=np.bool_) + fig_merit_detect_tb = 0 return ( rfi_cpi_flag_array, eig_vec_sort_array, + diag_power_array, + diag_valid_array, + fig_merit_detect_tb, ) - # Estimate a single threshold for all CPIs - detect_threshold = threshold_estimate_evd( - eig_val_sort_array, - num_max_trim, - num_min_trim, - max_num_rfi_ev, - min_ev_valid_idx, - threshold_params, - ) + # DUAL-CHECK: Condition number variability + Power spread + # Both checks must pass to skip TB (conservative approach) + if bright_target_check: + # Check 1: Condition number variability + condition_number_std_db = compute_tb_condition_number_std( + eig_val_sort_array, + min_ev_valid_idx, + ) + cond_num_passes = (np.isfinite(condition_number_std_db) and + condition_number_std_db <= bt_condition_num_thresh_db) + + # Check 2: Power spread (upper tail spread) + upper_tail_spread_db = compute_tb_upper_tail( + diag_power_array, + diag_valid_array, + pwr_ref_percentile=pwr_ref_percentile, + pwr_upper_percentile=pwr_upper_percentile, + ) + power_stat_passes = (np.isfinite(upper_tail_spread_db) and + upper_tail_spread_db <= bt_pwr_spread_thresh_db) + + # DECISION: SKIP only if BOTH checks pass + if cond_num_passes and power_stat_passes: + # Clean or Bright Target: stable structure + stationary power + num_cpi = eig_val_sort_array.shape[0] + rfi_cpi_flag_array = np.zeros((num_cpi, cpi_len), dtype=np.bool_) + fig_merit_detect_tb = 0 + + return ( + rfi_cpi_flag_array, + eig_vec_sort_array, + diag_power_array, + diag_valid_array, + fig_merit_detect_tb, + ) + + # TB has potential RFI (at least one check failed): proceed with detection + if threshold_method == 'ev_slope': + # Estimate a single threshold for all CPIs + detect_threshold, fig_merit_detect_tb = threshold_estimate_evd( + eig_val_sort_array, + num_max_trim, + num_min_trim, + max_num_rfi_ev, + min_ev_valid_idx, + threshold_params, + ) + elif threshold_method == 'max_ev': + detect_threshold, fig_merit_detect_tb = threshold_estimate_max_ev( + eig_val_sort_array, + max_deg_freedom=max_deg_freedom, + num_rfi_buffer=num_rfi_buffer, + num_max_trim=num_max_trim, + num_min_trim=num_min_trim, + max_num_rfi_ev=max_num_rfi_ev, + min_ev_valid_idx=min_ev_valid_idx, + rfi_fig_merit_thresh=rfi_fig_merit_thresh, + sig_ev_margin_upper_db=sig_ev_margin_upper_db, + sig_ev_margin_lower_db=sig_ev_margin_lower_db, + pcr_range=pcr_range, + ) + + # FoM below threshold: eigenvalue structure does not indicate RFI presence. + # Return early with no detections — distinct from invalid TB and bright target skip. + if not np.isfinite(detect_threshold).all(): + num_cpi = eig_val_sort_array.shape[0] + rfi_cpi_flag_array = np.zeros((num_cpi, cpi_len), dtype=np.bool_) + return ( + rfi_cpi_flag_array, + eig_vec_sort_array, + diag_power_array, + diag_valid_array, + fig_merit_detect_tb, + ) + else: + raise ValueError(f"Unsupported threshold method: {threshold_method}") # Detect RFI Eigenvalues of each CPI based on input detection threshold rfi_cpi_flag_array = rfi_detect_evd_tb( - eig_val_sort_array, detect_threshold, max_deg_freedom + eig_val_sort_array, detect_threshold, max_deg_freedom, threshold_method, + ) + + return rfi_cpi_flag_array, eig_vec_sort_array, diag_power_array, diag_valid_array, fig_merit_detect_tb + +def threshold_estimate_max_ev( + eig_val_sort_array, + max_deg_freedom=4, + num_rfi_buffer=3, + num_max_trim=0, + num_min_trim=0, + max_num_rfi_ev=2, + min_ev_valid_idx=10, + rfi_fig_merit_thresh=1.0, + sig_ev_margin_upper_db=3.0, + sig_ev_margin_lower_db=0.0, + pcr_range=[0.1, 0.4], +): + """Estimate signal max eigenvalue threshold using max_ev detection algorithm. + + Algorithm: + 1. Compute FoM using EST method + 2. Check if FoM indicates potential RFI + 3. If RFI confirmed (condition number check already performed in rfi_detect), + compute signal_ev_max_estimate by: + - Computing adaptive margin per CPI based on Principal Component Ratio (PCR) + - PCR = ev[0] / sum(ev[1:]) measures dominance of first eigenvalue + - High PCR (strong RFI) -> reduced margin for aggressive detection + - Low PCR (weak RFI) -> full margin for conservative detection + - Using clean eigenvalue segment from (max_deg_freedom+1) to min_ev_valid_idx + - Performing linear regression across this segment for each CPI independently + - Extrapolating each CPI's fitted line to index 0 to estimate maximum signal eigenvalue + - Adding adaptive margin (varies per CPI based on PCR) + - Returning per-CPI thresholds (no averaging across CPIs) + + Parameters + ---------- + eig_val_sort_array : 2D array of float, [num_cpi x cpi_len] + Sorted eigenvalues in descending order + max_deg_freedom : int, default=4 + Maximum number of independent RFI emitters designed to be detected. + Eigenvalue indices [0, ..., max_deg_freedom-1] can potentially be RFI. + num_rfi_buffer : int, default=3 + Number of buffer eigenvalue indices to skip after last possible RFI EV + before starting clean segment interpolation. + start_idx = max_deg_freedom + num_rfi_buffer - 1. + num_max_trim : int, default=0 + Number of large-value outliers to trim for FoM computation + num_min_trim : int, default=0 + Number of small-value outliers to trim for FoM computation + max_num_rfi_ev : int, default=2 + Maximum number of dominant EVs to check for FoM computation + min_ev_valid_idx : int, default=10 + Eigenvalue index for noise floor estimation + rfi_fig_merit_thresh : float, default=1.0 + Figure of merit threshold. TB proceeds to threshold estimation only if + fig_merit >= rfi_fig_merit_thresh. Lower values allow more TBs to be checked. + sig_ev_margin_upper_db : float, default=3.0 + Conservative safety margin in dB applied when PCR is low (weak RFI). + Upper bound of the adaptive margin interpolation range. + sig_ev_margin_lower_db : float, default=0.0 + Aggressive safety margin in dB applied when PCR is high (strong RFI). + Lower bound of the adaptive margin interpolation range. + Higher PCR indicates stronger RFI dominance, allowing more aggressive detection. + pcr_range : list of 2 floats, default=[0.1, 0.4] + (lower, upper) bounds of the PCR interpolation range. + CPIs with PCR <= lower use sig_ev_margin_upper_db (conservative). + CPIs with PCR >= upper use sig_ev_margin_lower_db (aggressive). + + Returns + ------- + detect_threshold : float or ndarray + Per-CPI signal max eigenvalue estimates in dB if RFI detected, else np.inf. + Returns a 1D array of length num_cpi, where each element is the threshold + for the corresponding CPI. + fig_merit : float + Figure of merit value computed + """ + + # Step 1: Compute FoM using EST method + _, fig_merit = threshold_estimate_evd( + eig_val_sort_array, + num_max_trim=num_max_trim, + num_min_trim=num_min_trim, + max_num_rfi_ev=max_num_rfi_ev, + min_ev_valid_idx=min_ev_valid_idx, + threshold_params=ThresholdParams(x=[2.0, 20.0], y=[5.0, 2.0]), ) - return rfi_cpi_flag_array, eig_vec_sort_array + # Step 2: Check if FoM indicates potential RFI + if fig_merit < rfi_fig_merit_thresh: + # Low FoM, no RFI detected + return np.inf, fig_merit + + # Step 3: RFI confirmed (high FoM, and power spread already checked in rfi_detect) + # Compute signal_ev_max_estimate using clean eigenvalue segment + eps = np.finfo(np.float64).tiny + eig_val = np.maximum(np.real(eig_val_sort_array), eps) + eig_val_sort_db_array = 10.0 * np.log10(eig_val) + + # Define clean eigenvalue segment: from (max_deg_freedom + num_rfi_buffer - 1) to min_ev_valid_idx + # max_deg_freedom indices [0, ..., max_deg_freedom-1] can be RFI + # Add num_rfi_buffer indices as buffer, then start at max_deg_freedom + num_rfi_buffer - 1 + # This segment should be free from RFI contamination and usable for interpolation + start_idx = max_deg_freedom + num_rfi_buffer - 1 + end_idx = min_ev_valid_idx + + # Validate that clean segment has sufficient points + if start_idx >= end_idx: + raise ValueError( + f"Invalid clean eigenvalue segment: start_idx ({start_idx}) >= end_idx ({end_idx}). " + f"Require min_ev_valid_idx ({min_ev_valid_idx}) > max_deg_freedom + num_rfi_buffer - 1 " + f"({max_deg_freedom} + {num_rfi_buffer} - 1 = {start_idx})." + ) + + # Ensure we have at least 2 points for linear fit + if end_idx <= start_idx + 1: + # Fallback: segment too short for regression (< 2 points) + # Use simple slope estimate from last 2 eigenvalues before min_ev_valid_idx + # Returns a scalar threshold (averaged across all CPIs) + slopes = (eig_val_sort_db_array[:, min_ev_valid_idx - 1] - + eig_val_sort_db_array[:, min_ev_valid_idx]) + slope_avg = np.mean(slopes) + ev_min_avg = np.mean(eig_val_sort_db_array[:, min_ev_valid_idx]) + signal_ev_max_estimate = slope_avg * min_ev_valid_idx + ev_min_avg + sig_ev_margin_upper_db + else: + # Use linear regression across clean eigenvalue segment for each CPI + # Generate one threshold per CPI (no averaging) + num_cpi = eig_val_sort_db_array.shape[0] + signal_ev_max_estimate = np.zeros(num_cpi) + + for i in range(num_cpi): + # Compute adaptive margin based on Principal Component Ratio (PCR) + # PCR = ev[0] / sum(ev[1:]) measures dominance of first eigenvalue + # High PCR (strong RFI) -> reduce margin for more aggressive detection + # Low PCR (weak RFI or clutter) -> keep full margin for conservative detection + ev_dominant = eig_val[i, 0] + ev_sum_rest = np.sum(eig_val[i, 1:]) + + if ev_sum_rest > eps: + pcr = ev_dominant / ev_sum_rest + else: + pcr = 0.0 # Fallback if sum is zero + + # Adaptive margin interpolation: + # PCR < pcr_range[0]: use upper margin (conservative) + # PCR in pcr_range: interpolate between upper and lower margin + # PCR > pcr_range[1]: use lower margin (aggressive) + pcr_clamped = np.clip(pcr, pcr_range[0], pcr_range[1]) + adaptive_margin_db = np.interp( + pcr_clamped, + pcr_range, # PCR range + [sig_ev_margin_upper_db, sig_ev_margin_lower_db] # Margin range + ) + + # Extract clean eigenvalue segment for this CPI + ev_segment = eig_val_sort_db_array[i, start_idx:end_idx] + indices = np.arange(start_idx, end_idx) + + # Linear regression: ev_db = slope * index + intercept + # Using polyfit (degree=1 for linear) + coeffs = np.polyfit(indices, ev_segment, deg=1) + intercept = coeffs[1] # dB at index 0 + + # Extrapolate to index 0 for this CPI with adaptive margin + signal_ev_max_estimate[i] = intercept + adaptive_margin_db + + return signal_ev_max_estimate, fig_merit def threshold_estimate_evd( eig_val_sort_array, @@ -169,12 +480,12 @@ def threshold_estimate_evd( threshold_params: ThresholdParams = ThresholdParams(), ): """Perform data-centric thresholding algorithm: "Slow-Time Eigenvalue Slope - Thresholding "(ST-EST)"[1]_ based on the assumption that first difference of + Thresholding "(ST-EST)"[1]_ based on the assumption that first difference of minimum Eigenvalue across slow time is an estimate of signal power variation as if there is no RFI" Algorithm Overview for applying ST-EST on one raw data block: - 1. Remove a specified number of outliers in maximum and minimum Eigenvalues + 1. Remove a specified number of outliers in maximum and minimum Eigenvalues 2. Compute slow-time standard deviation (STD) of maximum Eigenvalue slope. 3. Compute slow-time standard deviation (STD) of minimum Eigenvalue slope. 4. Compute STD ratio of maximum and minimum Eigenvalue slopes (SRMMES). @@ -196,9 +507,9 @@ def threshold_estimate_evd( num_min_trim: int, default = 0 Number of small value outlliers to be trimmed in slow-time minimum Eigenvalues max_num_rfi_ev: int, default = 2 - A detection error (miss) happens when a maximum power RFI emitter contaminates - multiple consecutive CPIs, resulting in a flat maximum Eigenvalue slope in slow - time. Hence the standard (STD) deviation of multiple dominant EVs across slow time + A detection error (miss) happens when a maximum power RFI emitter contaminates + multiple consecutive CPIs, resulting in a flat maximum Eigenvalue slope in slow + time. Hence the standard (STD) deviation of multiple dominant EVs across slow time defined by this parameter are compared. The one with the maximum STD is used for RFI Eigenvalue first difference computation. min_ev_valid_idx: int @@ -244,8 +555,8 @@ def threshold_estimate_evd( ev_min_db = np.delete(ev_min_db, ev_min_trim_idx) if num_max_trim > 0: - ev_max_trim_idx = np.argsort(ev_min_db)[-num_max_trim:] - ev_min_db = np.delete(ev_min_db, ev_max_trim_idx) + ev_max_trim_idx = np.argsort(ev_max_db)[-num_max_trim:] + ev_max_db = np.delete(ev_max_db, ev_max_trim_idx) # Compute STD of the slope of max and min Eigenvalues ev_slope_max = np.diff(ev_max_db) @@ -267,7 +578,7 @@ def threshold_estimate_evd( detect_threshold = ev_slope_min_mean + num_sigma * ev_slope_min_std - return detect_threshold + return detect_threshold, std_ratio_ev_slope def rfi_detect_evd( @@ -275,9 +586,9 @@ def rfi_detect_evd( detect_threshold, max_deg_freedom=8, ): - """Perform RFI detection of Eigenvalues within a CPI based on input detection - threshold in dB/Eigenvalue index. The threshold is set to be a negative value. - If the magnitude of an Eigenvalue slope exceeds this threshold, then it is identified + """Perform RFI detection of Eigenvalues within a CPI based on input detection + threshold in dB/Eigenvalue index. The threshold is set to be a negative value. + If the magnitude of an Eigenvalue slope exceeds this threshold, then it is identified as RFI. Parameters @@ -287,7 +598,7 @@ def rfi_detect_evd( detect_threshold: float A positive RFI detection threshold used by EVD detection algorithm to identify RFI Eigenvalue slope valules. - max_deg_freedom: int, default = 12 + max_deg_freedom: int, default=8 Max number of independent RFI emitters designed to be detected and mitigated. This number should be less than cpi_len. @@ -309,38 +620,63 @@ def rfi_detect_evd( def rfi_detect_evd_tb( eig_val_sort_array, detect_threshold, - max_deg_freedom=8, + max_deg_freedom=4, + threshold_method='max_ev', ): - """Wrapper function which performs RFI detection of data within a Threshold Block (TB) - one CPI at a time base don input detection threshold in dB/Eigenvalue index. + """Wrapper function which performs RFI detection of data within a Threshold Block (TB) + one CPI at a time based on input detection threshold. + + Supports two detection methods: + - 'ev_slope': Eigenvalue Slope Thresholding (slope-based detection) + - 'max_ev': EV Maximum Estimation (absolute eigenvalue threshold from noise extrapolation) Parameters ---------- eig_val_sort_array: 2D array of float, [num_cpi x cpi_len] Sorted Eigenvalues in descending order of all CPIs in raw data - detect_threshold: float - A positive RFI detection threshold used by EVD detection algorithm to identify - RFI Eigenvalues in a CPI in dB/Eigenvalue index. All CPIs of the input data - share a common detection threshold. - max_deg_freedom: int, default = 12 + detect_threshold: float or ndarray + RFI detection threshold. For 'ev_slope': positive slope threshold in dB/index. + For 'max_ev': absolute eigenvalue threshold in dB. + Can be a scalar (same threshold for all CPIs) or a 1D array of length num_cpi + (one threshold per CPI). + max_deg_freedom: int, default = 4 Max number of independent RFI emitters designed to be detected and mitigated. This number should be less than cpi_len. + threshold_method: str, default='max_ev' + Detection method: 'ev_slope' or 'max_ev' Returns ------- rfi_cpi_flag_array: 2D array of bool, [num_cpi x cpi_len] RFI flag array that marks each Eigenvalue index in a CPI as either RFI or signal. - 1 = RFI Eigenvalue index; 0 = Signal Eigenvalue index + True = RFI Eigenvalue index; False = Signal Eigenvalue index """ # Number of pulses, CPI length, and number of range blocks in a CPI num_cpi, cpi_len = eig_val_sort_array.shape - # Ensure detection threshold is a positive value - if (not np.isfinite(detect_threshold)) or (detect_threshold <= 0): - # Add a CPI-Level rank check: detection threshold must be a positive value - warnings.warn("Warning: Non-positive detection threshold. Skipping TB detection.") - return np.zeros((num_cpi, cpi_len), dtype=np.bool_) + # Convert threshold to array for uniform handling + detect_threshold_arr = np.atleast_1d(detect_threshold) + + # Validate threshold array shape + if detect_threshold_arr.size == 1: + # Scalar threshold - broadcast to all CPIs + detect_threshold_arr = np.full(num_cpi, detect_threshold_arr[0]) + elif detect_threshold_arr.size != num_cpi: + raise ValueError(f"Threshold array size ({detect_threshold_arr.size}) must match num_cpi ({num_cpi})") + + # Validate threshold values + if threshold_method == 'ev_slope': + # Ensure detection threshold is a positive value + if np.any(~np.isfinite(detect_threshold_arr)) or np.any(detect_threshold_arr <= 0): + warnings.warn("Warning: Non-positive detection threshold. Skipping TB detection.") + return np.zeros((num_cpi, cpi_len), dtype=np.bool_) + elif threshold_method == 'max_ev': + if np.any(~np.isfinite(detect_threshold_arr)): + warnings.warn(f"Warning: Non-finite {threshold_method.upper()} detection threshold. Skipping TB detection.") + return np.zeros((num_cpi, cpi_len), dtype=np.bool_) + else: + raise ValueError(f"Unsupported threshold_method: {threshold_method}") # Maximum number of degrees of freedom must be less than cpi_len if max_deg_freedom >= cpi_len: @@ -348,21 +684,103 @@ def rfi_detect_evd_tb( "Max number of deg. of freedom must be less than number of pulses in a CPI." ) - # RFI flag for each eigenvalue index in all CPIs: RFI=1, signal=0:w + # RFI flag for each eigenvalue index in all CPIs: RFI=1, signal=0 rfi_cpi_flag_array = np.ones((num_cpi, cpi_len), dtype=np.bool_) - # Compute Eigenvalue Slope or first difference of all CPIs + # Compute Eigenvalue in dB for all CPIs eig_val_sort_db_array = 10 * np.log10(np.abs(eig_val_sort_array)) - eig_val_db_slope_array = np.diff(eig_val_sort_db_array, axis=1) - for idx_cpi in range(num_cpi): - eig_val_db = eig_val_sort_db_array[idx_cpi] - eig_val_db_slope = eig_val_db_slope_array[idx_cpi] + if threshold_method == 'ev_slope': + eig_val_db_slope_array = np.diff(eig_val_sort_db_array, axis=1) + for idx_cpi in range(num_cpi): + eig_val_db_slope = eig_val_db_slope_array[idx_cpi] + + # Use per-CPI threshold + sig_ev_idx_start = rfi_detect_evd(eig_val_db_slope, detect_threshold_arr[idx_cpi], max_deg_freedom) + + # Sets signal eigenvalue indices to False + rfi_cpi_flag_array[idx_cpi, sig_ev_idx_start:] = False + else: + # 'max_ev' uses absolute eigenvalue comparison + for idx_cpi in range(num_cpi): + eig_val_db_valid = eig_val_sort_db_array[idx_cpi, :max_deg_freedom] + + # Use per-CPI threshold + rfi_ev_idx = np.where(eig_val_db_valid > detect_threshold_arr[idx_cpi])[0] - # Determine starting index of signal EV - sig_ev_idx_start = rfi_detect_evd(eig_val_db_slope, detect_threshold, max_deg_freedom) + if rfi_ev_idx.size: + sig_ev_idx_start = rfi_ev_idx[-1] + 1 - # Sets signal eigenvalue indices to zero - rfi_cpi_flag_array[idx_cpi, sig_ev_idx_start:] = 0 + # Sets signal eigenvalue indices to zero + rfi_cpi_flag_array[idx_cpi, sig_ev_idx_start:] = False + else: + rfi_cpi_flag_array[idx_cpi, :] = False return rfi_cpi_flag_array + +def compute_tb_upper_tail( + diag_power_array, + diag_valid_array, + pwr_ref_percentile=50, + pwr_upper_percentile=99.5, + eps=1e-12, +): + pwr = np.asarray(diag_power_array)[diag_valid_array] + + pwr = pwr[np.isfinite(pwr)] + + if pwr.size == 0: + return np.nan + + # Convert to dB + pwr_db = 10 * np.log10(np.maximum(pwr, eps)) + + # Remove non-finite values + pwr_db = pwr_db[np.isfinite(pwr_db)] + + # Compute upper-tail spread + baseline = np.percentile(pwr_db, pwr_ref_percentile) + upper = np.percentile(pwr_db, pwr_upper_percentile) + + upper_tail_db = upper - baseline + + return upper_tail_db + +def compute_tb_condition_number_std( + eig_val_sort_array, + min_ev_valid_idx, + eps=1e-12, +): + """Compute standard deviation of condition numbers across CPIs in a TB. + + Parameters + ---------- + eig_val_sort_array : 2D array [num_cpi x cpi_len] + Sorted eigenvalues in descending order + min_ev_valid_idx : int + Index of minimum valid eigenvalue + eps : float + Small value to avoid log(0) + + Returns + ------- + condition_number_std_db : float + Standard deviation of condition numbers in dB + """ + num_cpi = eig_val_sort_array.shape[0] + condition_numbers_db = np.zeros(num_cpi) + + for i in range(num_cpi): + ev_max = np.maximum(eig_val_sort_array[i, 0], eps) + ev_min = np.maximum(eig_val_sort_array[i, min_ev_valid_idx], eps) + condition_numbers_db[i] = 10 * np.log10(ev_max) - 10 * np.log10(ev_min) + + # Remove non-finite values + condition_numbers_db = condition_numbers_db[np.isfinite(condition_numbers_db)] + + if condition_numbers_db.size == 0: + return np.nan + + condition_number_std_db = np.std(condition_numbers_db) + + return condition_number_std_db diff --git a/python/packages/isce3/signal/rfi_process_evd.py b/python/packages/isce3/signal/rfi_process_evd.py index db09e179f..d23574507 100644 --- a/python/packages/isce3/signal/rfi_process_evd.py +++ b/python/packages/isce3/signal/rfi_process_evd.py @@ -1,6 +1,8 @@ """ Perform RFI detection and mitigation of input raw data using Slow-Time Eigenvalue Decomposition (ST-EVD). + +DUAL-CHECK VERSION: Uses condition number check AND power spread check """ import numpy as np from isce3.signal.compute_evd_cpi import slice_gen @@ -13,19 +15,30 @@ def run_slow_time_evd( cpi_len, max_deg_freedom, *, + num_rfi_buffer=3, num_max_trim=0, num_min_trim=0, max_num_rfi_ev=2, num_samples_rng_blk=250, use_entire_pulse=False, threshold_params: ThresholdParams = ThresholdParams(), - num_cpi_tb=20, - off_diag_overlap_ratio=0.25, - diag_valid_ratio=0.20, + num_cpi_per_threshold_block=12, + off_diag_overlap_ratio=0.20, + diag_valid_ratio=0.15, mitigate_enable=False, min_rank_frac=0.70, rx_dynamic_range_db=50.0, swaths=None, + threshold_method='max_ev', + rfi_fig_merit_thresh=1.0, + bright_target_check=True, + bt_condition_num_thresh_db=2.0, + bt_pwr_spread_thresh_db=6.0, + pwr_ref_percentile=50, + pwr_upper_percentile=99.5, + sig_ev_margin_upper_db=3.0, + sig_ev_margin_lower_db=0.0, + pcr_range=[0.1, 0.4], raw_data_mitigated=None, ): @@ -34,6 +47,8 @@ def run_slow_time_evd( Each TB is consisted of M Coherent Processing Intervals (CPI) and N range samples Each CPI is consisted of K slow-time pulses. 2. Derive slow-time RFI detection threshold for each TB. + Note: Bright target check (dual condition number + power spread check) is performed + inside rfi_detect() before threshold computation for all methods. 3. Mitigate RFI of all CPIs above the detection threshold if mitigation is enabled. Parameters @@ -45,37 +60,39 @@ def run_slow_time_evd( max_deg_freedom: int Max number of independent RFI emitters designed to be detected and mitigated. This number should be less than cpi_len to avoid unintended removal of signal data. + num_rfi_buffer: int, default=3 + Number of buffer eigenvalue indices to skip after last possible RFI EV before + starting clean segment interpolation for 'max_ev' method. The clean segment + starts at index (max_deg_freedom + num_rfi_buffer - 1). num_max_trim: int, default=0 Number of large value outliers to be trimmed in slow-time minimum Eigenvalues. num_min_trim: int, default=0 Number of small value outliers to be trimmed in slow-time minimum Eigenvalues max_num_rfi_ev: int, default=2 - A detection error (miss) happens when a maximum power RFI emitter contaminates - multiple consecutive CPIs, resulting in a flat maximum Eigenvalue slope in slow - time. Hence the standard (STD) deviation of multiple dominant EVs across slow time + A detection error (miss) happens when a maximum power RFI emitter contaminates + multiple consecutive CPIs, resulting in a flat maximum Eigenvalue slope in slow + time. Hence the standard (STD) deviation of multiple dominant EVs across slow time defined by this parameter are compared. The one with the maximum STD is used for RFI Eigenvalue first difference computation. - num_samples_rng_blk: int + num_samples_rng_blk: int, default=250 Number of range samples per range block when data blockin is applied in range direction - for sample covariance matrix estimation. It is recommended that this parameter is at - least 5 x cpi_len to avoid discrepancy from true sample covaraince matrix. In addition, - in order to avoid a run-time error for ST-EVD, this parameter needs to be - at least 2 x cpi. - default = 256 - use_entire_pulse: bool - Ignore any value passed for num_samples_rng_blk and instead Use all samples - in the slow-time pulses for detection if this is True - default = False + for sample covariance matrix estimation. It is recommended that this parameter is at + least 5 x cpi_len to avoid discrepancy from true sample covaraince matrix. In addition, + in order to avoid a run-time error for ST-EVD, this parameter needs to be + at least 2 x cpi_len. + use_entire_pulse: bool, default=False + Ignore any value passed for num_samples_rng_blk and instead use all samples + in the slow-time pulses for detection if this is True. threshold_params: ThresholdParams object, default=ThresholdParams() RFI detection threshold interpolation parameters. The x field defines STD ratio between maximum and minimum Eigenvalue slopes (MMES) of the slow-time threshold interval. The y field defines the number of sigma (STD) from the mean of MMES. - num_cpi_tb: int, default=20 + num_cpi_per_threshold_block: int, default=12 Number of slow-time CPIs in a TB - off_diag_overlap_ratio : float, optional, default = 0.25 + off_diag_overlap_ratio : float, optional, default=0.20 Minimum overlap ratio used by gap exclusion covariance estimation - diag_valid_ratio : float, optional, default = 0.20 + diag_valid_ratio : float, optional, default=0.15 Minimum fraction of valid samples required to compute a diagonal term in the sample covariance matrix entry R_ii. mitigate_enable: bool, default=False @@ -94,6 +111,59 @@ def run_slow_time_evd( the [start, stop) indices of the sub-swath. It's recommended to supply this for modes with dithered PRI, where it will be used to normalize the sample covariance matrix. + threshold_method : str, default='max_ev' + RFI detection method: 'ev_slope' (Eigenvalue Slope Thresholding) or + 'max_ev' (EV maximum estimation). + DUAL-CHECK VERSION: Uses condition number check AND power spread check. + rfi_fig_merit_thresh : float, default=1.0 + Figure of merit threshold for 'max_ev' method. Default of 1.0 ensures + aggressive RFI detection - checks as many TBs as possible. + Only used when threshold_method='max_ev'. + bright_target_check : bool, default=True + Enable dual bright target rejection check using BOTH condition number variability + AND power stationarity. When True, TB is skipped only if BOTH checks pass: + - Condition number std <= bt_condition_num_thresh_db (stable eigenvalue structure) + - Upper tail spread <= bt_pwr_spread_thresh_db (low power spread) + When False, bright target checks are disabled (all TBs are processed). + bt_condition_num_thresh_db : float, default=2.0 + Condition number std threshold in dB. TBs with std(cond#) <= this value + pass the condition# check (stable eigenvalue structure). + Only used when bright_target_check=True. + bt_pwr_spread_thresh_db : float, default=6.0 + Power spread threshold in dB. TBs with upper tail spread <= this value + pass the power check (low power spread across slow time). + Upper tail spread = pwr_upper_percentile - pwr_ref_percentile of diagonal power. + Only used when bright_target_check=True. + pwr_ref_percentile : float, default=50 + Reference (baseline) percentile for power spread computation. + Represents the median power level of the diagonal covariance entries. + Only used when bright_target_check=True. + pwr_upper_percentile : float, default=99.5 + Upper percentile for power spread computation. + Recommended values: 98.0 (top 2%) or 99.5 (top 0.5%). + Only used when bright_target_check=True. + sig_ev_margin_upper_db : float, default=3.0 + Conservative safety margin in dB applied when PCR is low (weak RFI). + Upper bound of the adaptive margin range. Only used when threshold_method='max_ev'. + sig_ev_margin_lower_db : float, default=0.0 + Aggressive safety margin in dB applied when PCR is high (strong RFI). + Lower bound of the adaptive margin range. Only used when threshold_method='max_ev'. + pcr_range : list of 2 floats, default=[0.1, 0.4] + (lower, upper) bounds of the PCR interpolation range. + CPIs with PCR <= lower use sig_ev_margin_upper_db (conservative). + CPIs with PCR >= upper use sig_ev_margin_lower_db (aggressive). + Only used when threshold_method='max_ev'. + + DUAL-CHECK LOGIC (when bright_target_check=True): + -------------------------------------------------- + TB is skipped only if BOTH checks pass (AND logic): + IF (cond# std <= bt_condition_num_thresh_db) AND (upper tail spread <= bt_pwr_spread_thresh_db): + SKIP TB (Clean or Bright Target) + ELSE: + PROCEED with RFI detection on TB + + When bright_target_check=False: + All TBs proceed to RFI detection (no bright target filtering) raw_data_mitigated: array-like complex [num_pulses x num_rng_samples] or None, optional output array in which the mitigated data values is placed. It must be an array-like object supporting `multidimensional array access @@ -125,19 +195,38 @@ def run_slow_time_evd( if use_entire_pulse: num_samples_rng_blk = num_rng_samples - num_rng_blks = num_rng_samples // num_samples_rng_blk - # If the number of pulses is not an integer multiple of TB size, following # operations will ensue. If the number of remaining pulses is greater than # the CPI length, additional CPI(s) will be constructed, the very last TB will # include the additional CPI(s). The rest of the remaining pulses not enough to - # construct a full CPI will not be processed. If the number of remaining pulses - # is less than CPI length, then they will not be processed. In both cases, + # construct a full CPI will not be processed. If the number of remaining pulses + # is less than CPI length, then they will not be processed. In both cases, # at most cpi_len-1 number of pulses will be ignored. num_cpi = num_pulses // cpi_len num_pulses_proc = cpi_len * num_cpi - num_pulses_tb = cpi_len * num_cpi_tb + num_pulses_tb = cpi_len * num_cpi_per_threshold_block + num_tb = num_pulses_proc // num_pulses_tb + + # Figue out how many range slices are there + rng_slices = list( + slice_gen(num_rng_samples, num_samples_rng_blk, combine_rem=True) + ) + num_rng_blks = len(rng_slices) + + # RFI EV Count Map + rfi_ev_count_map = np.zeros( + (num_cpi, num_rng_blks), + dtype=np.int16, + ) + + # Boolean Detection Map + rfi_cpi_detection_map = np.zeros( + (num_cpi, num_rng_blks), + dtype=bool, + ) + + figure_merit_array = np.zeros((num_tb, num_rng_blks), dtype=np.float32) # Modify raw_data in-place if raw_data_mitigated is None: @@ -187,7 +276,7 @@ def run_slow_time_evd( ) # Run RFI Detection and Mitigation - for _, tb_slow_time in enumerate(slice_gen(num_pulses_proc, num_pulses_tb)): + for idx_tb, tb_slow_time in enumerate(slice_gen(num_pulses_proc, num_pulses_tb)): # Get valid data mask for all rows in current block. if swaths is not None: swaths_tb = swaths[:, tb_slow_time, :] @@ -196,20 +285,22 @@ def run_slow_time_evd( for start, end in swaths_tb[:, i, :]: mask_valid[i, start:end] = True - for _, tb_fast_time in enumerate( - slice_gen(num_rng_samples, num_samples_rng_blk, combine_rem=True) - ): + for idx_rng, tb_fast_time in enumerate(rng_slices): raw_tb_blk = raw_data[tb_slow_time, tb_fast_time] mask_valid_tb = None if swaths is None else mask_valid[:, tb_fast_time] ( - rfi_cpi_flag_tb, - evec_sort_tb, + rfi_cpi_flag_tb, + evec_sort_tb, + _, # diag_power_array (unused, power check now in rfi_detect) + _, # diag_valid_array (unused, power check now in rfi_detect) + figure_merit_tb, ) = rfi_detect( raw_tb_blk, cpi_len, max_deg_freedom, min_ev_valid_idx, + num_rfi_buffer=num_rfi_buffer, num_max_trim=num_max_trim, num_min_trim=num_min_trim, max_num_rfi_ev=max_num_rfi_ev, @@ -217,22 +308,46 @@ def run_slow_time_evd( diag_valid_ratio=diag_valid_ratio, rx_dynamic_range_db=rx_dynamic_range_db, mask_valid=mask_valid_tb, + threshold_method=threshold_method, threshold_params=threshold_params, + rfi_fig_merit_thresh=rfi_fig_merit_thresh, + bright_target_check=bright_target_check, + bt_condition_num_thresh_db=bt_condition_num_thresh_db, + bt_pwr_spread_thresh_db=bt_pwr_spread_thresh_db, + pwr_ref_percentile=pwr_ref_percentile, + pwr_upper_percentile=pwr_upper_percentile, + sig_ev_margin_upper_db=sig_ev_margin_upper_db, + sig_ev_margin_lower_db=sig_ev_margin_lower_db, + pcr_range=pcr_range, ) + # Global CPI indices for this threshold block + cpi_start = idx_tb * num_cpi_per_threshold_block + cpi_end = cpi_start + rfi_cpi_flag_tb.shape[0] + + # Power stationarity check is now handled inside rfi_detect() + # Detection returns no RFI (all zeros) for bright target TBs automatically + has_rfi = np.any(rfi_cpi_flag_tb) + # Compute number of CPIs detected with RFI presence - num_rfi_ev_cpi = np.sum(rfi_cpi_flag_tb, axis=1) + num_rfi_ev_cpi = np.sum(rfi_cpi_flag_tb, axis=1).astype(np.int16) rfi_cpi_count = np.sum(num_rfi_ev_cpi != 0) rfi_cpi_count_sum += rfi_cpi_count + figure_merit_array[idx_tb, idx_rng] = figure_merit_tb + rfi_ev_count_map[cpi_start:cpi_end, idx_rng] = num_rfi_ev_cpi + rfi_cpi_detection_map[cpi_start:cpi_end, idx_rng] = num_rfi_ev_cpi > 0 + # Run Mitigation: - if mitigate_enable: + if mitigate_enable and has_rfi: rfi_mitigate_tb( raw_tb_blk, evec_sort_tb, rfi_cpi_flag_tb, raw_data_mitigated[tb_slow_time, tb_fast_time], ) + else: + raw_data_mitigated[tb_slow_time, tb_fast_time] = raw_tb_blk # Percentage of RFI Eigenvalues based on slow-time min EV slope detection rfi_likelihood = rfi_cpi_count_sum / (num_cpi * num_rng_blks) diff --git a/python/packages/nisar/workflows/focus.py b/python/packages/nisar/workflows/focus.py index 4ec69e181..079ab98de 100644 --- a/python/packages/nisar/workflows/focus.py +++ b/python/packages/nisar/workflows/focus.py @@ -1163,19 +1163,30 @@ def process_rfi(cfg: Struct, raw_data: np.ndarray, raw_data, opt_evd.cpi_length, opt_evd.max_emitters, + num_rfi_buffer=opt_evd.num_rfi_buffer, num_max_trim=opt_evd.num_max_trim, num_min_trim=opt_evd.num_min_trim, max_num_rfi_ev=opt_evd.max_num_rfi_ev, num_samples_rng_blk=opt.num_samples_rng_blk, use_entire_pulse=opt.use_entire_pulse, threshold_params=threshold_params, - num_cpi_tb=opt_evd.num_cpi_per_threshold_block, + num_cpi_per_threshold_block=opt_evd.num_cpi_per_threshold_block, off_diag_overlap_ratio=opt_evd.off_diag_overlap_ratio, diag_valid_ratio=opt_evd.diag_valid_ratio, mitigate_enable=opt.mitigation_enabled, min_rank_frac=opt_evd.min_rank_frac, rx_dynamic_range_db=opt_evd.rx_dynamic_range_db, swaths=swaths, + threshold_method=opt_evd.threshold_method, + rfi_fig_merit_thresh=opt_evd.rfi_fig_merit_thresh, + bright_target_check=opt_evd.bright_target_check, + bt_condition_num_thresh_db=opt_evd.bt_condition_num_thresh_db, + bt_pwr_spread_thresh_db=opt_evd.bt_pwr_spread_thresh_db, + pwr_ref_percentile=opt_evd.pwr_ref_percentile, + pwr_upper_percentile=opt_evd.pwr_upper_percentile, + sig_ev_margin_upper_db=opt_evd.sig_ev_margin_upper_db, + sig_ev_margin_lower_db=opt_evd.sig_ev_margin_lower_db, + pcr_range=opt_evd.pcr_range, raw_data_mitigated=raw_data_mitigated) elif opt.mitigation_algorithm == "FDNF": opt_fnf = opt.freq_notch_filter diff --git a/share/nisar/defaults/focus.yaml b/share/nisar/defaults/focus.yaml index 8bbac52dc..37e9d925b 100644 --- a/share/nisar/defaults/focus.yaml +++ b/share/nisar/defaults/focus.yaml @@ -249,23 +249,65 @@ runconfig: # maximum STD is used for RFI Eigenvalue first difference # computation. max_num_rfi_ev: 2 + # Number of buffer CPIs added on each side of a detected RFI + # CPI to account for uncertainty in RFI onset and offset. + num_rfi_buffer: 3 # Number of slow-time CPIs grouped together for determining and # applying a common EV slope threshold. - num_cpi_per_threshold_block: 15 - # minimum fraction of samples required for cross-variance terms - # of sample covariance matrix - off_diag_overlap_ratio: 0.25 - # minimum fraction of samples required for covariance terms - # of sample covariance matrix - diag_valid_ratio: 0.20 + num_cpi_per_threshold_block: 12 + # Minimum fraction of samples required for cross-variance terms + # of sample covariance matrix. + off_diag_overlap_ratio: 0.20 + # Minimum fraction of samples required for diagonal covariance + # terms of sample covariance matrix. + diag_valid_ratio: 0.15 # This ratio will be used to determine the minimum number of valid Eigenvalues # required for a CPI. min_ev_valid_idx = int(np.round(min_rank_frac * cpi_len)) - 1 min_rank_frac: 0.70 - # radar platform receiver dynamic range, e.g. 50 dB. This is applied as a threshold - # to determine if the Eigenvalue under test is meaningfully signficant. + # Radar platform receiver dynamic range in dB. Applied as a threshold + # to determine if the Eigenvalue under test is meaningfully significant. rx_dynamic_range_db: 50.0 + # RFI detection threshold method. Supported values: + # - 'max_ev': Maximum Eigenvalue Estimation (default) + # - 'ev_slope': Eigenvalue slope (ST-EST) + threshold_method: max_ev + # Figure of merit threshold for RFI detection. The figure of merit is + # the STD ratio of the maximum vs. minimum Eigenvalue slopes across + # the threshold block. Values above this threshold trigger RFI detection. + rfi_fig_merit_thresh: 1.0 + # Enable bright target (e.g. strong point scatterer) rejection check. + # If enabled, threshold blocks where BOTH the condition number check + # AND the power spread check pass are skipped (not flagged as RFI). + bright_target_check: True + # Threshold (dB) on the standard deviation of the condition number + # (ratio of max to min diagonal power) across CPIs within a threshold + # block. A low value indicates a spatially stationary scene (possible + # bright target). + bt_condition_num_thresh_db: 2.0 + # Threshold (dB) on the power spread within a threshold block, + # defined as the difference between the upper and reference + # percentiles of the diagonal covariance power. A low value + # indicates a low-dynamic-range scene (possible bright target). + bt_pwr_spread_thresh_db: 6.0 + # Percentile used as the reference (baseline) power level when + # computing the power spread for bright target rejection. + pwr_ref_percentile: 50 + # Upper percentile used when computing the power spread for + # bright target rejection. + pwr_upper_percentile: 99.5 + # Upper bound of the adaptive detection margin (dB). Applied when + # the Principal Component Ratio (PCR) is below pcr_range[0]. + sig_ev_margin_upper_db: 3.0 + # Lower bound of the adaptive detection margin (dB). Applied when + # the Principal Component Ratio (PCR) is above pcr_range[1]. + sig_ev_margin_lower_db: 0.0 + # PCR range [lower, upper] used for adaptive margin interpolation. + # PCR = ev[0] / sum(ev[1:]) where ev are the sorted eigenvalues. + # The detection margin linearly interpolates between + # sig_ev_margin_upper_db and sig_ev_margin_lower_db over this range. + pcr_range: [0.1, 0.4] # Parameters used to adaptively determine RFI eigenvalue - # difference thresholds + # difference thresholds (used when threshold_method = 'max_ev') threshold_hyperparameters: # The computed sigma ratio of maximum and minimum Eigenvalue # first differences. It is a dimensionless figure of merit. diff --git a/share/nisar/schemas/focus.yaml b/share/nisar/schemas/focus.yaml index 76d5a64e6..2d947de1b 100644 --- a/share/nisar/schemas/focus.yaml +++ b/share/nisar/schemas/focus.yaml @@ -449,22 +449,64 @@ slow_time_evd_options: # compared. The one with the maximum STD is used for RFI Eigenvalue first # difference computation. max_num_rfi_ev: int(min=1, required=False) + # Number of buffer CPIs added on each side of a detected RFI CPI to account + # for uncertainty in RFI onset and offset. + num_rfi_buffer: int(min=0, required=False) # Number of slow-time CPIs grouped together for determining and applying a # common EV slope threshold. num_cpi_per_threshold_block: int(min=2, required=False) - # minimum fraction of samples required for cross-variance terms - # of sample covariance matrix - off_diag_overlap_ratio: num(min=0.0, required=False) - # minimum fraction of samples required for co-variance terms - # of sample covariance matrix - diag_valid_ratio: num(min=0.0, required=False) - #This ratio will be used to determine the minimum number of valid Eigenvalues - ##required for a CPI. min_ev_valid_idx = int(np.round(min_rank_frac * cpi_len)) - 1 + # Minimum fraction of samples required for cross-variance terms + # of sample covariance matrix. + off_diag_overlap_ratio: num(min=0.0, max=1.0, required=False) + # Minimum fraction of samples required for diagonal covariance terms + # of sample covariance matrix. + diag_valid_ratio: num(min=0.0, max=1.0, required=False) + # This ratio will be used to determine the minimum number of valid Eigenvalues + # required for a CPI. min_ev_valid_idx = int(np.round(min_rank_frac * cpi_len)) - 1 min_rank_frac: num(exclusiveMinimum=0.0, max=1.0, required=False) - # radar platform receiver dynamic range, e.g. 50 dB. This is applied as a threshold - # to determine if the Eigenvalue under test is meaningfully signficant. + # Radar platform receiver dynamic range in dB. Applied as a threshold to + # determine if the Eigenvalue under test is meaningfully significant. rx_dynamic_range_db: num(min=0.0, max=100.0, required=False) + # RFI detection threshold method: 'max_ev' (Maximum Eigenvalue Estimation, + # default) or 'ev_slope' (Eigenvalue slope / ST-EST). + threshold_method: enum('max_ev', 'ev_slope', required=False) + # Figure of merit threshold for RFI detection. The figure of merit is the + # STD ratio of the maximum vs. minimum Eigenvalue slopes across the threshold + # block. Values above this threshold trigger RFI detection. + rfi_fig_merit_thresh: num(min=0.0, required=False) + # Enable bright target (e.g. strong point scatterer) rejection check. + # Threshold blocks where BOTH the condition number check AND the power spread + # check pass are skipped (not flagged as RFI). + bright_target_check: bool(required=False) + # Threshold (dB) on the standard deviation of the condition number (ratio of + # max to min diagonal power) across CPIs within a threshold block. A low value + # indicates a spatially stationary scene (possible bright target). + bt_condition_num_thresh_db: num(min=0.0, required=False) + # Threshold (dB) on the power spread within a threshold block, defined as the + # difference between the upper and reference percentiles of the diagonal + # covariance power. A low value indicates a low-dynamic-range scene (possible + # bright target). + bt_pwr_spread_thresh_db: num(min=0.0, required=False) + # Percentile (in [0, 100]) used as the reference (baseline) power level when + # computing the power spread for bright target rejection. + pwr_ref_percentile: num(min=0.0, max=100.0, required=False) + # Upper percentile (in [0, 100]) used when computing the power spread for + # bright target rejection. Must be greater than pwr_ref_percentile. + pwr_upper_percentile: num(min=0.0, max=100.0, required=False) + # Upper bound of the adaptive detection margin in dB. Applied when the + # Principal Component Ratio (PCR) is below pcr_range[0]. + sig_ev_margin_upper_db: num(min=0.0, required=False) + # Lower bound of the adaptive detection margin in dB. Applied when the + # Principal Component Ratio (PCR) is above pcr_range[1]. + sig_ev_margin_lower_db: num(min=0.0, required=False) + # PCR range [lower, upper] for adaptive margin interpolation. + # PCR = ev[0] / sum(ev[1:]) where ev are the sorted eigenvalues. + # The detection margin interpolates between sig_ev_margin_upper_db and + # sig_ev_margin_lower_db over this range. Must be a list of 2 non-negative + # values with lower < upper. + pcr_range: list(num(min=0.0), min=2, max=2, required=False) # Parameters used to adaptively determine RFI eigenvalue difference thresholds + # (used when threshold_method = 'max_ev'). threshold_hyperparameters: include('rfi_threshold_options', required=False) freq_notch_filter_options: diff --git a/tests/python/packages/isce3/signal/rfi_process_evd.py b/tests/python/packages/isce3/signal/rfi_process_evd.py index d9c271839..334724aa8 100644 --- a/tests/python/packages/isce3/signal/rfi_process_evd.py +++ b/tests/python/packages/isce3/signal/rfi_process_evd.py @@ -170,8 +170,8 @@ def rfi_wb_gen( @pytest.mark.parametrize( - "cpi_len, num_samples_rng_blk, num_cpi_tb, max_deg_freedom, num_max_trim, num_min_trim, max_num_rfi_ev, use_entire_pulse, mitigate_enable, test_case", - [ + "cpi_len, num_samples_rng_blk, num_cpi_per_threshold_block, max_deg_freedom, num_max_trim, num_min_trim, max_num_rfi_ev, use_entire_pulse, mitigate_enable, test_case", + [ (32, 2148, 20, 10, 1, 1, 2, True, True, 'mitigate'), # No range blks (32, 256, 20, 10, 1, 1, 2, False, True, 'mitigate'), # 256 range samples / range blocks (32, 2148, 20, 10, 0, 0, 2, True, False, 'no-op'), # No-op: no rng blks @@ -182,7 +182,7 @@ def rfi_wb_gen( def test_slow_time_evd( cpi_len, num_samples_rng_blk, - num_cpi_tb, + num_cpi_per_threshold_block, max_deg_freedom, num_max_trim, num_min_trim, @@ -200,8 +200,8 @@ def test_slow_time_evd( num_max_trim=1, num_min_trim=1, max_num_rfi_ev=2, mitigate_enable=True, test_case=no-op 4: cpi_len=32, 8 range blocks, 20 cpi/thresh blk, max_deg_freedom=10, num_max_trim=1, num_min_trim=1, max_num_rfi_ev=2, mitigate_enable=True, test_case=no-op - 5: cpi_len=32, no range blocking, 20 cpi/thresh blk, max_deg_freedom=10, - num_max_trim=1, num_min_trim=1, max_num_rfi_ev=2, mitigate_enable=False, + 5: cpi_len=32, no range blocking, 20 cpi/thresh blk, max_deg_freedom=10, + num_max_trim=1, num_min_trim=1, max_num_rfi_ev=2, mitigate_enable=False, test_case=no-op, detection only """ @@ -297,11 +297,12 @@ def test_slow_time_evd( # Perform Slow-Time EVD Detection and Mitigation # Compute Eigenvalues and Eigenvectors - threshold_params = ThresholdParams([2, 10], [5, 2]) + threshold_params = ThresholdParams([2, 10], [4, 1.5]) off_diag_overlap_ratio = 0.25 diag_valid_ratio = 0.20 min_rank_frac = 0.8 rx_dynamic_range_db = 50 + bright_target_check = False swaths=None rfi_likelihood = run_slow_time_evd( @@ -314,13 +315,15 @@ def test_slow_time_evd( num_samples_rng_blk=num_samples_rng_blk, use_entire_pulse=use_entire_pulse, threshold_params=threshold_params, - num_cpi_tb=num_cpi_tb, + num_cpi_per_threshold_block=num_cpi_per_threshold_block, off_diag_overlap_ratio=off_diag_overlap_ratio, diag_valid_ratio=diag_valid_ratio, mitigate_enable=mitigate_enable, + bright_target_check=bright_target_check, min_rank_frac=min_rank_frac, rx_dynamic_range_db=rx_dynamic_range_db, swaths=swaths, + threshold_method='ev_slope', raw_data_mitigated=raw_data_mitigated, ) diff --git a/tools/imagesets/runconfigs/rslc_REE2.yaml b/tools/imagesets/runconfigs/rslc_REE2.yaml index ceeb379f6..91c6ffe7a 100644 --- a/tools/imagesets/runconfigs/rslc_REE2.yaml +++ b/tools/imagesets/runconfigs/rslc_REE2.yaml @@ -155,23 +155,65 @@ runconfig: # maximum STD is used for RFI Eigenvalue first difference # computation. max_num_rfi_ev: 2 + # Number of buffer CPIs added on each side of a detected RFI + # CPI to account for uncertainty in RFI onset and offset. + num_rfi_buffer: 3 # Number of slow-time CPIs grouped together for determining and # applying a common EV slope threshold. - num_cpi_per_threshold_block: 15 - # minimum fraction of samples required for cross-variance terms - # of sample covariance matrix - off_diag_overlap_ratio: 0.25 - # minimum fraction of samples required for covariance terms - # of sample covariance matrix - diag_valid_ratio: 0.20 + num_cpi_per_threshold_block: 12 + # Minimum fraction of samples required for cross-variance terms + # of sample covariance matrix. + off_diag_overlap_ratio: 0.20 + # Minimum fraction of samples required for diagonal covariance + # terms of sample covariance matrix. + diag_valid_ratio: 0.15 # This ratio will be used to determine the minimum number of valid Eigenvalues # required for a CPI. min_ev_valid_idx = int(np.round(min_rank_frac * cpi_len)) - 1 min_rank_frac: 0.70 - # radar platform receiver dynamic range, e.g. 50 dB. This is applied as a threshold - # to determine if the Eigenvalue under test is meaningfully signficant. + # Radar platform receiver dynamic range in dB. Applied as a threshold + # to determine if the Eigenvalue under test is meaningfully significant. rx_dynamic_range_db: 50.0 + # RFI detection threshold method. Supported values: + # - 'max_ev': Maximum Eigenvalue Estimation (default) + # - 'ev_slope': Eigenvalue slope (ST-EST) + threshold_method: max_ev + # Figure of merit threshold for RFI detection. The figure of merit is + # the STD ratio of the maximum vs. minimum Eigenvalue slopes across + # the threshold block. Values above this threshold trigger RFI detection. + rfi_fig_merit_thresh: 1.0 + # Enable bright target (e.g. strong point scatterer) rejection check. + # If enabled, threshold blocks where BOTH the condition number check + # AND the power spread check pass are skipped (not flagged as RFI). + bright_target_check: True + # Threshold (dB) on the standard deviation of the condition number + # (ratio of max to min diagonal power) across CPIs within a threshold + # block. A low value indicates a spatially stationary scene (possible + # bright target). + bt_condition_num_thresh_db: 2.0 + # Threshold (dB) on the power spread within a threshold block, + # defined as the difference between the upper and reference + # percentiles of the diagonal covariance power. A low value + # indicates a low-dynamic-range scene (possible bright target). + bt_pwr_spread_thresh_db: 6.0 + # Percentile used as the reference (baseline) power level when + # computing the power spread for bright target rejection. + pwr_ref_percentile: 50 + # Upper percentile used when computing the power spread for + # bright target rejection. + pwr_upper_percentile: 99.5 + # Upper bound of the adaptive detection margin (dB). Applied when + # the Principal Component Ratio (PCR) is below pcr_range[0]. + sig_ev_margin_upper_db: 3.0 + # Lower bound of the adaptive detection margin (dB). Applied when + # the Principal Component Ratio (PCR) is above pcr_range[1]. + sig_ev_margin_lower_db: 0.0 + # PCR range [lower, upper] used for adaptive margin interpolation. + # PCR = ev[0] / sum(ev[1:]) where ev are the sorted eigenvalues. + # The detection margin linearly interpolates between + # sig_ev_margin_upper_db and sig_ev_margin_lower_db over this range. + pcr_range: [0.1, 0.4] # Parameters used to adaptively determine RFI eigenvalue - # difference thresholds + # difference thresholds (used when threshold_method = 'max_ev') threshold_hyperparameters: # The computed sigma ratio of maximum and minimum Eigenvalue # first differences. It is a dimensionless figure of merit. diff --git a/tools/imagesets/runconfigs/rslc_REE3.yaml b/tools/imagesets/runconfigs/rslc_REE3.yaml index d2f919d89..e500130f8 100644 --- a/tools/imagesets/runconfigs/rslc_REE3.yaml +++ b/tools/imagesets/runconfigs/rslc_REE3.yaml @@ -95,23 +95,65 @@ runconfig: # maximum STD is used for RFI Eigenvalue first difference # computation. max_num_rfi_ev: 2 + # Number of buffer CPIs added on each side of a detected RFI + # CPI to account for uncertainty in RFI onset and offset. + num_rfi_buffer: 3 # Number of slow-time CPIs grouped together for determining and # applying a common EV slope threshold. - num_cpi_per_threshold_block: 15 - # minimum fraction of samples required for cross-variance terms - # of sample covariance matrix - off_diag_overlap_ratio: 0.25 - # minimum fraction of samples required for covariance terms - # of sample covariance matrix - diag_valid_ratio: 0.20 + num_cpi_per_threshold_block: 12 + # Minimum fraction of samples required for cross-variance terms + # of sample covariance matrix. + off_diag_overlap_ratio: 0.20 + # Minimum fraction of samples required for diagonal covariance + # terms of sample covariance matrix. + diag_valid_ratio: 0.15 # This ratio will be used to determine the minimum number of valid Eigenvalues # required for a CPI. min_ev_valid_idx = int(np.round(min_rank_frac * cpi_len)) - 1 min_rank_frac: 0.70 - # radar platform receiver dynamic range, e.g. 50 dB. This is applied as a threshold - # to determine if the Eigenvalue under test is meaningfully signficant. + # Radar platform receiver dynamic range in dB. Applied as a threshold + # to determine if the Eigenvalue under test is meaningfully significant. rx_dynamic_range_db: 50.0 + # RFI detection threshold method. Supported values: + # - 'max_ev': Maximum Eigenvalue Estimation (default) + # - 'ev_slope': Eigenvalue slope (ST-EST) + threshold_method: max_ev + # Figure of merit threshold for RFI detection. The figure of merit is + # the STD ratio of the maximum vs. minimum Eigenvalue slopes across + # the threshold block. Values above this threshold trigger RFI detection. + rfi_fig_merit_thresh: 1.0 + # Enable bright target (e.g. strong point scatterer) rejection check. + # If enabled, threshold blocks where BOTH the condition number check + # AND the power spread check pass are skipped (not flagged as RFI). + bright_target_check: True + # Threshold (dB) on the standard deviation of the condition number + # (ratio of max to min diagonal power) across CPIs within a threshold + # block. A low value indicates a spatially stationary scene (possible + # bright target). + bt_condition_num_thresh_db: 2.0 + # Threshold (dB) on the power spread within a threshold block, + # defined as the difference between the upper and reference + # percentiles of the diagonal covariance power. A low value + # indicates a low-dynamic-range scene (possible bright target). + bt_pwr_spread_thresh_db: 6.0 + # Percentile used as the reference (baseline) power level when + # computing the power spread for bright target rejection. + pwr_ref_percentile: 50 + # Upper percentile used when computing the power spread for + # bright target rejection. + pwr_upper_percentile: 99.5 + # Upper bound of the adaptive detection margin (dB). Applied when + # the Principal Component Ratio (PCR) is below pcr_range[0]. + sig_ev_margin_upper_db: 3.0 + # Lower bound of the adaptive detection margin (dB). Applied when + # the Principal Component Ratio (PCR) is above pcr_range[1]. + sig_ev_margin_lower_db: 0.0 + # PCR range [lower, upper] used for adaptive margin interpolation. + # PCR = ev[0] / sum(ev[1:]) where ev are the sorted eigenvalues. + # The detection margin linearly interpolates between + # sig_ev_margin_upper_db and sig_ev_margin_lower_db over this range. + pcr_range: [0.1, 0.4] # Parameters used to adaptively determine RFI eigenvalue - # difference thresholds + # difference thresholds (used when threshold_method = 'max_ev') threshold_hyperparameters: # The computed sigma ratio of maximum and minimum Eigenvalue # first differences. It is a dimensionless figure of merit.