Python scipy.signal.firwin() Examples

def low_cut_filter(x, fs, cutoff=70):
    """FUNCTION TO APPLY LOW CUT FILTER

    Args:
        x (ndarray): Waveform sequence
        fs (int): Sampling frequency
        cutoff (float): Cutoff frequency of low cut filter

    Return:
        (ndarray): Low cut filtered waveform sequence
    """

    nyquist = fs // 2
    norm_cutoff = cutoff / nyquist

    # low cut filter
    fil = firwin(255, norm_cutoff, pass_zero=False)
    lcf_x = lfilter(fil, 1, x)

    return lcf_x 

def mono_FM(x,fs=2.4e6,file_name='test.wav'):
    """
    Decimate complex baseband input by 10
    Design 1st decimation lowpass filter (f_c = 200 KHz)
    """
    b = signal.firwin(64,2*200e3/float(fs))
    # Filter and decimate (should be polyphase)
    y = signal.lfilter(b,1,x)
    z = ss.downsample(y,10)
    # Apply complex baseband discriminator
    z_bb = discrim(z)
    # Design 2nd decimation lowpass filter (fc = 12 KHz)
    bb = signal.firwin(64,2*12e3/(float(fs)/10))
    # Filter and decimate
    zz_bb = signal.lfilter(bb,1,z_bb)
    # Decimate by 5
    z_out = ss.downsample(zz_bb,5)
    # Save to wave file
    ss.to_wav(file_name, 48000, z_out/2)
    print('Done!')
    return z_bb, z_out 

def design(self, fs, fc, bandwidth, ripple_db=60.0):
        """
        Designs a FIR filter that is a low-pass filter.
        fs : sampling frequency (Hz)
        fc : cut-off frequency (Hz)
        bandwidth : transition bandwidth (Hz)s
        """
        # Compute the order and Kaiser parameter for the FIR filter.
        N, beta = signal.kaiserord(ripple_db, bandwidth / fs * 2)

        # Use firwin with a Kaiser window to create a lowpass FIR filter.
        fir = signal.firwin(N, fc / fs * 2, window=('kaiser', beta))

        # the filter must be symmetric, in order to be zero-phase
        assert np.all(np.abs(fir - fir[::-1]) < 1e-15)

        self.fir = fir / np.sum(fir)
        self.fs = fs
        return self 

def decimate(x,M,fs=2.4e6):
    b = signal.firwin(64,2*200e3/float(fs))
    # Filter and decimate (should be polyphase)
    y = signal.lfilter(b,1,x)
    z = ss.downsample(y,M)
    return z 

def __init__(self, id_mask='F1804', ntaps=4, debug=False):
        """Initialize the WeightSensor.

        Parameters
        ----------
        id_mask : str
            A template for the first n digits of the device IDs for valid load cells.
        ntaps : int
            Maximum number of samples to perform filtering over.
        debug : bool
            If True, have sensor seem to work normally but just return zeros.
        """
        self._id_mask = id_mask
        self._weight_buffers = []
        self._ntaps = ntaps
        self._debug = debug
        self._filter_coeffs = signal.firwin(ntaps, 0.1)
        self._running = False 

def __init__(self, fs, parameters):

        fs_min = 1000.0
        if (fs > fs_min):
            dec_factor = parameters['dec_factor']
        else:
            dec_factor = 1

        filter_order = parameters['bp_forder']
        f_hp = parameters['bp_low']
        f_lp = parameters['bp_high']

        f1 = f_hp/(fs/2)
        f2 = f_lp/(fs/2)

        self.b = firwin(filter_order+1, [f1, f2], pass_zero=False)
        self.a = 1
        self.dec_factor = dec_factor 

def apply(self, data):
        nyq = self.f / 2.0
        cutoff = min(self.f, nyq-1)
        h = signal.firwin(numtaps=101, cutoff=cutoff, nyq=nyq)

        # data[i][ch][dim0]
        for i in range(len(data)):
            data_point = data[i]
            for j in range(len(data_point)):
                data_point[j] = signal.lfilter(h, 1.0, data_point[j])

        return data 

def smooth_ppg_signal(signal, sample_rate=PPG_SAMPLE_RATE, numtaps=PPG_FIR_FILTER_TAP_NUM, cutoff=PPG_FILTER_CUTOFF):
    if numtaps % 2 == 0:
        numtaps += 1
    return convolve(signal, firwin(numtaps, [x*2/sample_rate for x in cutoff], pass_zero=False), mode='valid').tolist() 

