Python scipy.signal.freqz() Examples

def lsf2mfbe(lsf, mel_filters):

    NFFT = 512
    M = get_filterbank(n_filters=mel_filters, NFFT=NFFT, normalize=False, htk=True)

    mfbe = np.zeros(( len(lsf), mel_filters), dtype=np.float64)
    spec = np.zeros((len(lsf), NFFT/2+1,), dtype=np.float64)
   
    x = np.zeros((NFFT,), dtype=np.float64)
    x[0] = 1.0
    b = np.ones((1,), dtype=np.float64)
 
    for i, lsf_vec in enumerate(lsf):     
        #convert lsf to filter polynomial
        a_poly = lsf2poly(lsf_vec)
        # compute power spectrum
        w, H = freqz(b=1.0, a=a_poly, worN=NFFT, whole=True)
        spec_vec = np.abs(H[:(NFFT/2+1)])
        #spec_vec = np.square(spec_vec)
        # apply filterbank matrix
        mfbe[i,:] = np.log10( np.dot(M,spec_vec) )
        spec[i,:] = spec_vec
         
    return mfbe, spec 

def test_plot(self):

        def plot(w, h):
            assert_array_almost_equal(w, np.pi * np.arange(8.0) / 8)
            assert_array_almost_equal(h, np.ones(8))

        assert_raises(ZeroDivisionError,
                      freqz, [1.0], worN=8, plot=lambda w, h: 1 / 0)
        freqz([1.0], worN=8, plot=plot) 

def freqz_(sys, w, dt=8e-9):
    """
    This function computes the frequency response of a zpk system at an
    array of frequencies.

    It loosely mimicks 'scipy.signal.frequresp'.

    Parameters
    ----------
    system: (zeros, poles, k)
        zeros and poles both in rad/s, k is the actual coefficient, not DC gain
    w: np.array
        frequencies in rad/s
    dt: sampling time

    Returns
    -------
    np.array(..., dtype=np.complex) with the response
    """
    z, p, k = sys
    b, a = sig.zpk2tf(z, p, k)
    _, h = sig.freqz(b, a, worN=w*dt)
    return h 

def test_ticket1441(self):
        """Regression test for ticket 1441."""
        # Because freqz previously used arange instead of linspace,
        # when N was large, it would return one more point than
        # requested.
        N = 100000
        w, h = freqz([1.0], worN=N)
        assert_equal(w.shape, (N,)) 

def spec2lsf(spec, lsf_order=30):

    NFFT = 2*(spec.shape[0]-1)
    n_frames = spec.shape[1]

    p = lsf_order

    lsf = np.zeros(( n_frames, lsf_order), dtype=np.float64)
    spec_rec = np.zeros(spec.shape)

    for i, spec_vec in enumerate(spec.T):
    
        # floor reconstructed spectrum
        spec_vec = np.maximum(spec_vec, 1e-9)
 
        # squared magnitude 2-sided spectrum
        twoside = np.r_[spec_vec, np.flipud(spec_vec[1:-1])]
        twoside = np.square(twoside) 
        r = np.fft.ifft(twoside)
        r = r.real
  
        # levinson-durbin
        a = LA.solve_toeplitz(r[0:p],r[1:p+1])
        a = np.r_[1.0, -1.0*a]
   
        lsf[i,:] = poly2lsf(a)
   
        # reconstructed all-pole spectrum
        w, H = freqz(b=1.0, a=a, worN=NFFT, whole=True)
        spec_rec[:,i] = np.abs(H[:(NFFT/2+1)])
            
    return lsf, spec_rec 

