def freqz(b, a=1, fs=1, xlim=None, N=1000, xlog=False):
    """Calculate the frequency response of a discrete time system over the
    range xlim, over a log or linear interval.

    Parameters
    ----------

    b : array-like
        Numerator coefficients of discrete time system
    a : array-like, optional
        Denominator coefficients of discrete time system
    fs : float, optional
        Sampling frequency; use to scale the output frequency array
    xlim : tuple of (x_min, x_max), optional
        Calculate the response from x_min to x_max. If omitted, the entire
        digital frequency axis is used
    N : int, optional
        The number of points to calculate the system response at
    xlog : bool, optional
        Calculate the frequency response at a log (True) or linearly spaced
        set of points"""
    # Squeeze arrays to deal with cont2discrete array issues
    b = np.squeeze(b)
    a = np.squeeze(a)
    if xlim is None:
        w, resp = signal.freqz(b, a)
        w = lin_or_logspace(w[w > 0][0], w[-1], N, True)
        _, resp = signal.freqz(b, a, w)
    else:
        w = 2 * np.pi * lin_or_logspace(xlim[0], xlim[1], N, xlog) / fs
        _, resp = signal.freqz(b, a, worN=w)

    freq = w * fs / (2 * np.pi)
    return freq, resp

        def _find_min_max(self, f_start, f_stop, unit = 'dB'):
            """
            Find minimum and maximum magnitude and the corresponding frequencies
            for the filter defined in the filter dict in a given frequency band
            [f_start, f_stop].
            """
            w = np.linspace(f_start, f_stop, params['N_FFT'])*2*np.pi
            [w, H] = sig.freqz(bb, aa, worN = w)

            # add antiCausals if we have them
            if (antiC):
               #
               # Evaluate transfer function of anticausal half on the same freq grid.
               #
               wa, ha = sig.freqz(bbA, aaA, worN = w)
               ha = ha.conjugate()
               #
               # Total transfer function is the product
               #
               H = H*ha

            f = w / (2.0 * pi) # frequency normalized to f_S
            H_abs = abs(H)
            H_max = max(H_abs)
            H_min = min(H_abs)
            F_max = f[np.argmax(H_abs)] # find the frequency where H_abs 
            F_min = f[np.argmin(H_abs)] # becomes max resp. min
            if unit == 'dB':
                H_max = 20*log10(H_max)
                H_min = 20*log10(H_min)
            return F_min, H_min, F_max, H_max

    def plot_response(self, ax):
        # Plot the designed filter response
        if self.n_section == 1:
            fw, fh = signal.freqz(self.b, self.a, worN=self.Nfft)
            ax.plot(fw, 20*log10(np.abs(fh)), linewidth=2, alpha=.75)
            fxw, fxh = signal.freqz(self.fxb, self.fxa, worN=self.Nfft)
            ax.plot(fxw, 20*log10(np.abs(fxh)), linestyle='--', linewidth=3,
                    alpha=.75)

        # plot the simulated response, if simulated data exists
        if self.xfavg is None or self.yfavg is None or self.pfavg is None:
            pass
        else:
            #  -- Fixed Point Sim --
            xa = 2*pi * np.arange(self.Nfft)/self.Nfft
            h = self.yfavg / self.xfavg
            ax.plot(xa, 20*log10(h), linewidth=4, alpha=.5)
            #  -- Floating Point Sim --
            hp = self.pfavg / self.xfavg
            ax.plot(xa, 20*log10(hp), color='k', linestyle=':',
                    linewidth=2, alpha=.75)

        ax.set_ylim(-80, 10)
        ax.set_xlim(0, np.pi)

        ax.set_ylabel('Magnitude dB');
        ax.set_xlabel('Frequency Normalized Radians')
        ax.legend(('Ideal', 'Quant. Coeff.',
                      'Fixed-P. Sim', 'Floating-P. Sim'))

def diffplot(freq, B, A, B2, A2):
	w, h = sps.freqz(B, A)
	w2, h2 = sps.freqz(B2, A2)