def triangular(ar, header, freqmin, freqmax, zerophase=True, verbose=False):
    """
    Vertical triangular FIR bandpass. This filter is designed to closely emulate that of RADAN.

    Filter design is implemented by :py:func:`scipy.signal.firwin` with :code:`numtaps=25` and implemented with :py:func:`scipy.signal.lfilter`.

    .. note:: This function is not compatible with scipy versions prior to 1.3.0.

    :param np.ndarray ar: The radar array
    :param dict header: The file header dictionary
    :param int freqmin: The lower corner of the bandpass
    :param int freqmax: The upper corner of the bandpass
    :param bool zerophase: Whether to run the filter forwards and backwards in order to counteract the phase shift
    :param bool verbose: Verbose, defaults to False
    :rtype: :py:class:`numpy.ndarray`
    """
    if verbose:
        fx.printmsg('vertical triangular FIR bandpass filter')
    #samp_freq = 1 / ((header['rhf_depth'] * 2) / header['cr'] / header['rh_nsamp'])
    samp_freq = header['samp_freq']
    freqmin = freqmin * 10 ** 6
    freqmax = freqmax * 10 ** 6
    
    numtaps = 25

    if verbose:
        fx.printmsg('sampling frequency:       %.2E Hz' % samp_freq)
        fx.printmsg('minimum filter frequency: %.2E Hz' % freqmin)
        fx.printmsg('maximum filter frequency: %.2E Hz' % freqmax)
        fx.printmsg('numtaps: %s, zerophase: %s' % (numtaps, zerophase))

    filt = firwin(numtaps=numtaps, cutoff=[freqmin, freqmax], window='triangle', pass_zero='bandpass', fs=samp_freq)

    far = lfilter(filt, 1.0, ar, axis=0).copy()
    if zerophase:
        far = lfilter(filt, 1.0, far[::-1], axis=0)[::-1]

    return far 

def getFilterParameters(self, fc_lp, fc_hp, n=511):
        # get lowpass and highpass filters
        n = 511
        fc = 1.5e3
        hp = firwin(n, fc_lp, fs=self.fs, pass_zero=False)
        lp = firwin(n, fc_hp, fs=self.fs, pass_zero=True)
        return hp, lp

    # define the transmit chirp signal sent over the speakers 

def polyphase_lowpass(arr, downsample=2, n_taps=50, filter_pad=1.1):
    filt = firwin(downsample * n_taps, 1 / (downsample * filter_pad))
    filtered = polyphase_single_filter(arr, downsample, filt)
    return filtered 

def polyphase_lowpass(arr, downsample=2, n_taps=50, filter_pad=1.1):
    filt = firwin(downsample * n_taps, 1 / (downsample * filter_pad))
    filtered = polyphase_single_filter(arr, downsample, filt)
    return filtered 

def low_cut_filter(x, fs, cutoff=70):
    nyquist = fs // 2
    norm_cutoff = cutoff / nyquist

    # low cut filter
    fil = firwin(255, norm_cutoff, pass_zero=False)
    lcf_x = lfilter(fil, 1, x)

    return lcf_x 

def _design(self):
        # Compute the order and Kaiser parameter for the FIR filter.
        N, beta = signal.kaiserord(self.ripple_db,
                                   self.bandwidth / self.fs * 2)

        # Use firwin with a Kaiser window to create a lowpass FIR filter.
        fir = signal.firwin(N, self.fc / self.fs * 2, window=('kaiser', beta))

        # the filter must be symmetric, in order to be zero-phase
        assert np.all(np.abs(fir - fir[::-1]) < 1e-15)

        self.fir = fir / np.sum(fir)
        return self 

