"""
Functions that act on the Spatial Fourier Coefficients
`BEMA`
   Bandwidth Extension for Microphone Arrays

`FFT`
   (Fast) Fourier Transform
`iFFT`
   Inverse (Fast) Fourier Transform

`spatFT`
   Spatial Fourier Transform
`spatFT_RT`
   Spatial Fourier Transform for real-time application
`spatFT_LSF`
   Spatial Fourier Transform by least square fit to provided data
`iSpatFT`
   Fast Inverse Spatial Fourier Transform

`rfi`
   Radial Filter Improvement

`plane_wave_decomp`
   Plane Wave Decomposition
"""
import sys

import numpy as _np
from scipy.linalg import lstsq
from scipy.signal import butter, fftconvolve, sosfreqz

from .io import SphericalGrid
from .sph import besselj, hankel1, sph_harm_all


# noinspection PyUnusedLocal
def BEMA(Pnm, center_sig, dn, transition, avg_band_width=1, fade=True, max_order=None):
    """BEMA Spatial Anti-Aliasing

    Parameters
    ----------
    Pnm : array_like
       Sound field SH spatial Fourier coefficients
    center_sig : array_like
      center microphone in shape [0, NFFT]
    dn : array_like
       Radial filters for the current array configuration
    transition : int
       Highest stable bin, approx: transition = (NFFT/fs) * (N*c)/(2*pi*r)
    avg_band_width : int, optional
       Averaging Bandwidth in oct [Default: 1]
    fade : bool, optional
       Fade over if True, else hard cut [Default: True]
    max_order : int, optional
       Maximum transform order [Default: highest available order]

    Returns
    -------
    Pnm : array_like
       BEMA optimized sound field SH coefficients

    References
    ----------
    .. [1] B. Bernsch├╝tz, "Bandwidth Extension for Microphone Arrays",
       AES Convention 2012, Convention Paper 8751, 2012. http://www.aes.org/e-lib/browse.cfm?elib=16493
    """

    if not max_order:
        max_order = int(_np.sqrt(Pnm.shape[0] - 1))  # SH order

    transition = int(_np.floor(transition))

    # computing spatio-temporal Image Imn
    Imn = _np.zeros((Pnm.shape[0], 1), dtype=complex)  # pre-allocating the spectral image matrix
    start_avg_bin = int(_np.floor(transition / (_np.power(2, avg_band_width))))  # first bin for averaging

    modeCnt = 0
    avgPower = 0
    for n in range(0, max_order + 1):  # sh orders
        for m in range(0, 2 * n + 1):  # modes
            for inBin in range(start_avg_bin - 1, transition):  # bins
                synthBin = Pnm[modeCnt, inBin] * dn[n, inBin]
                Imn[modeCnt, 0] += synthBin
                avgPower += _np.abs(synthBin) ** 2
            modeCnt += 1

    # energy normalization (Eq. 9)
    energySum = _np.sum(_np.power(_np.abs(Imn), 2))
    normFactor = avgPower / (transition - start_avg_bin + 1)
    Imn *= _np.sqrt(normFactor / energySum)

    # normalize center signal to its own average level in the source band. (Eq. 13)
    center_sig[_np.where(center_sig == 0)] = 1e-12
    sq_avg = _np.sqrt(_np.mean(_np.power(_np.abs(center_sig[:, (start_avg_bin - 1):transition]), 2)))
    center_sig = _np.multiply(center_sig, (1 / sq_avg))

    Pnm_synth = _np.zeros(Pnm.shape, dtype=complex)
    modeCnt = 0
    for n in range(0, max_order + 1):
        for m in range(0, 2 * n + 1):
            for inBin in range(start_avg_bin - 1, Pnm_synth.shape[1]):
                Pnm_synth[modeCnt, inBin] = Imn[modeCnt, 0] * (1 / dn[n, inBin]) * center_sig[0, inBin]
            modeCnt += 1

    # Phase correction (Eq. 16)
    phaseOffset = _np.angle(Pnm[0, transition] / Pnm_synth[0, transition])
    Pnm_synth *= _np.exp(1j * phaseOffset)

    # Merge original and synthetic SH coefficients. (Eq. 14)
    Pnm = _np.concatenate((Pnm[:, 0:(transition - 1)], Pnm_synth[:, (transition - 1):]), axis=1)

    # Fade in of synthetic SH coefficients. (Eq. 17)
    if fade:
        Pnm_fade0 = Pnm[:, (start_avg_bin - 1):(transition - 1)]
        Pnm_fadeS = Pnm_synth[:, (start_avg_bin - 1):(transition - 1)]

        fadeUp = _np.linspace(0, 1, Pnm_fade0.shape[1])
        fadeDown = _np.flip(fadeUp)

        Pnm_fade0 *= fadeDown
        Pnm_fadeS *= fadeUp
        Pnm[:, start_avg_bin - 1:transition - 1] = Pnm_fade0 + Pnm_fadeS

    return Pnm


