Source code for ocmw.core.spectral

# -*- coding: utf-8 -*-

"""
Basic timeseries spectral analysis functions.
"""

# ----------------------------------------------------------------------------
#   IMPORTS
# ----------------------------------------------------------------------------
# Standard Python Dependencies
import numpy as np
from scipy.fftpack import fft, fftfreq
# Non-Standard Python Dependencies
# Local Module Dependencies
# Other Dependencies


#--------------------------------------------------------------------------
# GLOBAL CONSTANTS
#--------------------------------------------------------------------------


#--------------------------------------------------------------------------
# FUNCTION DEFINITIONS
#--------------------------------------------------------------------------



def signal_from_fft(times,freq,components,spectrum):
    """
    Reconstruct a times series from the frequency spectrum, with the option to
    separate out frequency dependent processes using a subset of the spectral 
    components.

    Parameters
    ----------
    times : numpy.array, float
        time stamps (in seconds) to reconstruct signal for.
    freq : numpy.array, float
        bin frequencies (in Hz).
    components : tuple, int
        set of frequency components to include in reconstruction.
    spectrum : numpy.array, complex, float
        complex spectral data at defined frequencies.

    Returns
    -------
    signal : numpy.array, float
        reconstructed signal at the specifcied times.

    """
    n_comp = len(components)

    for idx in np.arange(n_comp):
        index = components[idx];
        basis = np.exp(1j*2*np.pi*freq[index]*times);
        if idx == 0:
            sig_cmplx = basis * spectrum[index];
        else:
            sig_cmplx = sig_cmplx + basis * spectrum[index];
    
    signal = np.real(sig_cmplx);

    return signal




def nextpow2(n):
    """
    Find the power of two the gives n or the next highest power of 2 beyond 
    if n is not a power of 2.

    Parameters
    ----------
    n : int
        number of points to find the next power of 2 for.

    Returns
    -------
    power : int
        exponent in base two that gives n of the next highest power of 2.

    """
    power = np.ceil(np.log2(abs(n)))
    return power




def ssas_fft(X,Fs): 
    """Generates a single-sided amplitude spectrum.
    """
    L = len(X)                      # Number of samples
    T = 1.0/Fs                      # Sample period
    t = np.arange(L)*T              # Time array
    n = int(np.power(2.0,nextpow2(L)))
    Y = fft(X,n)
    f = Fs*np.arange((n/2)+1)/n
    P = np.abs(Y/n)
    a = P[0:((n//2)+1)]
    a[1:] = 2.0*a[1:]
    return f, a