def load_filterbank(bandwidth,fs, order = 4, max_freq = 40,ftype = 'butter'): 
	'''	Calculate Filters bank with Butterworth filter  

	Keyword arguments:
	bandwith -- numpy array containing bandwiths ex. [2,4,8,16,32]
	f_s -- sampling frequency

	Return:	numpy array containing filters coefficients dimesnions 'butter': [N_bands,order,6] 'fir': [N_bands,order]
	'''
	
	f_band_nom = load_bands(bandwidth,fs,max_freq) # get normalized bands 
	n_bands = f_band_nom.shape[0]
	
	if ftype == 'butter': 
		filter_bank = np.zeros((n_bands,order,6))
	elif ftype == 'fir':
		filter_bank = np.zeros((n_bands,order))



	for band_idx in range(n_bands):
		if ftype == 'butter': 
			filter_bank[band_idx] = butter(order, f_band_nom[band_idx], analog=False, btype='band', output='sos')
		elif ftype == 'fir':
			
			
			filter_bank[band_idx] = signal.firwin(order,f_band_nom[band_idx],pass_zero=False)
	return filter_bank 

def polyphase_lowpass(arr, downsample=2, n_taps=50, filter_pad=1.1):
    filt = firwin(downsample * n_taps, 1 / (downsample * filter_pad))
    filtered = polyphase_single_filter(arr, downsample, filt)
    return filtered 

def low_cut_filter(x, fs, cutoff=70):
    """Low cut filter

    Parameters
    ---------
    x : array, shape(`samples`)
        Waveform sequence
    fs: array, int
        Sampling frequency
    cutoff : float, optional
        Cutoff frequency of low cut filter
        Default set to 70 [Hz]

    Returns
    ---------
    lcf_x : array, shape(`samples`)
        Low cut filtered waveform sequence
    """

    nyquist = fs // 2
    norm_cutoff = cutoff / nyquist

    # low cut filter
    fil = firwin(255, norm_cutoff, pass_zero=False)
    lcf_x = lfilter(fil, 1, x)

    return lcf_x 

def stereo_FM(x,fs=2.4e6,file_name='test.wav'):
    """
    Stereo demod from complex baseband at sampling rate fs.
    Assume fs is 2400 ksps
    
    Mark Wickert July 2017
    """
    N1 = 10
    b = signal.firwin(64,2*200e3/float(fs))
    # Filter and decimate (should be polyphase)
    y = signal.lfilter(b,1,x)
    z = ss.downsample(y,N1)
    # Apply complex baseband discriminator
    z_bb = discrim(z)
    # Work with the (3) stereo multiplex signals:
    # Begin by designing a lowpass filter for L+R and DSP demoded (L-R)
    # (fc = 12 KHz)
    b12 = signal.firwin(128,2*12e3/(float(fs)/N1))
    # The L + R term is at baseband, we just lowpass filter to remove 
    # other terms above 12 kHz.
    y_lpr = signal.lfilter(b12,1,z_bb)
    b19 = signal.firwin(128,2*1e3*np.array([19-5,19+5])/(float(fs)/N1),
                        pass_zero=False);
    z_bb19 = signal.lfilter(b19,1,z_bb)
    # Lock PLL to 19 kHz pilot
    # A type 2 loop with bandwidth Bn = 10 Hz and damping zeta = 0.707 
    # The VCO quiescent frequency is set to 19000 Hz.
    theta, phi_error = pilot_PLL(z_bb19,19000,fs/N1,2,10,0.707)
    # Coherently demodulate the L - R subcarrier at 38 kHz.
    # theta is the PLL output phase at 19 kHz, so to double multiply 
    # by 2 and wrap with cos() or sin().
    # First bandpass filter
    b38 = signal.firwin(128,2*1e3*np.array([38-5,38+5])/(float(fs)/N1),
                        pass_zero=False);
    x_lmr = signal.lfilter(b38,1,z_bb)
    # Coherently demodulate using the PLL output phase
    x_lmr = 2*np.sqrt(2)*np.cos(2*theta)*x_lmr
    # Lowpass at 12 kHz to recover the desired DSB demod term
    y_lmr = signal.lfilter(b12,1,x_lmr)
    # Matrix the y_lmr and y_lpr for form right and left channels:
    y_left = y_lpr + y_lmr
    y_right = y_lpr - y_lmr
    
    # Decimate by N2 (nominally 5)
    N2 = 5
    fs2 = float(fs)/(N1*N2) # (nominally 48 ksps)
    y_left_DN2 = ss.downsample(y_left,N2)
    y_right_DN2 = ss.downsample(y_right,N2)
    # Deemphasize with 75 us time constant to 'undo' the preemphasis 
    # applied at the transmitter in broadcast FM.
    # A 1-pole digital lowpass works well here.
    a_de = np.exp(-2.1*1e3*2*np.pi/fs2)
    z_left = signal.lfilter([1-a_de],[1, -a_de],y_left_DN2)
    z_right = signal.lfilter([1-a_de],[1, -a_de],y_right_DN2)
    # Place left and righ channels as side-by-side columns in a 2D array
    z_out = np.hstack((np.array([z_left]).T,(np.array([z_right]).T)))
    
    ss.to_wav(file_name, 48000, z_out/2)
    print('Done!')
    #return z_bb, z_out
    return z_bb, theta, y_lpr, y_lmr, z_out 

