def calculate_psd(container, twin, fargs, nperseg, normalize):
    """
    Calculate power spectral density and return as numpy arrays

    Parameters
    ----------
    container : dict
        TimeSeries objects
    twin : tuple
        Time window. Time series are cropped to time window before estimating psd.
    fargs : tuple
        Filter arguments. Time series are filtered before estimating psd.
    nperseg : int
        Size of segments used for ustimating PSD using Welch's method.
    normalize : bool
        Normalize power spectral density on maximum density.

    Returns
    -------
    dict
        Frequency versus power spectral density
    """
    container_out = dict()

    for name, ts in container.items():
        # resampling to average time step for robustness (necessary for series with varying time step)
        f, s = ts.psd(twin=twin, filterargs=fargs, resample=ts.dt, taperfrac=0.1, nperseg=nperseg, normalize=normalize)
        container_out[name] = tuple([f, s])

    return container_out