def mfbe2lsf(mfbe, lsf_order):

    NFFT = 512
    M = get_filterbank(n_filters=mfbe.shape[1], NFFT=NFFT, normalize=False, htk=True)

    M_inv = pinv(M)
    p = lsf_order

    lsf = np.zeros(( len(mfbe), lsf_order), dtype=np.float64)
    spec = np.zeros((len(mfbe), NFFT/2+1), dtype=np.float64)

    for i, mfbe_vec in enumerate(mfbe):
    
        # invert mel filterbank
        spec_vec = np.dot(M_inv, np.power(10, mfbe_vec))

        # floor reconstructed spectrum
        spec_vec = np.maximum(spec_vec, 1e-9)
 
        # squared magnitude 2-sided spectrum
        twoside = np.r_[spec_vec, np.flipud(spec_vec[1:-1])]
        twoside = np.square(twoside) 
        r = np.fft.ifft(twoside)
        r = r.real

        # reference from talkbox
        # a,_,_ = TB.levinson(r, order=p)
  
        # levinson-durbin
        a = LA.solve_toeplitz(r[0:p],r[1:p+1])
        a = np.r_[1.0, -1.0*a]
   
        lsf[i,:] = poly2lsf(a)
   
        # reconstructed all-pole spectrum
        w, H = freqz(b=1.0, a=a, worN=NFFT, whole=True)
        spec[i,:] = np.abs(H[:(NFFT/2+1)])
            
    return lsf, spec 

def arma_periodogram(ar, ma, worN=None, whole=0):
    '''periodogram for ARMA process given by lag-polynomials ar and ma

    Parameters
    ----------
    ar : array_like
        autoregressive lag-polynomial with leading 1 and lhs sign
    ma : array_like
        moving average lag-polynomial with leading 1
    worN : {None, int}, optional
        option for scipy.signal.freqz   (read "w or N")
        If None, then compute at 512 frequencies around the unit circle.
        If a single integer, the compute at that many frequencies.
        Otherwise, compute the response at frequencies given in worN
    whole : {0,1}, optional
        options for scipy.signal.freqz
        Normally, frequencies are computed from 0 to pi (upper-half of
        unit-circle.  If whole is non-zero compute frequencies from 0 to 2*pi.

    Returns
    -------
    w : array
        frequencies
    sd : array
        periodogram, spectral density

    Notes
    -----
    Normalization ?

    This uses signal.freqz, which does not use fft. There is a fft version
    somewhere.

    '''
    w, h = signal.freqz(ma, ar, worN=worN, whole=whole)
    sd = np.abs(h)**2/np.sqrt(2*np.pi)
    if np.sum(np.isnan(h)) > 0:
        # this happens with unit root or seasonal unit root'
        print('Warning: nan in frequency response h, maybe a unit root')
    return w, sd 

def MOEar(self, correctionType = 'ELC'):
        """ Method to approximate middle-outer ear transfer function for linearly scaled
            frequency representations, using an FIR approximation of order 600 taps.
            As appears in :
            - A. Härmä, and K. Palomäki, ''HUTear – a free Matlab toolbox for modeling of human hearing'',
            in Proceedings of the Matlab DSP Conference, pp 96-99, Espoo, Finland 1999.
        Arguments          :
            correctionType : (string)     String which specifies the type of correction :
                                          'ELC' - Equal Loudness Curves at 60 dB (default)
                                          'MAP' - Minimum Audible Pressure at ear canal
                                          'MAF' - Minimum Audible Field
        Returns            :
            LTq            : (ndarray)    1D Array containing the transfer function, without the DC sub-band.
        """
        # Parameters
        firOrd = self.nfft
        Cr, fr, Crdb = self.OutMidCorrection(correctionType, firOrd, self.fs)
        Cr[self.nfft - 1] = 0.

        # FIR Design
        A = firwin2(firOrd, fr, Cr, nyq = self.fs/2)
        B = 1
        _, LTq = freqz(A, B, firOrd, self.fs)

        LTq = 20. * np.log10(np.abs(LTq))
        LTq -= max(LTq)
        return LTq[:self.nfft/2 + 1] 

def freqz_cas(sos,w):
    """
    Cascade frequency response
    
    Mark Wickert October 2016
    """
    Ns,Mcol = sos.shape
    w,Hcas = signal.freqz(sos[0,:3],sos[0,3:],w)
    for k in range(1,Ns):
        w,Htemp = signal.freqz(sos[k,:3],sos[k,3:],w)
        Hcas *= Htemp
    return w, Hcas 