def firwin_bpf(N_taps, f1, f2, fs = 1.0, pass_zero=False):
    """
    Design a windowed FIR bandpass filter in terms of passband
    critical frequencies f1 < f2 in Hz relative to sampling rate
    fs in Hz. The number of taps must be provided.

    Mark Wickert October 2016
    """
    return signal.firwin(N_taps,2*(f1,f2)/fs,pass_zero=pass_zero) 

def firwin_lpf(N_taps, fc, fs = 1.0):
    """
    Design a windowed FIR lowpass filter in terms of passband
    critical frequencies f1 < f2 in Hz relative to sampling rate
    fs in Hz. The number of taps must be provided.
    
    Mark Wickert October 2016
    """
    return signal.firwin(N_taps,2*fc/fs) 

def save_wav(wav, path, hparams):
	wav = wav / np.abs(wav).max() * 0.999
	f1 = 0.5 * 32767 / max(0.01, np.max(np.abs(wav)))
	f2 = np.sign(wav) * np.power(np.abs(wav), 0.95)
	wav = f1 * f2
	wav = signal.convolve(wav, signal.firwin(hparams.num_freq, [hparams.fmin, hparams.fmax], pass_zero=False, fs=hparams.sample_rate))
	#proposed by @dsmiller
	wavfile.write(path, hparams.sample_rate, wav.astype(np.int16)) 

def pb2bb(x, fs, fc, fd=None, flen=127, cutoff=None, axis=-1):
    """Convert passband signal to baseband.

    The baseband conversion uses a low-pass filter after downconversion, with a
    default cutoff frequency of `0.6*fd`, if `fd` is specified, or `1.1*fc` if `fd`
    is not specified. Alternatively, the user may specify the cutoff frequency
    explicitly.

    For communication applications, one may wish to use :func:`arlpy.comms.downconvert` instead,
    as that function supports matched filtering with a pulse shape rather than a generic
    low-pass filter.

    The convention used in that exp(2j*pi*fc*t) is a positive frequency carrier.

    :param x: passband signal
    :param fs: sampling rate of passband signal in Hz
    :param fc: carrier frequency in passband in Hz
    :param fd: sampling rate of baseband signal in Hz (``None`` => same as `fs`)
    :param flen: number of taps in the low-pass FIR filter
    :param cutoff: cutoff frequency in Hz (``None`` means auto-select)
    :param axis: axis of the signal, if multiple signals specified
    :returns: complex baseband signal, sampled at `fd`
    """
    if cutoff is None:
        cutoff = 0.6*fd if fd is not None else 1.1*_np.abs(fc)
    osc = _np.sqrt(2)*_np.exp(-2j*_np.pi*fc*time(x.shape[axis],fs))
    y = x * _utils.broadcastable_to(osc, x.shape, axis)
    hb = _sig.firwin(flen, cutoff=cutoff, nyq=fs/2.0)
    y = _sig.filtfilt(hb, 1, y, axis=axis)
    if fd is not None and fd != fs:
        y = _sig.resample_poly(y, 2*fd, fs, axis=axis)
        y = _np.apply_along_axis(lambda a: a[::2], axis, y)
    return y 

def test_poly_vs_filtfilt(self):
        # Check that up=1.0 gives same answer as filtfilt + slicing
        random_state = np.random.RandomState(17)
        try_types = (int, np.float32, np.complex64, float, complex)
        size = 10000
        down_factors = [2, 11, 79]

        for dtype in try_types:
            x = random_state.randn(size).astype(dtype)
            if dtype in (np.complex64, np.complex128):
                x += 1j * random_state.randn(size)

            # resample_poly assumes zeros outside of signl, whereas filtfilt
            # can only constant-pad. Make them equivalent:
            x[0] = 0
            x[-1] = 0

            for down in down_factors:
                h = signal.firwin(31, 1. / down, window='hamming')
                yf = filtfilt(h, 1.0, x, padtype='constant')[::down]

                # Need to pass convolved version of filter to resample_poly,
                # since filtfilt does forward and backward, but resample_poly
                # only goes forward
                hc = convolve(h, h[::-1])
                y = signal.resample_poly(x, 1, down, window=hc)
                assert_allclose(yf, y, atol=1e-7, rtol=1e-7) 