#	h = h - h2
	dabs = abs(h2) / abs(h)
	dphase = np.unwrap(np.angle(h2)) - np.unwrap(np.angle(h))  

	fig = plt.figure()
	plt.title('Difference between digital filter frequency responses')

	ax1 = fig.add_subplot(111)

	plt.plot(w * (freq/np.pi) / 2.0, 20 * np.log10(dabs), 'b')
	plt.ylabel('Amplitude [dB]', color='b')
	plt.xlabel('Frequency [rad/sample]')

	ax2 = ax1.twinx()
	angles = np.unwrap(np.angle(h))
	angles = dphase
	plt.plot(w * (freq/np.pi) / 2.0, angles, 'g')
	plt.ylabel('Angle (radians)', color='g')
	plt.grid()
	plt.axis('tight')
	plt.show()

def estimate(rg, filt, para, impulse, fs):
    aa = zeros(len(rg))
    Q = zeros(len(rg))
    maxh = zeros(len(rg))
    for i, freq in enumerate(rg):
        filt[para] = freq
        filt.compute(zeros(10000, dtype=float32))
        w, h = freqz(filt.compute(impulse), worN=15000)
        k1 = int(round(len(h) * 100.0 / fs))
        k2 = int(round(len(h) * 10000.0 / fs))
        B, A = invfreqz(w[k1 : k2 + 1], h[k1 : k2 + 1])

        R = sqrt(A[2])
        theta = math.acos(A[1] / (-2 * R))
        aa[i] = fs * theta
        frn = theta / (2 * pi)
        Q[i] = (pi * frn) / (1 - R)

        print "Pole frequency = %f Hz" % (aa[i] / (2 * pi))
        # print 'Q = %f' % Q[i]
        # print "R =", R
        # print "frn =", frn, "theta =", theta

        A = array((1, -2 * R * cos(theta), R * R))
        w1, h1 = freqz(Bconst, A, worN=15000)
        maxh[i] = max(abs(h))
        # print "gain =", gain[i]
    return aa, Q, maxh

def plot_filter(h, title, freq, gain, show=True):
    if h.ndim == 2:  # second-order sections
        sos = h
        n = mne.filter.estimate_ringing_samples(sos)
        h = np.zeros(n)
        h[0] = 1
        h = signal.sosfilt(sos, h)
        H = np.ones(512, np.complex128)
        for section in sos:
            f, this_H = signal.freqz(section[:3], section[3:])
            H *= this_H
    else:
        f, H = signal.freqz(h)
    fig, axs = plt.subplots(2)
    t = np.arange(len(h)) / sfreq
    axs[0].plot(t, h, color=blue)
    axs[0].set(xlim=t[[0, -1]], xlabel='Time (sec)',
               ylabel='Amplitude h(n)', title=title)
    box_off(axs[0])
    f *= sfreq / (2 * np.pi)
    axs[1].semilogx(f, 10 * np.log10((H * H.conj()).real), color=blue,
                    linewidth=2, zorder=4)
    plot_ideal(freq, gain, axs[1])
    mne.viz.tight_layout()
    if show:
        plt.show()

 def compute_frequencies(self, N=None):
     if hasattr(self, 'sample_rate'):
         try:
             self.W, self.H = signal.freqz(self.B, self.A)
         except:
             self.W, self.H = signal.freqz(self.B)
     else:
         self.W, self.H = signal.freqs(self.B, self.A, N)

    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 doplot2(B, A, B2, A2, freq = (315.0/88.0) * 8.0):
	w, h = sps.freqz(B, A)
	w2, h2 = sps.freqz(B2, A2)

#	h.real /= C
#	h2.real /= C2

	begin = 0
	end = len(w)
#	end = int(len(w) * (12 / freq))

#	chop = len(w) / 20
	chop = 0
	w = w[begin:end]
	w2 = w2[begin:end]
	h = h[begin:end]
	h2 = h2[begin:end]

	v = np.empty(len(w))
	
#	print len(w)

	hm = np.absolute(h)
	hm2 = np.absolute(h2)

	v0 = hm[0] / hm2[0]
	for i in range(0, len(w)):
#		print i, freq / 2 * (w[i] / pi), hm[i], hm2[i], hm[i] / hm2[i], (hm[i] / hm2[i]) / v0
		v[i] = (hm[i] / hm2[i]) / v0

	fig = plt.figure()
	plt.title('Digital filter frequency response')

	ax1 = fig.add_subplot(111)

	v  = 20 * np.log10(v )

