def filter_convolve(data, filters, filter_rot=False, method='scipy'):
    r"""Filter convolve

    This method convolves the input image with the wavelet filters

    Parameters
    ----------
    data : np.ndarray
        Input data, 2D array
    filters : np.ndarray
        Wavelet filters, 3D array
    filter_rot : bool, optional
        Option to rotate wavelet filters (default is 'False')
    method : str {'astropy', 'scipy'}, optional
        Convolution method (default is 'scipy')

    Returns
    -------
    np.ndarray convolved data

    Examples
    --------
    >>> from modopt.signal.wavelet import filter_convolve
    >>> x = np.arange(9).reshape(3, 3).astype(float)
    >>> y = np.arange(36).reshape(4, 3, 3).astype(float)
    >>> filter_convolve(x, y)
    array([[[  174.,   165.,   174.],
            [   93.,    84.,    93.],
            [  174.,   165.,   174.]],

           [[  498.,   489.,   498.],
            [  417.,   408.,   417.],
            [  498.,   489.,   498.]],

           [[  822.,   813.,   822.],
            [  741.,   732.,   741.],
            [  822.,   813.,   822.]],

           [[ 1146.,  1137.,  1146.],
            [ 1065.,  1056.,  1065.],
            [ 1146.,  1137.,  1146.]]])

    >>> filter_convolve(y, y, filter_rot=True)
    array([[ 14550.,  14586.,  14550.],
           [ 14874.,  14910.,  14874.],
           [ 14550.,  14586.,  14550.]])

    """

    if filter_rot:
        return np.sum([convolve(coef, f, method=method) for coef, f in
                       zip(data, rotate_stack(filters))], axis=0)

    else:
        return np.array([convolve(data, f, method=method) for f in filters])