def FFT(time_signals, fs=None, NFFT=None, oversampling=1, first_sample=0, last_sample=None, calculate_freqs=True):
    """Real-valued Fast Fourier Transform.

    Parameters
    ----------
    time_signals : TimeSignal/tuple/object
       Time-domain signals to be transformed.
    fs : int, optional
       Sampling frequency - only optional no frequency vector should be calculated or if a TimeSignal or tuple/array
       containing fs is passed
    NFFT : int, optional
       Number of frequency bins. Resulting array will have size NFFT//2+1 Default: Next power of 2
    oversampling : int, optional
       Oversample the incoming signal to increase frequency resolution [Default: 1]
    first_sample : int, optional
       First time domain sample to be included. [Default: 0]
    last_sample : int, optional
       Last time domain sample to be included. [Default: -1]
    calculate_freqs : bool, optional
       Calculate frequency scale if True, else return only spectrum. [Default: True]

    Returns
    -------
    fftData : array_like
       One-sided frequency domain spectrum
    f : array_like, optional
       Vector of frequency bins of one-sided spectrum, in case of calculate_freqs

    Notes
    -----
    An oversampling*NFFT point Fourier Transform is applied to the time domain data, where NFFT is the next power of
    two of the number of samples. Time-windowing can be used by providing a first_sample and last_sample index.
    """
    try:
        signals = time_signals.signal
        fs = time_signals.fs
    except AttributeError:
        if fs is None and calculate_freqs:
            raise ValueError('No valid signal found. Either pass an io.TimeSignal, a tuple/array containing the signal '
                             'and the sampling frequency or use the fs argument.')
        else:
            signals = time_signals

    signals = _np.atleast_2d(signals)

    if last_sample is None:  # assign lastSample to length of signals if not provided
        last_sample = signals.shape[1]

    if oversampling < 1:
        raise ValueError('oversampling must be >= 1.')

    if last_sample < first_sample or last_sample > signals.shape[1]:
        raise ValueError('lastSample must be between firstSample and nSamples.')

    if first_sample < 0 or first_sample > last_sample:
        raise ValueError('firstSample must be between 0 and lastSample.')

    total_samples = last_sample - first_sample
    signals = signals[:, first_sample:last_sample]

    if not NFFT:
        NFFT = int(2 ** _np.ceil(_np.log2(total_samples)))

    fftData = _np.fft.rfft(signals, NFFT * oversampling, 1)
    if not calculate_freqs:
        return fftData

    f = _np.fft.rfftfreq(NFFT * oversampling, d=1 / fs)
    return fftData, f


def spatFT(data, position_grid, order_max=10, spherical_harmonic_bases=None):
    """Spatial Fourier Transform

    Parameters
    ----------
    data : array_like
       Data to be transformed, with signals in rows and frequency bins in columns
    position_grid : array_like or io.SphericalGrid
       Azimuths/Colatitudes/Gridweights of spatial sampling points
    order_max : int, optional
       Maximum transform order [Default: 10]
    spherical_harmonic_bases : array_like, optional
       Spherical harmonic base coefficients (not yet weighted by spatial sampling grid) [Default: None]

    Returns
    -------
    Pnm : array_like
       Spatial Fourier Coefficients with nm coeffs in rows and FFT bins in columns

    Notes
    -----
    In case no weights in spatial sampling grid are given, the pseudo inverse of the SH bases is computed according to
    Eq. 3.34 in [1].

    References
    ----------
    .. [1] Boaz Rafaely: Fundamentals of spherical array processing. In. Springer topics in signal processing.
           Benesty, J.; Kellermann, W. (Eds.), Springer, Heidelberg et al. (2015).
    """
    data = _np.atleast_2d(data)
    position_grid = SphericalGrid(*position_grid)

    # Re-generate spherical harmonic bases if they were not provided or their order is too small
    if (spherical_harmonic_bases is None or
            spherical_harmonic_bases.shape[0] < data.shape[0] or
            spherical_harmonic_bases.shape[1] < (order_max + 1) ** 2):
        spherical_harmonic_bases = sph_harm_all(order_max, position_grid.azimuth, position_grid.colatitude)

    if position_grid.weight is None:
        # calculate pseudo inverse in case no spatial sampling point weights are given
        spherical_harmonic_weighted = _np.linalg.pinv(spherical_harmonic_bases)
    else:
        # apply spatial sampling point weights in case they are given
        spherical_harmonic_weighted = (_np.conj(spherical_harmonic_bases).T * (4 * _np.pi * position_grid.weight))

    return spatFT_RT(data, spherical_harmonic_weighted)