#	plt.plot(w * (freq/pi) / 2.0, v)
#	plt.show()
#	exit()

	plt.plot(w * (freq/pi) / 2.0, 20 * np.log10(abs(h)), 'r')
	plt.plot(w * (freq/pi) / 2.0, 20 * np.log10(abs(h2)), 'b')
	plt.ylabel('Amplitude [dB]', color='b')
	plt.xlabel('Frequency [rad/sample]')
	
	ax2 = ax1.twinx()
	angles = np.unwrap(np.angle(h))
	angles2 = np.unwrap(np.angle(h2))
	plt.plot(w * (freq/pi) / 2.0, angles, 'g')
	plt.plot(w * (freq/pi) / 2.0, angles2, 'y')
	plt.ylabel('Angle (radians)', color='g')

	plt.grid()
	plt.axis('tight')
	plt.show()

def lab5_ex4():

    # define coefficient arrays
    b = firwin(20, 0.3)  # given 10th degree polynomial
    a = array([1.0])  # no feedback

    # define quantized versions of 10th order array at Q3, Q6, Q8, Q16 (second array will always be 1.0)
    q3 = fix(b * 2 ** 3) * 2 ** (
        -3
    )  # (convert decimal into integer at N precision, then convert back into float decimal at same precision)
    q6 = fix(b * 2 ** 6) * 2 ** (-6)
    q8 = fix(b * 2 ** 8) * 2 ** (-8)
    q16 = fix(b * 2 ** 16) * 2 ** (-16)

    # define and plot full precision and quantized frequency responses on linear and dB scales
    w, h = freqz(b, a)
    subplot(211)
    plot(w / pi, abs(h), "k-")
    subplot(212)
    plot(w / pi, 20 * log10(abs(h)), "k-")

    w, h = freqz(q3, a)
    subplot(211)
    plot(w / pi, abs(h), "b-")
    subplot(212)
    plot(w / pi, 20 * log10(abs(h)), "b-")

    w, h = freqz(q6, a)
    subplot(211)
    plot(w / pi, abs(h), "g-")
    subplot(212)
    plot(w / pi, 20 * log10(abs(h)), "g-")

    w, h = freqz(q8, a)
    subplot(211)
    plot(w / pi, abs(h), "c-")
    subplot(212)
    plot(w / pi, 20 * log10(abs(h)), "c-")

    w, h = freqz(q16, a)
    subplot(211)
    plot(w / pi, abs(h), "r-")
    legend(("full precision", "Q3", "Q6", "Q8", "Q16"))
    title("freq responses of 10th order filter by quantization and scale")
    ylabel("linear scale")
    subplot(212)
    plot(w / pi, 20 * log10(abs(h)), "r-")
    legend(("full precision", "Q3", "Q6", "Q8", "Q16"), loc=3)
    ylabel("dB scale")
    xlabel("normalized frequency")
    show()

    print (
        "\nFilter's responses are quite different depending on the number of bits used for quantization. When using 3 bits especially, the filter's response is markedly different from it's response to the full precision signal. We can clearly see this effect on the dB scale, where it looks as if the number of bits available for quantization corresponds to the number of zeros/poles in the transfer function."
    )