def filter2spectrum(self, param):
        return np.abs(freqz([1], param)[1]) 

def test_absorption_filter(self):
        b = uwa.absorption_filter(200000)
        w, h = sp.freqz(b, 1, 4)
        h = utils.mag2db(np.abs(h))
        self.assertEqual(list(np.round(h)), [0.0, -3.0, -11.0, -22.0]) 

def test_basic_whole(self):
        w, h = freqz([1.0], worN=8, whole=True)
        assert_array_almost_equal(w, 2 * np.pi * np.arange(8.0) / 8)
        assert_array_almost_equal(h, np.ones(8)) 

def test_basic(self):
        w, h = freqz([1.0], worN=8)
        assert_array_almost_equal(w, np.pi * np.arange(8.0) / 8)
        assert_array_almost_equal(h, np.ones(8)) 

def ten_band_eq_resp(GdB,Q=3.5):
    """
    Create a frequency response magnitude plot in dB of a ten band equalizer
    using a semilogplot (semilogx()) type plot
    
    
    Parameters
    ----------
    GdB : Gain vector for 10 peaking filters [G0,...,G9]
    Q : Quality factor for each peaking filter (default 3.5)
    
    Returns
    -------
    Nothing : two plots are created
    
    Examples
    --------
    >>> import matplotlib.pyplot as plt
    >>> from sk_dsp_comm import sigsys as ss
    >>> ss.ten_band_eq_resp([0,10.0,0,0,-1,0,5,0,-4,0])
    >>> plt.show()
    """
    fs = 44100.0 # Hz
    NB = len(GdB)
    if not NB == 10:
        raise ValueError("GdB length not equal to ten")
    Fc = 31.25*2**np.arange(NB)
    B = np.zeros((NB,3));
    A = np.zeros((NB,3));
    
    # Create matrix of cascade coefficients
    for k in range(NB):
        b,a = peaking(GdB[k],Fc[k],Q,fs)
        B[k,:] = b
        A[k,:] = a
    # Create the cascade frequency response
    F = np.logspace(1,np.log10(20e3),1000)
    H = np.ones(len(F))*np.complex(1.0,0.0)
    for k in range(NB):
       w,Htemp = signal.freqz(B[k,:],A[k,:],2*np.pi*F/fs)
       H *= Htemp
    plt.figure(figsize=(6,4))
    plt.subplot(211)
    plt.semilogx(F,20*np.log10(abs(H)))
    plt.axis([10, fs/2, -12, 12])
    plt.grid()
    plt.title('Ten-Band Equalizer Frequency Response')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain (dB)')
    plt.subplot(212)
    plt.stem(np.arange(NB),GdB,'b','bs')
    #plt.bar(np.arange(NB)-.1,GdB,0.2)
    plt.axis([0, NB-1, -12, 12])
    plt.xlabel('Equalizer Band Number')
    plt.ylabel('Gain Set (dB)')
    plt.grid() 