def spatFT_RT(data, spherical_harmonic_weighted):
    """Spatial Fourier Transform for real-time application, otherwise use `spatFT()` for more more convenience and
    flexibility.

    Parameters
    ----------
    data : array_like
       Data to be transformed, with signals in rows and frequency bins in columns
    spherical_harmonic_weighted : array_like
       Spherical harmonic base coefficients (already weighted by spatial sampling grid)

    Returns
    -------
    Pnm : array_like
       Spatial Fourier Coefficients with nm coeffs in rows and FFT bins in columns
    """
    return _np.dot(spherical_harmonic_weighted, data)


def iSpatFT(spherical_coefficients, position_grid, order_max=None, spherical_harmonic_bases=None):
    """Inverse spatial Fourier Transform

    Parameters
    ----------
    spherical_coefficients : array_like
       Spatial Fourier coefficients with columns representing frequency bins
    position_grid : array_like or io.SphericalGrid
       Azimuth/Colatitude angles of spherical coefficients
    order_max : int, optional
       Maximum transform order [Default: highest available order]
    spherical_harmonic_bases : array_like, optional
       Spherical harmonic base coefficients (not yet weighted by spatial sampling grid) [Default: None]

    Returns
    -------
    P : array_like
       Sound pressures with frequency bins in columns and angles in rows
    """
    position_grid = SphericalGrid(*position_grid)
    spherical_coefficients = _np.atleast_2d(spherical_coefficients)
    number_of_coefficients = spherical_coefficients.shape[0]

    # todo: check coeffs and bases length correspond with order_max
    if order_max is None:
        order_max = int(_np.sqrt(number_of_coefficients) - 1)

    # Re-generate spherical harmonic bases if they were not provided or their order is too small
    if (spherical_harmonic_bases is None
            or spherical_harmonic_bases.shape[1] < number_of_coefficients
            or spherical_harmonic_bases.shape[1] != position_grid.azimuths.size):
        spherical_harmonic_bases = sph_harm_all(order_max, position_grid.azimuth, position_grid.colatitude)

    return _np.dot(spherical_harmonic_bases, spherical_coefficients)


def spatFT_LSF(data, position_grid, order_max=10, spherical_harmonic_bases=None):
    """Spatial Fourier Transform by least square fit to provided data

    Parameters
    ----------
    data : array_like, complex
       Data to be fitted to
    position_grid : array_like, or io.SphericalGrid
       Azimuth / colatitude data locations
    order_max: int, optional
       Maximum transform order [Default: 10]
    spherical_harmonic_bases : array_like, optional
       Spherical harmonic base coefficients (not yet weighted by spatial sampling grid) [Default: None]

    Returns
    -------
    coefficients: array_like, float
       Fitted spherical harmonic coefficients (indexing: n**2 + n + m + 1)
    """
    # Re-generate spherical harmonic bases if they were not provided or their order is too small
    if (spherical_harmonic_bases is None or
            spherical_harmonic_bases.shape[0] < data.shape[0] or
            spherical_harmonic_bases.shape[1] < (order_max + 1) ** 2):
        position_grid = SphericalGrid(*position_grid)
        spherical_harmonic_bases = sph_harm_all(order_max, position_grid.azimuth, position_grid.colatitude)

    return lstsq(spherical_harmonic_bases, data)[0]