def butterfilt(data, sr, passband=(1.0, 50.0), stopband=(59., 61.),
               order=(5, 5), plotfreqz=False, plotlim=100):
    """Wrapper for Butterworth filter of sample*channel*trial EEG data

    Inputs:
    data - sample*channel*trial EEG data
    sr - sample rate (samps/sec)

    Optional Inputs:
    passband - low and high cutoffs for Butterworth passband
    stopband - low and high cuttoffs for Butterworth stopband
             Note that this should be used as a Notch filter
             For countries w/ 60 Hz power lines (Americas etc.) use (59., 61.)
             For countries w/ 50 Hz power lines (everywhere else) use (49., 51.)
    order - Order of both Butterworth filters: (passband, stopband)
    plotfreqz - Flag for plotting frequency responses of both filters
    plotlim - Upper limit of frequencies to plot

    Outputs:
    filtdata - Filtered sample*channel*trial EEG data
    """

    nyquist = .5 * float(sr)
    b, a = butter(order[0], [float(passband[0]) / nyquist,
                             float(passband[1]) / nyquist], btype='bandpass')
    filtdata = filtfilt(b, a, data, axis=0)

    if plotfreqz:
        w, h = freqz(b, a)
        plt.plot((nyquist / np.pi) * w, abs(h))
        plt.setp(plt.gca(),XLim=[0,plotlim],YLim=[0,1.1])
        plt.plot([0, nyquist], [np.sqrt(0.5), np.sqrt(0.5)], '--')
        plt.title(
            'Butterworth Passband Frequency Response, Order = %d' % order[1])

    if not not stopband:
        B, A = butter(order[1], [float(stopband[0]) / nyquist,
                      float(stopband[1]) / nyquist],
                      btype='bandstop')
        filtdata = filtfilt(B, A, data, axis=0)
        if plotfreqz:
            W, H = freqz(B, A)
            plt.figure()
            plt.plot((nyquist / np.pi) * W, abs(H))
            plt.setp(plt.gca(),XLim=[0,plotlim],YLim=[0,1.1])
            plt.plot([0, nyquist], [np.sqrt(0.5), np.sqrt(0.5)], '--')
            plt.title(
                'Butterworth Stopband Frequency Response, Order = %d'
                % order[1])

    return filtdata