def test_vs_lfilter(self):
        # Check that up=1.0 gives same answer as lfilter + slicing
        random_state = np.random.RandomState(17)
        try_types = (int, np.float32, np.complex64, float, complex)
        size = 10000
        down_factors = [2, 11, 79]

        for dtype in try_types:
            x = random_state.randn(size).astype(dtype)
            if dtype in (np.complex64, np.complex128):
                x += 1j * random_state.randn(size)

            for down in down_factors:
                h = firwin(31, 1. / down, window='hamming')
                yl = lfilter(h, 1.0, x)[::down]
                y = upfirdn(h, x, up=1, down=down)
                assert_allclose(yl, y[:yl.size], atol=1e-7, rtol=1e-7) 

def polyphase_lowpass(arr, downsample=2, n_taps=50, filter_pad=1.1):
    filt = firwin(downsample * n_taps, 1 / (downsample * filter_pad))
    filtered = polyphase_single_filter(arr, downsample, filt)
    return filtered 

def high_pass_filter(data, cutoff, fs, n_taps=255):
    """Apply high-pass filter

    Parameters
    ----------
    data : array, shape (`T`, `dim`)
        Array of sequence.
    cutoff : int,
        Cutoff frequency
    fs : int,
        Sampling frequency
    n_taps : int, optional
        Tap number

    Returns
    -------
    modified data: array, shape (`T`, `dim`)
        Array of modified sequence.
    """
    if data.shape[0] < n_taps * 3:
        raise ValueError(
            'Length of data should be three times longer than n_taps.')

    fil = firwin(n_taps, cutoff, pass_zero=False, nyq=fs//2)
    modified_data = filtfilt(fil, 1, data, axis=0)
    return modified_data 

def low_pass_filter(data, cutoff, fs, n_taps=255):
    """Apply low-pass filter

    Parameters
    ----------
    data : array, shape (`T`, `dim`)
        Array of sequence.
    cutoff : int,
        Cutoff frequency
    fs : int,
        Sampling frequency
    n_taps : int, optional
        Tap number

    Returns
    -------
    modified data: array, shape (`T`, `dim`)
        Array of modified sequence.
    """
    if data.shape[0] < n_taps * 3:
        raise ValueError(
            'Length of data should be three times longer than n_taps.')

    fil = firwin(n_taps, cutoff, pass_zero=True, nyq=fs//2)
    modified_data = filtfilt(fil, 1, data, axis=0)
    return modified_data 

def _high_frequency_completion(self, x, transformed):
        """
        Please see Sect. 3.2 and 3.3 in the following paper to know why we complete the
        unvoiced synthesized voice of the original voice into high frequency range
        of F0 transformed voice.

        - K. Kobayashi et al., "F0 transformation techniques for statistical voice
        conversion with direct waveform modification with spectral differential,"
        Proc. IEEE SLT 2016, pp. 693-700. 2016.
        """
        # construct feature extractor and synthesis
        feat = FeatureExtractor(fs=self.fs)
        f0, spc, ap = feat.analyze(x)
        uf0 = np.zeros(len(f0))

        # synthesis
        synth = Synthesizer(fs=self.fs)
        unvoice_anasyn = synth.synthesis_spc(uf0, spc, ap)

        # HPF for synthesized speech
        fil = firwin(255, self.f0rate, pass_zero=False)
        HPFed_unvoice_anasyn = filtfilt(fil, 1, unvoice_anasyn)

        if len(HPFed_unvoice_anasyn) > len(transformed):
            return transformed + HPFed_unvoice_anasyn[:len(transformed)]
        else:
            transformed[:len(HPFed_unvoice_anasyn)] += HPFed_unvoice_anasyn
            return transformed 

def fast_xambg(refChannel, srvChannel, rangeBins, freqBins, shortFilt=True):
    ''' Fast Cross-Ambiguity Fuction (frequency domain method)
    
    Parameters:
        refChannel: reference channel data
        srvChannel: surveillance channel data
        rangeBins:  number of range bins to compute
        freqBins:   number of doppler bins to compute (should be power of 2)
        shortFilt:  (bool) chooses the type of decimation filter to use.
                    If True, uses an all-ones filter of length 1*(decimation factor)
                    If False, uses a flat-top window of length 10*(decimation factor)+1
    Returns:
        xambg: the cross-ambiguity surface. Dimensions are (nf, nlag+1, 1)
        third dimension added for easy stacking in dask

    '''
    if refChannel.shape != srvChannel.shape:
        raise ValueError('Input vectors must have the same length')

    # calculate decimation factor
    ndecim = int(refChannel.shape[0]/freqBins)

    # pre-allocate space for the result
    xambg = np.zeros((freqBins, rangeBins+1, 1), dtype=np.complex64)

    # complex conjugate of the second input vector
    srvChannelConj = np.conj(srvChannel)    

    if shortFilt:
        # precompute short FIR filter for decimation (all ones filter with length
        # equal to the decimation factor)
        dtaps = np.ones((ndecim + 1,))
    else:
        # precompute long FIR filter for decimation. (flat top filter of length
        # 10*decimation factor).  
        dtaps = signal.firwin(10*ndecim + 1, 1. / ndecim, window='flattop')

    dfilt = signal.dlti(dtaps, 1)

    # loop over range bins 
    for k, lag in enumerate(np.arange(-1*rangeBins, 1)):
        channelProduct = np.roll(srvChannelConj, lag)*refChannel
        #decimate the product of the reference channel and the delayed surveillance channel
        xambg[:,k,0] = signal.decimate(channelProduct, ndecim, ftype=dfilt)[0:freqBins]

    # take the FFT along the first axis (Doppler)
    # xambg = np.fft.fftshift(np.fft.fft(xambg, axis=0), axes=0)
    xambg = np.fft.fftshift(fft(xambg, axis=0), axes=0)
    return xambg 

def initialize_online_filter(fsample, highpass, lowpass, order, x, axis=-1):
    # boxcar, triang, blackman, hamming, hann, bartlett, flattop, parzen, bohman, blackmanharris, nuttall, barthann
    filtwin = 'nuttall'
    nyquist = fsample / 2.
    ndim = len(x.shape)
    axis = axis % ndim

    if highpass != None:
        highpass = highpass/nyquist
        if highpass < 0.001:
            print('Warning: highpass is too low, disabling')
            highpass = None
        elif highpass > 0.999:
            print('Warning: highpass is too high, disabling')
            highpass = None

    if lowpass != None:
        lowpass = lowpass/nyquist
        if lowpass < 0.001:
            print('Warning: lowpass is too low, disabling')
            lowpass = None
        elif lowpass > 0.999:
            print('Warning: lowpass is too low, disabling')
            lowpass = None

    if not(highpass is None) and not(lowpass is None) and highpass>=lowpass:
        # totally blocking all signal
        print('using NULL filter', [highpass, lowpass, order])
        b = np.zeros(order)
        a = np.ones(1)
    elif not(lowpass is None) and (highpass is None):
        print('using lowpass filter', [highpass, lowpass, order])
        b = firwin(order, cutoff = lowpass, window = filtwin, pass_zero = True)
        a = np.ones(1)
    elif not(highpass is None) and (lowpass is None):
        print('using highpass filter', [highpass, lowpass, order])
        b = firwin(order, cutoff = highpass, window = filtwin, pass_zero = False)
        a = np.ones(1)
    elif not(highpass is None) and not(lowpass is None):
        print('using bandpass filter', [highpass, lowpass, order])
        b = firwin(order, cutoff = [highpass, lowpass], window = filtwin, pass_zero = False)
        a = np.ones(1)
    else:
        # no filtering at all
        print('using IDENTITY filter', [highpass, lowpass, order])
        b = np.ones(1)
        a = np.ones(1)

    # initialize the state for the filtering based on the previous data
    if ndim == 1:
        zi = zi = lfiltic(b, a, x, x)
    elif ndim == 2:
        f = lambda x : lfiltic(b, a, x, x)
        zi = np.apply_along_axis(f, axis, x)

    return b, a, zi

####################################################################