def plane_wave_decomp(order, wave_direction, field_coeffs, radial_filter, weights=None):
    """Plane wave decomposition

    Parameters
    ----------
    order : int
       Decomposition order
    wave_direction : array_like
       Direction of plane wave as [azimuth, colatitude] pair. io.SphericalGrid is used internally
    field_coeffs : array_like
       Spatial fourier coefficients
    radial_filter : array_like
       Radial filters
    weights : array_like, optional
       Weighting function. Either scalar, one per directions or of dimension (nKR_bins x  nDirections). [Default: None]

    Returns
    -------
    Y : matrix of floats
       Matrix of the decomposed wavefield with kr bins in rows
    """
    wave_direction = SphericalGrid(*wave_direction)
    number_of_angles = wave_direction.azimuth.size
    field_coeffs = _np.atleast_2d(field_coeffs)
    radial_filter = _np.atleast_2d(radial_filter)

    if field_coeffs.shape[0] == 1:
        field_coeffs = field_coeffs.T
    if radial_filter.shape[0] == 1:
        radial_filter = radial_filter.T

    NMDeliveredSize, FFTBlocklengthPnm = field_coeffs.shape
    Ndn, FFTBlocklengthdn = radial_filter.shape

    if FFTBlocklengthdn != FFTBlocklengthPnm:
        raise ValueError('FFT Blocksizes of field coefficients (Pnm) and radial filter (dn) are not consistent.')

    max_order = int(_np.floor(_np.sqrt(NMDeliveredSize) - 1))
    if order > max_order:
        raise ValueError(
            f'The provided coefficients deliver a maximum order of {max_order} but order {order} was requested.')

    gaincorrection = 4 * _np.pi / ((order + 1) ** 2)

    if weights is not None:
        weights = _np.asarray(weights)
        if weights.ndim == 1:
            number_of_weights = weights.size
            if (number_of_weights != FFTBlocklengthPnm and number_of_weights != number_of_angles
                    and number_of_weights != 1):
                raise ValueError('Weights is not a scalar nor consistent with shape of the field coefficients (Pnm).')
            if number_of_weights == number_of_angles:
                weights = _np.broadcast_to(weights, (FFTBlocklengthPnm, number_of_angles)).T
        else:
            if weights.shape != (number_of_angles, FFTBlocklengthPnm):
                raise ValueError('Weights is not a scalar nor consistent with shape of field coefficients (Pnm) and '
                                 'radial filter (dn).')
    else:
        weights = 1

    sph_harms = sph_harm_all(order, wave_direction.azimuth, wave_direction.colatitude)
    filtered_coeffs = field_coeffs * _np.repeat(radial_filter, _np.r_[:order + 1] * 2 + 1, axis=0)
    OutputArray = _np.dot(sph_harms, filtered_coeffs)
    OutputArray = _np.multiply(OutputArray, weights)
    return OutputArray * gaincorrection