def output_freqz_libcrybaby():
    from pluginloader import Plugin
    filt = Plugin("../build/default/src/plugins/libcrybaby.so")
    filt['crybaby2.refvolt'] = 0.1
    para = 'crybaby2.hotpotz'
    name = filt.get_var_attr(para)[0]
    fs = 44100
    filt.init(fs)
    impulse = zeros(10000,dtype=float32)
    impulse[0] = 1
    rg = log10(linspace(*np.power(10,filt.get_range(para)), num=20))
    if False:
        aa, Q, maxh, p1, Bconst, freq_const = estimate(rg, filt, para, impulse, fs)
        pickle.dump((aa,Q,maxh,p1,Bconst,freq_const), file("p.out","w"))
    else:
        aa, Q, maxh, p1, Bconst, freq_const = pickle.load(file("p.out"))
    z1, z2A, z3A, gcA = estimate_SR(filt, para, freq_const, impulse)
    off, a1A, qA = fit_param(rg, aa, Q, fs)
    gain, g_off, gA = make_gain(rg, off, a1A, qA, maxh, p1, Bconst, fs)
    #show_param(rg, off, g_off, aa, a1A, Q, qA, gain, gA)
    rg = log10(linspace(*np.power(10,filt.get_range(para)), num=10))
    filt.init(fs)
    for i, freq in enumerate(rg):
        filt[para] = freq
        filt.compute(zeros(10000,dtype=float32))
        w, h = freqz(filt.compute(impulse), worN=15000)

        q  = polyval(qA, freq)
        a1 = (off - 1 / polyval(a1A, freq)) / fs
        g = g_off - 1 / polyval(gA, freq)
        gc = polyval(gcA, fs)

        r = 1 - a1/(2*q)
        A = array((1, -2*r*cos(a1), r*r))
        A = poly(list(roots(A))+[p1])
        B = poly([polyval(z2A,fs),polyval(z3A,fs),z1])
        w1, h1 = freqz(B*g*gc, A, worN=15000)

        semilogx((w*fs/(2*pi))[1:],  20*log10(abs(h[1:])), "b")
        semilogx((w1*fs/(2*pi))[1:], 20*log10(abs(h1[1:])), "r")

    print "theta2pi = (%g - 1000 / (%s)) / SR;" % (off, string_polyval(a1A*1000, "wah"))
    print "Q = %s;" % string_polyval(qA, "wah")
    print "g = %g - 1 / (%s);" % (g_off, string_polyval(gA, "wah"))
    print "gc = %s;" % string_polyval(gcA, "SR")
    print "p1 = exp(-1000/(%g*SR));" % (-1000/(fs*log(p1)))
    print "z1 = %g;" % z1
    print "z2 = %s;" % string_polyval(z2A, "SR")
    print "z3 = %s;" % string_polyval(z3A, "SR")

    show()

    def fourier(self,List=None,fs=1):
        if type(List) is list:
            indices=self.getInds(List=List)
            if 0 not in indices:
                List.append((0,0))
            domain=[int(self.getDom(List=List)[0]),int(self.getDom(List=List)[1])]

            for i in range(domain[0],domain[1]+1,1):
                if i not in  indices:
                    List.append((i,0))

            indices=self.getInds(List=List)
            values=self.getVals(List=List)
            f_indices,f_values=freqz(values,indices,worN=2000)
            f_indices=f_indices*(fs * 0.5 / np.pi)
            f_List=[]
            f_indices=list(f_indices)
            values_copy=f_values.copy()
            values_copy=list(values_copy)
            f_values=list(f_values)
            i=0
            for value in f_values:
                f_List.append((f_indices[i],float(np.real(f_values[i]))))
                i=i+1
            return f_List
        else:
            List=self.points
            indices=self.getInds(List=List)
            if 0 not in indices:
                List.append((0,0))
            domain=[int(self.getDom(List=List)[0]),int(self.getDom(List=List)[1])]
            for i in range(domain[0],domain[1]+1,1):
                if i not in  indices:
                    List.append((i,0))

            indices=self.getInds(List=List)
            values=self.getVals(List=List)
            f_indices,f_values=freqz(values,indices,worN=2000)
            f_indices=f_indices*(fs * 0.5 / np.pi)
            f_List=[]
            f_indices=list(f_indices)
            values_copy=f_values.copy()
            values_copy=list(values_copy)
            f_values=list(f_values)
            i=0
            for value in f_values:
                f_List.append((f_indices[i],float(np.real(f_values[i]))))
                i=i+1
            return f_List

    def plot_filter_characteristics(self):
        w, h = freqz(self.freq_filter.num, self.freq_filter.denom)
        plt.figure(1)
        plt.subplot(2,1,1)
        plt.hold(True)
        powa = plt.plot((self.filter_parameters.sample_rate*0.5/pi)*w, abs(h),'b-', label = 'Char. amplitudowa')
        plt.title('Charakterystyki filtru')
        plt.xlabel('Czestotliwosc [Hz]')
        plt.ylabel('Amplituda')


        plt.twinx(ax=None)
        angles = unwrap(angle(h))
        plt.znie = plot((self.filter_parameters.sample_rate*0.5/pi)*w,angles, 'g-', label = 'Char. fazowa')
        plt.ylabel('Faza')

        plt.grid()
        tekst = powa + znie
        wybierz = [l.get_label() for l in tekst]

        plt.legend(tekst, wybierz, loc='best')
    ########################################################################################################################
        plt.subplot(2,1,2)

        w2, gd = group_delay((num, denom))

        plt.plot((sample_rate*0.5/pi)*w2, gd)
        plt.grid()
        plt.xlabel('Czestotliwosc [Hz]')
        plt.ylabel('Opoznienie grupowe [probki]')
        plt.title('Opoznienie grupowe filtru')

        plt.show()

    def __init__(self, fir_len, fir_bits, tb_width):
        tb_ctr = 1/(2*4)
        pass_corner = tb_ctr - (tb_ctr*tb_width/2)
        stop_corner = tb_ctr + (tb_ctr*tb_width/2)
        fir_bands = [0, pass_corner, stop_corner, 0.5]
        b = signal.remez(fir_len, fir_bands, [1, 0])
        coeff_scl = 2**(fir_bits-1)
        self.fir_coeff = np.floor(b*coeff_scl + 0.5)
        
        # Dump Coefficients?
        if 1:
            write_meminit("fir4dec_coeff.v", self.fir_coeff)

        self.fir_coeff = self.fir_coeff/coeff_scl;
        
        # plot FIR response?
        if 1:
            W, H = signal.freqz(self.fir_coeff)
            plt.figure()
            plt.plot(W/(2*np.pi), 20*np.log10(np.abs(H)))
            plt.grid()
            plt.xlabel("Freq (normalized)")
            plt.ylabel("dB")
            plt.title("fir4dec response (close to continue sim)")
            plt.show()

 def frequency_response(self, n=10000):
     """
     Generate a filter frequency response from a set of filter taps.
     Returns plottable (x, y) with respect to an actual sampling rate
     """
     w, h = signal.freqz(self.b, self.a, worN=n)
     return (0.5 * self.fs * w / np.pi, np.abs(h))

