Untitled

def compute_features(I, Q, fs=2e6):
    """
    Computes various statistical, spectral, and time-frequency features from I/Q signals.
    fs: Sampling frequency (default 2 MHz for ADS-B)
    """
    features = {}
    
    # Compute amplitude and phase
    amplitude = np.sqrt(I**2 + Q**2)
    phase = np.arctan2(Q, I)
    
    # === 1. Statistical Features ===
    features['mean_I'] = np.mean(I)
    features['mean_Q'] = np.mean(Q)
    features['var_I'] = np.var(I)
    features['var_Q'] = np.var(Q)
    features['skew_I'] = stats.skew(I)
    features['skew_Q'] = stats.skew(Q)
    features['kurt_I'] = stats.kurtosis(I)
    features['kurt_Q'] = stats.kurtosis(Q)
    features['iq_correlation'] = np.corrcoef(I, Q)[0, 1]
    
    # === 2. Amplitude & Phase Features ===
    features['mean_amplitude'] = np.mean(amplitude)
    features['var_amplitude'] = np.var(amplitude)
    features['mean_phase'] = np.mean(phase)
    features['var_phase'] = np.var(phase)
    
    # === 3. Spectral Features (FFT) ===
    fft_I = np.fft.fft(I)
    fft_Q = np.fft.fft(Q)
    freq = np.fft.fftfreq(len(I), d=1/fs)
    
    power_spectrum = np.abs(fft_I)**2 + np.abs(fft_Q)**2
    features['spectral_centroid'] = np.sum(freq * power_spectrum) / np.sum(power_spectrum)
    features['spectral_bandwidth'] = np.sqrt(np.sum(((freq - features['spectral_centroid'])**2) * power_spectrum) / np.sum(power_spectrum))
    
    # === 4. Time-Frequency Representations (Spectrogram) ===
    f, t, Sxx = signal.spectrogram(amplitude, fs=fs, nperseg=256)
    features['spectrogram_mean'] = np.mean(Sxx)
    features['spectrogram_variance'] = np.var(Sxx)

    # === 5. Wavelet Transform Features ===
    coeffs, _ = pywt.cwt(amplitude, np.arange(1, 128), 'morl')
    features['wavelet_energy'] = np.sum(coeffs**2)

    # === 6. MFCCs (Mel Frequency Cepstral Coefficients) ===
    mfccs = librosa.feature.mfcc(y=amplitude, sr=int(fs), n_mfcc=13)
    features['mfcc_mean'] = np.mean(mfccs)
    features['mfcc_var'] = np.var(mfccs)
    
    return features