def peaking(GdB, fc, Q=3.5, fs=44100.):
    """
    A second-order peaking filter having GdB gain at fc and approximately
    and 0 dB otherwise.
    
    The filter coefficients returns correspond to a biquadratic system function
    containing five parameters.
    
    Parameters
    ----------
    GdB : Lowpass gain in dB
    fc : Center frequency in Hz
    Q : Filter Q which is inversely proportional to bandwidth
    fs : Sampling frquency in Hz
    
    Returns
    -------
    b : ndarray containing the numerator filter coefficients
    a : ndarray containing the denominator filter coefficients
    
    Examples
    --------
    >>> import matplotlib.pyplot as plt
    >>> import numpy as np
    >>> from sk_dsp_comm.sigsys import peaking
    >>> from scipy import signal
    >>> b,a = peaking(2.0,500)
    >>> f = np.logspace(1,5,400)
    >>> w,H = signal.freqz(b,a,2*np.pi*f/44100)
    >>> plt.semilogx(f,20*np.log10(abs(H)))
    >>> plt.ylabel("Power Spectral Density (dB)")
    >>> plt.xlabel("Frequency (Hz)")
    >>> plt.show()

    >>> b,a = peaking(-5.0,500,4)
    >>> w,H = signal.freqz(b,a,2*np.pi*f/44100)
    >>> plt.semilogx(f,20*np.log10(abs(H)))
    >>> plt.ylabel("Power Spectral Density (dB)")
    >>> plt.xlabel("Frequency (Hz)")
    """
    mu = 10**(GdB/20.)
    kq = 4/(1 + mu)*np.tan(2*np.pi*fc/fs/(2*Q))
    Cpk = (1 + kq *mu)/(1 + kq)
    b1 = -2*np.cos(2*np.pi*fc/fs)/(1 + kq*mu)
    b2 = (1 - kq*mu)/(1 + kq*mu)
    a1 = -2*np.cos(2*np.pi*fc/fs)/(1 + kq)
    a2 = (1 - kq)/(1 + kq)
    b = Cpk*np.array([1, b1, b2])
    a = np.array([1, a1, a2])
    return b,a 

def am_rx_BPF(N_order = 7, ripple_dB = 1, B = 10e3, fs = 192e3):
    """
    Bandpass filter design for the AM receiver Case Study of Chapter 17.

    Design a 7th-order Chebyshev type 1 bandpass filter to remove/reduce
    adjacent channel intereference at the envelope detector input.
    
    Parameters
    ----------
    N_order : the filter order (default = 7)
    ripple_dB : the passband ripple in dB (default = 1)
    B : the RF bandwidth (default = 10e3)
    fs : the sampling frequency 

    Returns
    -------
    b_bpf : ndarray of the numerator filter coefficients
    a_bpf : ndarray of the denominator filter coefficients

    Examples
    --------
    >>> from scipy import signal
    >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> import sk_dsp_comm.sigsys as ss
    >>> # Use the default values
    >>> b_bpf,a_bpf = ss.am_rx_BPF()

    Pole-zero plot of the filter.

    >>> ss.zplane(b_bpf,a_bpf)
    >>> plt.show()

    Plot of the frequency response.

    >>> f = np.arange(0,192/2.,.1)
    >>> w, Hbpf = signal.freqz(b_bpf,a_bpf,2*np.pi*f/192)
    >>> plt.plot(f*10,20*np.log10(abs(Hbpf)))
    >>> plt.axis([0,1920/2.,-80,10])
    >>> plt.ylabel("Power Spectral Density (dB)")
    >>> plt.xlabel("Frequency (kHz)")
    >>> plt.show()
    """
    b_bpf,a_bpf = signal.cheby1(N_order,ripple_dB,2*np.array([75e3-B/2.,75e3+B/2.])/fs,'bandpass')
    return b_bpf,a_bpf 

