Commit d379154c authored by Andrew Quinn's avatar Andrew Quinn
Browse files

Add start to multi-channel wavelet and wavelet-csd function

parent 8ec9ee8b
Loading
Loading
Loading
Loading
+24 −5
Original line number Diff line number Diff line
@@ -28,16 +28,20 @@ def morlet(x, freqs, sample_rate, window_len=4, ncycles=5, ret_basis=False,
    2D array
        Array containing morlet wavelet transformed data [nfreqs x nsamples]
    """
    orig_dim = x.ndim
    if orig_dim == 1:
        x = x[np.newaxis, :]

    cwt = np.zeros((len(freqs), *x.shape[:]), dtype=complex)
    cwt = np.zeros((x.shape[0], len(freqs), x.shape[1]), dtype=complex)

    # Get morlet basis
    mlt = get_morlet_basis(freqs, ncycles, sample_rate, normalise=normalise)

    for jj in range(x.shape[0]):
        for ii in range(len(freqs)):
        a = signal.convolve(x, mlt[ii].real, mode='same', method='fft')
        b = signal.convolve(x, mlt[ii].imag, mode='same', method='fft')
        cwt[ii, ...] = a+1j*b
            a = signal.convolve(x[jj, :], mlt[ii].real, mode='same', method='fft')
            b = signal.convolve(x[jj, :], mlt[ii].imag, mode='same', method='fft')
            cwt[jj, ii, :] = a+1j*b

    if ret_mode == 'power':
        cwt = np.power(np.abs(cwt), 2)
@@ -46,12 +50,27 @@ def morlet(x, freqs, sample_rate, window_len=4, ncycles=5, ret_basis=False,
    elif ret_mode != 'complex':
        raise ValueError("'ret_mode not recognised, please use one of {'power','amplitude','complex'}")

    if orig_dim == 1:
        cwt = cwt[0, ...]

    if ret_basis:
        return cwt, mlt
    else:
        return cwt


def wavelet_csd(x, freqs, sample_rate):

    wt = morlet(x, freqs, sample_rate, ret_mode='complex')
    S = np.zeros((wt.shape[0], wt.shape[0], wt.shape[1], wt.shape[2]), dtype=complex)
    for ii in range(wt.shape[1]):
        for jj in range(wt.shape[2]):
            S[:, :, ii, jj] = np.dot(wt[:, ii, jj, np.newaxis], wt[np.newaxis, :, ii, jj].conj())

    return S



def get_morlet_basis(freq, ncycles, sample_rate, normalise=False, win_len=5):
    """Compute a morlet wavelet basis set based on specified parameters.