import numpy as np
from scipy.fft import fft, fftfreq
from scipy.stats import skew, kurtosis, entropy
from scipy.signal import hamming

def calculate_frequency_features(segmented_signals, fs):
    # Initialize arrays to store frequency-domain features
    num_segments = len(segmented_signals)
    mean_frequency = np.zeros(num_segments)
    median_frequency = np.zeros(num_segments)
    peak_frequency = np.zeros(num_segments)
    bandwidth = np.zeros(num_segments)
    power_spectral_density = [None] * num_segments
    total_power = np.zeros(num_segments)
    spectral_entropy = np.zeros(num_segments)
    spectral_flatness = np.zeros(num_segments)
    spectral_centroid = np.zeros(num_segments)
    spectral_skewness = np.zeros(num_segments)
    spectral_kurtosis = np.zeros(num_segments)

    # Compute frequency-domain features for each segment
    for i, signal in enumerate(segmented_signals):
        # Apply windowing to the segment
        windowed_signal = signal * hamming(len(signal))
        
        # Compute the Fast Fourier Transform (FFT)
        fft_signal = fft(windowed_signal)
        freqs = fftfreq(len(windowed_signal), d=1/fs)

        # Calculate the amplitude spectrum (absolute values)
        amplitude_spectrum = np.abs(fft_signal)

        # Find peak frequency
        peak_frequency[i] = freqs[np.argmax(amplitude_spectrum)]

        # Calculate mean frequency
        mean_frequency[i] = np.sum(freqs * amplitude_spectrum) / np.sum(amplitude_spectrum)

        # Calculate median frequency
        sorted_indices = np.argsort(amplitude_spectrum)
        sorted_amp_spectrum = amplitude_spectrum[sorted_indices]
        sorted_freqs = freqs[sorted_indices]
        cumulative_sum = np.cumsum(sorted_amp_spectrum)
        median_frequency_index = np.argmax(cumulative_sum >= 0.5 * np.sum(sorted_amp_spectrum))
        median_frequency[i] = sorted_freqs[median_frequency_index]

        # Calculate bandwidth (using Full Width at Half Maximum, FWHM)
        half_max_amplitude = max(amplitude_spectrum) / 2
        mask = amplitude_spectrum > half_max_amplitude
        fwhm_points = np.where(mask)[0]
        bandwidth[i] = freqs[fwhm_points[-1]] - freqs[fwhm_points[0]]

        # Calculate Power Spectral Density (PSD)
        psd = (amplitude_spectrum ** 2) / (2 * np.pi)
        power_spectral_density[i] = psd

        # Calculate total power
        total_power[i] = np.sum(psd)

        # Calculate spectral entropy
        spectral_entropy[i] = entropy(amplitude_spectrum)

        # Calculate spectral flatness
        spectral_flatness[i] = np.exp(np.mean(np.log(amplitude_spectrum))) / np.mean(amplitude_spectrum)

        # Calculate spectral centroid
        spectral_centroid[i] = np.sum(freqs * amplitude_spectrum) / np.sum(amplitude_spectrum)

        # Calculate spectral skewness
        spectral_skewness[i] = skew(amplitude_spectrum)

        # Calculate spectral kurtosis
        spectral_kurtosis[i] = kurtosis(amplitude_spectrum)

    return mean_frequency, median_frequency, peak_frequency, bandwidth, power_spectral_density, total_power, spectral_entropy, spectral_flatness, spectral_centroid, spectral_skewness, spectral_kurtosis

# Given time signal data
time_signal = [
    0.0546595813338489, 0.0580966797834281, 0.0615337782330085, 0.0649708766825892,
    0.0684079751321685, 0.0718450735817489, 0.0752821720313284, 0.0787192704809090,
    0.0821563689304883, 0.0855934673800687, 0.0890305658296480, 0.0924676642792286,
    0.0959047627288092, 0.0993418611783886, 0.102778959627969, 0.106216058077549,
    0.109653156527129, 0.113090254976708, 0.116527353426289, 0.119964451875868,
    0.123401550325449, 0.126838648775029, 0.130275747224608, 0.133712845674189,
    0.137149944123768, 0.140587042573349, 0.144024141022928, 0.147461239472509,
    0.150898337922090, 0.154335436371669, 0.157772534821249, 0.161209633270829,
    0.164646731720409, 0.168083830169988, 0.171520928619569, 0.174958027069148,
    0.178395125518729, 0.181832223968310, 0.185269322417889, 0.188706420867469,
    0.192143519317049, 0.195580617766629, 0.199017716216209, 0.202454814665789,
    0.205891913115368, 0.209329011564949, 0.212766110014530, 0.216203208464109,
    0.219640306913689, 0.223077405363269, 0.226514503812849
]

# Sampling frequency
fs = 50  # Assuming 1 Hz

# Calculate frequency-domain features
mean_frequency, median_frequency, peak_frequency, bandwidth, power_spectral_density, total_power, spectral_entropy, spectral_flatness, spectral_centroid, spectral_skewness, spectral_kurtosis = calculate_frequency_features([time_signal], fs)

# Output results
print("Mean Frequency:", mean_frequency[0])
print("Median Frequency:", median_frequency[0])
print("Peak Frequency:", peak_frequency[0])
print("Bandwidth:", bandwidth[0])
print("Total Power:", total_power[0])
print("Spectral Entropy:", spectral_entropy[0])
print("Spectral Flatness:", spectral_flatness[0])
print("Spectral Centroid:", spectral_centroid[0])
print("Spectral Skewness:", spectral_skewness[0])
print("Spectral Kurtosis:", spectral_kurtosis[0])