def main():
    edges = [30, 60, 120, 240]
    corners = zip(edges[:-1], edges[1:])
    centres = [(a + b) / 2 for a, b in corners]

#c = [get_linkwitz_riley_coeffs(1, b, a, edges[-1] * 2) for b, a in corners]
    sr = 2000
    c = [get_peak_coeffs(-24, i, sr, 1) for i in centres]
    c.append([[1, 0, 0], [1, 0, 0]])

    bm = [BiquadMemory(0, 0) for _ in c]
    bc = [BiquadCoefficients(b0, b1, b2, a1, a2)
          for [b0, b1, b2], [a0, a1, a2] in c]

    c.append(series_coeffs(c))

#    c.append(impedance_filter(c[-1]))

    wh = [signal.freqz(b, a) for b, a in c]

    plt.subplot(111)
    plt.title("Frequency response - reflection filter")
    for w, h in wh:
        plt.semilogx(w, 20 * np.log10(np.abs(h)))
    plt.ylabel('Amplitude Response (dB)')
    plt.xlabel('Frequency (rad/sample)')
    plt.grid()

    plt.show()

def filter_VE(data):
    f=np.linspace(0,1,50000)
    a = np.ones(50000)
    
    #Switching power supply frequency (~290 KHz)
    a[2850:2950] =0
    
    #1.49 MHz
    a[14850:14950]=0
    
    #AM 1.37 MHz
    a[13650:13750] = 0
    
    #80 m Ham band (3.97 MHz)
    a[39650:39750] = 0
    
    #80 m Ham band (4 MHz)
    a[39950:40050]= 0
    
    a[-1]=0
    
    b=sigs.firwin2(3000,f,a)
    [h, fpoints] = sigs.freqz(b, 1, 50000,10E6)
    
    #Run the FIR filter
    vec = sigs.lfilter(b,1,data)
    return vec

    def test_firls(self):
        N = 11  # number of taps in the filter
        a = 0.1  # width of the transition band

        # design a halfband symmetric low-pass filter
        h = firls(11, [0, a, 0.5-a, 0.5], [1, 1, 0, 0], fs=1.0)

        # make sure the filter has correct # of taps
        assert_equal(len(h), N)

        # make sure it is symmetric
        midx = (N-1) // 2
        assert_array_almost_equal(h[:midx], h[:-midx-1:-1])

        # make sure the center tap is 0.5
        assert_almost_equal(h[midx], 0.5)

        # For halfband symmetric, odd coefficients (except the center)
        # should be zero (really small)
        hodd = np.hstack((h[1:midx:2], h[-midx+1::2]))
        assert_array_almost_equal(hodd, 0)

        # now check the frequency response
        w, H = freqz(h, 1)
        f = w/2/np.pi
        Hmag = np.abs(H)

        # check that the pass band is close to unity
        idx = np.logical_and(f > 0, f < a)
        assert_array_almost_equal(Hmag[idx], 1, decimal=3)

        # check that the stop band is close to zero
        idx = np.logical_and(f > 0.5-a, f < 0.5)
        assert_array_almost_equal(Hmag[idx], 0, decimal=3)

    def test_hilbert(self):
        N = 11  # number of taps in the filter
        a = 0.1  # width of the transition band

        # design an unity gain hilbert bandpass filter from w to 0.5-w
        h = remez(11, [a, 0.5-a], [1], type='hilbert')

        # make sure the filter has correct # of taps
        assert_(len(h) == N, "Number of Taps")

        # make sure it is type III (anti-symmetric tap coefficients)
        assert_array_almost_equal(h[:(N-1)//2], -h[:-(N-1)//2-1:-1])

        # Since the requested response is symmetric, all even coeffcients
        # should be zero (or in this case really small)
        assert_((abs(h[1::2]) < 1e-15).all(), "Even Coefficients Equal Zero")

        # now check the frequency response
        w, H = freqz(h, 1)
        f = w/2/np.pi
        Hmag = abs(H)

        # should have a zero at 0 and pi (in this case close to zero)
        assert_((Hmag[[0, -1]] < 0.02).all(), "Zero at zero and pi")

        # check that the pass band is close to unity
        idx = np.logical_and(f > a, f < 0.5-a)
        assert_((abs(Hmag[idx] - 1) < 0.015).all(), "Pass Band Close To Unity")