powerspectra.py

"""
Estimates power spectra from proxy record dataframes
====================================================

"""
# Last modified: Wed Nov 04, 2020  01:48pm
#
# Copyright (C) 2020  Bedartha Goswami <bedartha.goswami@uni-tuebingen.de>
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.

# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.

# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
# -----------------------------------------------------------------------------


import pandas as pd
import numpy as np

from scipy.interpolate import interp1d
from scipy.signal import welch



def get_powerspectra(proxyrecords):
    """
    Computes power spectra for COPRA ensemble and saves to file
    """
    pra = proxyrecord["Age"]
    prp = proxyrecord[proxyrecord.columns[1:]].to_numpy().T
    ES = prp.shape[0]

    # loop over proxy records and estimate psd
    freq, psd = [], []
    for i in range(ES):
        f, p = welch(prp[i])
        freq.append(f)
        psd.append(p)
    freq = np.array(freq)
    psd = np.array(psd)

    # create DataFrame
    names = ["Frequency"]
    names.extend(["PSD %d" % (i + 1) for i in range(ES)])
    powerspectra = pd.DataFrame(np.c_[freq[0].T, psd.T], columns=names)

    return powerspectra