def freqz(b, a=1, fs=2.0, worN=None, whole=False, degrees=True, style='solid', thickness=1, title=None, xlabel='Frequency (Hz)', xlim=None, ylim=None, width=None, height=None, hold=False, interactive=None):
    """Plot frequency response of a filter.

    This is a convenience function to plot frequency response, and internally uses
    :func:`scipy.signal.freqz` to estimate the response. For further details, see the
    documentation for :func:`scipy.signal.freqz`.

    :param b: numerator of a linear filter
    :param a: denominator of a linear filter
    :param fs: sampling rate in Hz (optional, normalized frequency if not specified)
    :param worN: see :func:`scipy.signal.freqz`
    :param whole: see :func:`scipy.signal.freqz`
    :param degrees: True to display phase in degrees, False for radians
    :param style: line style ('solid', 'dashed', 'dotted', 'dotdash', 'dashdot')
    :param thickness: line width in pixels
    :param title: figure title
    :param xlabel: x-axis label
    :param ylabel1: y-axis label for magnitude
    :param ylabel2: y-axis label for phase
    :param xlim: x-axis limits (min, max)
    :param ylim: y-axis limits (min, max)
    :param width: figure width in pixels
    :param height: figure height in pixels
    :param interactive: enable interactive tools (pan, zoom, etc) for plot
    :param hold: if set to True, output is not plotted immediately, but combined with the next plot

    >>> import arlpy
    >>> arlpy.plot.freqz([1,1,1,1,1], fs=120000);
    """
    w, h = _sig.freqz(b, a, worN, whole)
    Hxx = 20*_np.log10(abs(h)+_np.finfo(float).eps)
    f = w*fs/(2*_np.pi)
    if xlim is None:
        xlim = (0, fs/2)
    if ylim is None:
        ylim = (_np.max(Hxx)-50, _np.max(Hxx)+10)
    figure(title=title, xlabel=xlabel, ylabel='Amplitude (dB)', xlim=xlim, ylim=ylim, width=width, height=height, interactive=interactive)
    _hold_enable(True)
    plot(f, Hxx, color=color(0), style=style, thickness=thickness, legend='Magnitude')
    fig = gcf()
    units = 180/_np.pi if degrees else 1
    fig.extra_y_ranges = {'phase': _bmodels.Range1d(start=-_np.pi*units, end=_np.pi*units)}
    fig.add_layout(_bmodels.LinearAxis(y_range_name='phase', axis_label='Phase (degrees)' if degrees else 'Phase (radians)'), 'right')
    phase = _np.angle(h)*units
    fig.line(f, phase, line_color=color(1), line_dash=style, line_width=thickness, legend_label='Phase', y_range_name='phase')
    _hold_enable(hold) 

def highpass(signal, Fs, fc=None, plot=False):
    """ Filter out the really low frequencies, default is below 50Hz """

    if fc is None:
        fc = constants.get("fc_hp")

    # have some predefined parameters
    rp = 5  # minimum ripple in dB in pass-band
    rs = 60  # minimum attenuation in dB in stop-band
    n = 4  # order of the filter
    type = "butter"

    # normalized cut-off frequency
    wc = 2.0 * fc / Fs

    # design the filter
    from scipy.signal import iirfilter, lfilter, freqz

    b, a = iirfilter(n, Wn=wc, rp=rp, rs=rs, btype="highpass", ftype=type)

    # plot frequency response of filter if requested
    if plot:
        try:
            import matplotlib.pyplot as plt
        except ImportError:
            import warnings

            warnings.warn("Matplotlib is required for plotting")
            return

        w, h = freqz(b, a)

        plt.figure()
        plt.title("Digital filter frequency response")
        plt.plot(w, 20 * np.log10(np.abs(h)))
        plt.title("Digital filter frequency response")
        plt.ylabel("Amplitude Response [dB]")
        plt.xlabel("Frequency (rad/sample)")
        plt.grid()

    # apply the filter
    signal = lfilter(b, a, signal.copy())

    return signal 

def _test_phaseshift(self, method, zero_phase):
        rate = 120
        rates_to = [15, 20, 30, 40]  # q = 8, 6, 4, 3

        t_tot = int(100)  # Need to let antialiasing filters settle
        t = np.arange(rate*t_tot+1) / float(rate)

        # Sinusoids at 0.8*nyquist, windowed to avoid edge artifacts
        freqs = np.array(rates_to) * 0.8 / 2
        d = (np.exp(1j * 2 * np.pi * freqs[:, np.newaxis] * t)
             * signal.windows.tukey(t.size, 0.1))

        for rate_to in rates_to:
            q = rate // rate_to
            t_to = np.arange(rate_to*t_tot+1) / float(rate_to)
            d_tos = (np.exp(1j * 2 * np.pi * freqs[:, np.newaxis] * t_to)
                     * signal.windows.tukey(t_to.size, 0.1))

            # Set up downsampling filters, match v0.17 defaults
            if method == 'fir':
                n = 30
                system = signal.dlti(signal.firwin(n + 1, 1. / q,
                                                   window='hamming'), 1.)
            elif method == 'iir':
                n = 8
                wc = 0.8*np.pi/q
                system = signal.dlti(*signal.cheby1(n, 0.05, wc/np.pi))

            # Calculate expected phase response, as unit complex vector
            if zero_phase is False:
                _, h_resps = signal.freqz(system.num, system.den,
                                          freqs/rate*2*np.pi)
                h_resps /= np.abs(h_resps)
            else:
                h_resps = np.ones_like(freqs)

            y_resamps = signal.decimate(d.real, q, n, ftype=system,
                                        zero_phase=zero_phase)

            # Get phase from complex inner product, like CSD
            h_resamps = np.sum(d_tos.conj() * y_resamps, axis=-1)
            h_resamps /= np.abs(h_resamps)
            subnyq = freqs < 0.5*rate_to

            # Complex vectors should be aligned, only compare below nyquist
            assert_allclose(np.angle(h_resps.conj()*h_resamps)[subnyq], 0,
                            atol=1e-3, rtol=1e-3) 