# noinspection PyUnusedLocal
def rfi(dn, kernelSize=512, highPass=0.0):
    """R/F/I Radial Filter Improvement

    Parameters
    ----------
    dn : array_like
       Analytical frequency domain radial filters (e.g. gen.radial_filter_fullspec())
    kernelSize : int, optional
       Target filter kernel size [Default: 512]
    highPass : float, optional
       Highpass Filter from 0.0 (off) to 1.0 (maximum kr) [Default: 0.0]

    Returns
    -------
    dn : array_like
       Improved radial filters
    kernelSize : int
       Filter kernel size (total)
    latency : float
       Approximate signal latency due to the filters

    Note
    ----
    This function improves the FIR radial filters from gen.radial_filter_fullspec(). The filters are made causal and are
    windowed in time domain. The DC components are estimated. The R/F/I module should always be inserted to the filter
    path when treating measured data even if no use is made of the included kernel downscaling or highpass filters.

    Do NOT use R/F/I for single open sphere filters (e.g.simulations).

    IMPORTANT
       Remember to choose a kernel size being large enough to cover all filter latencies and response slopes. Otherwise
       undesired cyclic convolution artifacts may appear in the output signal.

    HIGHPASS
       If HPF is on (0<highPass<=1) the radial filters and HPF share the available taps and the latency keeps constant.
       Be careful when using very small kernel sizes because since there might be too few taps. Observe the filters by
       plotting their spectra and impulse responses!
       > Be very careful if NFFT/max(kr) < 25
       > Do not use R/F/I if NFFT/max(kr) < 15
    """
    dn = _np.atleast_2d(dn)
    sourceKernelSize = (dn.shape[-1] - 1) * 2

    if kernelSize > sourceKernelSize:
        raise ValueError('Kernelsize greater than radial filters. Extension of kernelsize not yet implemented.')
        # TODO: Implement kernelsize extension

    # in case highpass should be applied, half both individual kernel sizes
    endKernelSize = kernelSize
    if highPass:
        kernelSize //= 2

    # DC-component estimation
    dn_diff = _np.abs(dn[:, 1] / dn[:, 2])
    oddOrders = range(1, dn.shape[0], 2)
    dn[oddOrders, 0] = -1j * dn[oddOrders, 1] * 2 * (sourceKernelSize / kernelSize) * dn_diff[oddOrders]

    # transform into time domain
    dn_ir = iFFT(dn)

    # make filters causal by circular shift
    dn_ir = _np.roll(dn_ir, dn_ir.shape[-1] // 2, axis=-1)

    # downsize kernel
    latency = kernelSize / 2
    mid = dn_ir.shape[-1] / 2
    dn_ir = dn_ir[:, int(mid - latency):int(mid + latency)]

    # apply Hanning / cosine window
    # 0.5 + 0.5 * _np.cos(2 * _np.pi * (_np.arange(0, kernelSize) - ((kernelSize - 1) / 2)) / (kernelSize - 1))
    dn_ir *= _np.hanning(kernelSize)

    # transform into one-sided spectrum
    dn = FFT(dn_ir, NFFT=endKernelSize, calculate_freqs=False)

    # calculate high pass (need to be zero phase since radial filters are already linear phase)
    if 0 < highPass <= 1:
        # by designing an IIR filter and taking the magnitude response
        # calculate IIR filter coefficients
        hp_sos = butter(8, highPass, btype='highpass', output='sos', analog=False)
        # calculate "equivalent" zero phase FIR filter one-sided spectrum
        _, HP = sosfreqz(hp_sos, worN=kernelSize // 2 + 1, whole=False)
        HP[_np.isnan(HP)] = 0  # prevent NaNs
        # make filter zero phase by circular shift
        hp = _np.roll(iFFT(HP), kernelSize // 2, axis=-1)
        latency += kernelSize / 2
        # append zeros in time domain
        HP = FFT(hp, NFFT=endKernelSize, calculate_freqs=False)

        # # by designing an FIR filter with FIRLS
        # from scipy.signal import firls
        # STOP_BAND_ATTENUATION = -30
        # # calculate FIR linear phase filter coefficients
        # bands = [0, highPass/2, highPass, 1]
        # desired = 10 ** (_np.array([STOP_BAND_ATTENUATION, STOP_BAND_ATTENUATION, 0, 0]) / 20)  # magnitudes as linear
        # hp = firls(kernelSize - 1, bands=bands, desired=desired)
        # # make filter zero phase by circular shift
        # hp = _np.roll(hp, hp.shape[-1] // 2, axis=-1)
        # # transform into one-sided spectrum
        # HP = FFT(hp, calculate_freqs=False)

        # # by designing an FIR filter with FIRWIN
        # from scipy.signal import firwin
        # # calculate FIR linear phase filter coefficients
        # hp = firwin(kernelSize - 1, highPass, pass_zero=False, window='hamming')
        # # make filter zero phase by circular shift
        # hp = _np.roll(hp, hp.shape[-1] // 2, axis=-1)
        # # transform into one-sided spectrum
        # HP = FFT(hp, calculate_freqs=False)

        # apply high pass
        dn *= HP

    return dn, kernelSize, latency


def sfe(Pnm_kra, kra, krb, problem='interior'):
    """ S/F/E Sound Field Extrapolation. CURRENTLY WIP

    Parameters
    ----------
    Pnm_kra : array_like
       Spatial Fourier Coefficients (e.g. from spatFT())
    kra,krb : array_like
       k * ra/rb vector
    problem : string{'interior', 'exterior'}
       Select between interior and exterior problem [Default: interior]
    """
    if kra.shape[1] != Pnm_kra.shape[1] or kra.shape[1] != krb.shape[1]:
        raise ValueError('FFTData: Complex Input Data expected.')

    FCoeff = Pnm_kra.shape[0]
    N = int(_np.floor(_np.sqrt(FCoeff) - 1))

    nvector = _np.zeros(FCoeff)
    IDX = 1

    for n in range(0, N + 1):
        for m in range(-n, n + 1):
            nvector[IDX] = n
            IDX += 1

    nvector = _np.tile(nvector, (1, Pnm_kra.shape[1]))
    kra = _np.tile(kra, (FCoeff, 1))
    krb = _np.tile(krb, (FCoeff, 1))

    if problem == 'interior':
        jn_kra = _np.sqrt(_np.pi / (2 * kra)) * besselj(nvector + 5, kra)
        jn_krb = _np.sqrt(_np.pi / (2 * krb)) * besselj(nvector + 5, krb)
        exp = jn_krb / jn_kra

        if _np.any(_np.abs(exp) > 1e2):  # 40dB
            print('WARNING: Extrapolation might be unstable for one or more frequencies/orders!', file=sys.stderr)

    elif problem == 'exterior':
        hn_kra = _np.sqrt(_np.pi / (2 * kra)) * hankel1(nvector + 0.5, kra)
        hn_krb = _np.sqrt(_np.pi / (2 * krb)) * hankel1(nvector + 0.5, krb)
        exp = hn_krb / hn_kra
    else:
        raise ValueError(
            f'Problem selector {problem} not recognized. Please either choose "interior" [Default] or "exterior".')

    return Pnm_kra * exp.T


def iFFT(Y, output_length=None, window=False):
    """ Inverse real-valued Fourier Transform

    Parameters
    ----------
    Y : array_like
       Frequency domain data [Nsignals x Nbins]
    output_length : int, optional
       Length of returned time-domain signal (Default: 2 x len(Y) + 1)
    window : boolean, optional
       Window applied to the resulting time-domain signal

    Returns
    -------
    y : array_like
       Reconstructed time-domain signal
    """
    Y = _np.atleast_2d(Y)
    y = _np.fft.irfft(Y, n=output_length)

    if window:
        no_of_samples = y.shape[-1]

        if window == 'hann':
            window_array = _np.hanning(no_of_samples)
        elif window == 'hamming':
            window_array = _np.hamming(no_of_samples)
        elif window == 'blackman':
            window_array = _np.blackman(no_of_samples)
        elif window == 'kaiser':
            window_array = _np.kaiser(no_of_samples, 3)
        else:
            raise ValueError('Selected window must be one of hann, hamming, blackman or kaiser')

        y *= window_array

    return y


# noinspection PyUnusedLocal
def wdr(Pnm, xAngle, yAngle, zAngle):
    """W/D/R Wigner-D Rotation - NOT YET IMPLEMENTED

    Parameters
    ----------
    Pnm : array_like
       Spatial Fourier coefficients
    xAngle, yAngle, zAngle : float
       Rotation angle around the x/y/z-Axis

    Returns
    -------
    PnmRot: array_like
       Rotated spatial Fourier coefficients
    """
    print('!WARNING. Wigner-D Rotation is not yet implemented. Continuing with un-rotated coefficients!',
          file=sys.stderr)
    return Pnm


def convolve(A, B, FFT=None):
    """ Convolve two arrays A & B row-wise. One or both can be one-dimensional for SIMO/SISO convolution

    Parameters
    ----------
    A, B: array_like
       Data to perform the convolution on of shape [Nsignals x NSamples]
    FFT: bool, optional
       Selects whether time or frequency domain convolution is applied. Default: On if Nsamples > 500 for both

    Returns
    -------
    out: array
       Array containing row-wise, linear convolution of A and B
    """
    A = _np.atleast_2d(A)
    B = _np.atleast_2d(B)

    N_sigA, L_sigA = A.shape
    N_sigB, L_sigB = B.shape

    if FFT is None and (L_sigA > 500 and L_sigB > 500):
        FFT = True
    else:
        FFT = False

    if (N_sigA != N_sigB) and not (N_sigA == 1 or N_sigB == 1):
        raise ValueError('Number of rows must either match or at least one must be one-dimensional.')

    if N_sigA == 1 and N_sigB != 1:
        A = _np.broadcast_to(A, (N_sigB, L_sigA))
    elif N_sigA != 1 and N_sigB == 1:
        B = _np.broadcast_to(B, (N_sigA, L_sigB))

    out = []

    for IDX, cur_row in enumerate(A):
        if FFT:
            out.append(fftconvolve(cur_row, B[IDX]))
        else:
            out.append(_np.convolve(cur_row, B[IDX]))

    return _np.array(out)