Bandpass filter the input time series image

    ANTsR function: `frequencyFilterfMRI`

    Arguments
    ---------

    image: input time series image

    tr:    sampling time interval (inverse of sampling rate)

    lowf:  low frequency cutoff

    highf: high frequency cutoff

    order: order of the butterworth filter run using `filtfilt`

    Returns
    -------
    filtered matrix

    Example
    -------

    >>> import numpy as np
    >>> import ants
    >>> import matplotlib.pyplot as plt
    >>> brainSignal = np.random.randn( 400, 1000 )
    >>> tr = 1
    >>> filtered = ants.bandpass_filter_matrix( brainSignal, tr = tr )
    >>> nsamples = brainSignal.shape[0]
    >>> t = np.linspace(0, tr*nsamples, nsamples, endpoint=False)
    >>> k = 20
    >>> plt.plot(t, brainSignal[:,k], label='Noisy signal')
    >>> plt.plot(t, filtered[:,k], label='Filtered signal')
    >>> plt.xlabel('time (seconds)')
    >>> plt.grid(True)
    >>> plt.axis('tight')
    >>> plt.legend(loc='upper left')
    >>> plt.show()

def bandpass_filter_matrix( matrix,
    tr=1, lowf=0.01, highf=0.1, order = 3):
    """
    Bandpass filter the input time series image

    ANTsR function: `frequencyFilterfMRI`

    Arguments
    ---------

    image: input time series image

    tr:    sampling time interval (inverse of sampling rate)

    lowf:  low frequency cutoff

    highf: high frequency cutoff

    order: order of the butterworth filter run using `filtfilt`

    Returns
    -------
    filtered matrix

    Example
    -------

    >>> import numpy as np
    >>> import ants
    >>> import matplotlib.pyplot as plt
    >>> brainSignal = np.random.randn( 400, 1000 )
    >>> tr = 1
    >>> filtered = ants.bandpass_filter_matrix( brainSignal, tr = tr )
    >>> nsamples = brainSignal.shape[0]
    >>> t = np.linspace(0, tr*nsamples, nsamples, endpoint=False)
    >>> k = 20
    >>> plt.plot(t, brainSignal[:,k], label='Noisy signal')
    >>> plt.plot(t, filtered[:,k], label='Filtered signal')
    >>> plt.xlabel('time (seconds)')
    >>> plt.grid(True)
    >>> plt.axis('tight')
    >>> plt.legend(loc='upper left')
    >>> plt.show()
    """
    from scipy.signal import butter, filtfilt

    def butter_bandpass(lowcut, highcut, fs, order ):
        nyq = 0.5 * fs
        low = lowcut / nyq
        high = highcut / nyq
        b, a = butter(order, [low, high], btype='band')
        return b, a

    def butter_bandpass_filter(data, lowcut, highcut, fs, order ):
        b, a = butter_bandpass(lowcut, highcut, fs, order=order)
        y = filtfilt(b, a, data)
        return y

    fs = 1/tr   # sampling rate based on tr
    nsamples = matrix.shape[0]
    ncolumns = matrix.shape[1]
    matrixOut = matrix.copy()
    for k in range( ncolumns ):
        matrixOut[:,k] = butter_bandpass_filter(
            matrix[:,k], lowf, highf, fs, order=order )
    return matrixOut