def compute_frequency_response(filter_coefs, a_vals, fs):
    """Compute the frequency response of a filter.

    Parameters
    ----------
    filter_coefs : 1d or 2d array
        If 1d, interpreted as the B-value filter coefficients.
        If 2d, interpreted as the second-order (sos) filter coefficients.
    a_vals : 1d array or None
        The A-value filter coefficients for a filter.
        If second-order filter coefficients are provided in `filter_coefs`, must be None.
    fs : float
        Sampling rate, in Hz.

    Returns
    -------
    f_db : 1d array
        Frequency vector corresponding to attenuation decibels, in Hz.
    db : 1d array
        Degree of attenuation for each frequency specified in `f_db`, in dB.

    Examples
    --------
    Compute the frequency response for an FIR filter:

    >>> from neurodsp.filt.fir import design_fir_filter
    >>> filter_coefs = design_fir_filter(fs=500, pass_type='bandpass', f_range=(8, 12))
    >>> f_db, db = compute_frequency_response(filter_coefs, 1, fs=500)

    Compute the frequency response for an IIR filter, which uses SOS coefficients:

    >>> from neurodsp.filt.iir import design_iir_filter
    >>> sos_coefs = design_iir_filter(fs=500, pass_type='bandpass',
    ...                               f_range=(8, 12), butterworth_order=3)
    >>> f_db, db = compute_frequency_response(sos_coefs, None, fs=500)
    """

    if filter_coefs.ndim == 1 and a_vals is not None:
        # Compute response for B & A value filter coefficient inputs
        w_vals, h_vals = freqz(filter_coefs, a_vals, worN=int(fs * 2))
    elif filter_coefs.ndim == 2 and a_vals is None:
        # Compute response for sos filter coefficient inputs
        w_vals, h_vals = sosfreqz(filter_coefs, worN=int(fs * 2))
    else:
        raise ValueError("The organization of the filter coefficient inputs is not understood.")

    f_db = w_vals * fs / (2. * np.pi)
    db = 20 * np.log10(abs(h_vals))

    return f_db, db 

def tf_coefficients(self, frequencies=None, coefficients=None,
                                delay=False):
        """
        computes implemented transfer function - assuming no delay and
        infinite precision (actually floating-point precision)
        Returns the discrete transfer function realized by coefficients at
        frequencies.

        Parameters
        ----------
        coefficients: np.array
            coefficients as returned from iir module

        frequencies: np.array
            frequencies to compute the transfer function for

        dt: float
            discrete sampling time (seconds)

        zoh: bool
            If true, zero-order hold implementation is assumed. Otherwise,
            the delay is expected to depend on the index of biquad.

        Returns
        -------
        np.array(..., dtype=np.complex)
        """
        if frequencies is None:
            frequencies = self.frequencies
        frequencies = np.asarray(frequencies, dtype=np.float64)
        if coefficients is None:
            fcoefficients = self.coefficients
        else:
            fcoefficients = coefficients
        # discrete frequency
        w = frequencies * 2 * np.pi * self.dt * self.loops
        # the higher stages have progressively more delay to the output
        if delay:
            delay_per_cycle = np.exp(-1j * self.dt * frequencies * 2 * np.pi)
        h = np.zeros(len(w), dtype=np.complex128)
        for i in range(len(fcoefficients)):
            sos = np.asarray(fcoefficients[i], dtype=np.float64)
            # later we can use sig.sosfreqz (very recent function, dont want
            #  to update scipy now)
            ww, hh = sig.freqz(sos[:3], sos[3:], worN=np.asarray(w,
                                                           dtype=np.float64))
            if delay:
                hh *= delay_per_cycle ** i
            h += hh
        return h 

