#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 14 12:43:50 2025

@author: imp
"""
import numpy as np

def func_pwelch(x, hwin, N_overl, fs):
    # Compute normalization constants
    U = np.mean(hwin**2)
    L = len(hwin)
    
    # Number of segments
    ii_max = int((len(x) - L) / (L - N_overl) + 1)
    
    # Frequency bins calculation
    k = np.arange(L)
    omega = 2 * np.pi * np.fft.fftfreq(L, 1/fs)
    id_kp = np.where(k >= 0)[0]
    
    # Initialize arrays
    Imtx = np.zeros((len(id_kp), ii_max))
    Sxxhat = np.zeros(len(id_kp), dtype=complex)
    DFT2bar = np.zeros(len(id_kp), dtype=complex)
    
    for ii in range(ii_max):
        # Segment indices
        start_idx = int(ii * (L - N_overl))
        end_idx = start_idx + L
        xsegment = x[start_idx:end_idx]
        
        # Window and compute FFT
        V = np.fft.fft(xsegment * hwin)
        
        # Periodogram computation
        I = (2 / fs) * (1 / (L * U)) * np.abs(V[id_kp])**2
        
        # Store and accumulate
        Imtx[:, ii] = I
        Sxxhat += I / ii_max
        DFT2bar += np.abs(V[id_kp])**2 / ii_max
    
    # Adjust DC component
    Sxxhat[0] *= 0.5
    
    # Frequency vector
    vfp = omega[id_kp] * (fs / (2 * np.pi))
    
    # Compute final spectral estimates
    S1 = np.sum(hwin)
    S2 = np.sum(hwin**2)
    fN = fs / 2
    
    PSD = (1 / fN) * (1 / S2) * DFT2bar
    PSD[0] *= 0.5
    
    PS = 2 * (1 / S1**2) * DFT2bar
    PS[0] *= 0.5
    
    LSD = np.sqrt(PSD)
    LS = np.sqrt(PS)
    
    return {
        'vfp': vfp,
        'Sxxhat': Sxxhat,
        'PSD': PSD,
        'PS': PS,
        'LSD': LSD,
        'LS': LS
    }