def plot_filter_response(
    frequency,
    sample_rate,
    filter_type,
    filter_order=2,
    show_grid=True,
    fig_size=(10, 5),
):
    """Plot filter frequency response.

    Generate a plot showing the frequency response curve of a filter with the
    specified parameters.

    Parameters
    ----------
    frequency : int or tuple of ints
        The cutoff frequency for the filter. If `filter_type` is set as
        'bandpass' then this needs to be a tuple of integers representing
        the lower and upper bound frequencies. For example, for a bandpass
        filter with range of 2Hz and 10Hz, you would pass in the tuple (2, 10).
        For filter types with a single cutoff frequency then a single integer
        should be used.
    sample_rate : float
        Sampling rate of the signal.
    filter_type : {'lowpass', 'highpass', 'bandpass', 'bandstop'}
        Type of filter to build.
    filter_order: int, optional
        Order of the filter.
    show_grid : bool, optional
        Toggle to show grid lines on the plot.
    fig_size : tuple, optional
        Tuple containing the width and height of the resulting figure.
    """

    b, a = sensormotion.signal.build_filter(
        frequency, sample_rate, filter_type, filter_order
    )

    # Plot the frequency response
    w, h = freqz(b, a, worN=8000)
    f, axarr = plt.subplots(figsize=fig_size)
    axarr.plot(0.5 * sample_rate * w / np.pi, np.abs(h), "b")
    axarr.set_xlim(0, 0.5 * sample_rate)
    axarr.set_xlabel("Frequency (Hz)")
    axarr.grid(show_grid)

    # Add lines and markers at the cutoff frequency/frequencies
    if filter_type == "bandpass":
        for i in range(len(frequency)):
            axarr.axvline(frequency[i], color="k")
            axarr.plot(frequency[i], 0.5 * np.sqrt(2), "ko")
    else:
        axarr.axvline(frequency, color="k")
        axarr.plot(frequency, 0.5 * np.sqrt(2), "ko")

    plt.suptitle("Filter Frequency Response", size=16)
    plt.show() 

def arma_periodogram(ar, ma, worN=None, whole=0):
    """
    Periodogram for ARMA process given by lag-polynomials ar and ma

    Parameters
    ----------
    ar : array_like
        autoregressive lag-polynomial with leading 1 and lhs sign
    ma : array_like
        moving average lag-polynomial with leading 1
    worN : {None, int}, optional
        option for scipy.signal.freqz (read "w or N")
        If None, then compute at 512 frequencies around the unit circle.
        If a single integer, the compute at that many frequencies.
        Otherwise, compute the response at frequencies given in worN
    whole : {0,1}, optional
        options for scipy.signal.freqz
        Normally, frequencies are computed from 0 to pi (upper-half of
        unit-circle.  If whole is non-zero compute frequencies from 0 to 2*pi.

    Returns
    -------
    w : array
        frequencies
    sd : array
        periodogram, spectral density

    Notes
    -----
    Normalization ?

    This uses signal.freqz, which does not use fft. There is a fft version
    somewhere.
    """
    w, h = signal.freqz(ma, ar, worN=worN, whole=whole)
    sd = np.abs(h) ** 2 / np.sqrt(2 * np.pi)
    if np.any(np.isnan(h)):
        # this happens with unit root or seasonal unit root'
        import warnings
        warnings.warn('Warning: nan in frequency response h, maybe a unit '
                      'root', RuntimeWarning)
    return w, sd