def mgr_filter(mgr, wp, ws, gpass, gstop, analog=0, ftype='ellip', output='ba', unit='hz', use_filtfilt=False, meancorr=1.0):
    if unit == 'radians':
        b,a = signal.iirdesign(wp, ws, gpass, gstop, analog, ftype, output)
    elif unit == 'hz':
        nyquist = float(mgr.get_param('sampling_frequency'))/2.0
        try:
            wp = wp/nyquist
            ws = ws/nyquist
        except TypeError:
            wp = [i/nyquist for i in wp]
            ws = [i/nyquist for i in ws]

        b,a = signal.iirdesign(wp, ws, gpass, gstop, analog, ftype, output)
    if use_filtfilt:    
        from scipy.signal import filtfilt
        #samples_source = read_data_source.MemoryDataSource(mgr.get_samples(), False)
        for i in range(int(mgr.get_param('number_of_channels'))):
            print("FILT FILT CHANNEL "+str(i))
            mgr.get_samples()[i,:] = signal.filtfilt(b, a, mgr.get_samples()[i]-np.mean(mgr.get_samples()[i])*meancorr)
        samples_source = read_data_source.MemoryDataSource(mgr.get_samples(), False)
    else:
        print("FILTER CHANNELs")
        filtered = signal.lfilter(b, a, mgr.get_samples())
        print("FILTER CHANNELs finished")
        samples_source = read_data_source.MemoryDataSource(filtered, True)

    info_source = copy.deepcopy(mgr.info_source)
    tags_source = copy.deepcopy(mgr.tags_source)
    new_mgr = read_manager.ReadManager(info_source, samples_source, tags_source)
    return new_mgr

    def __init__( self, band_start, band_stop ):
        nyquist_frequency = float(SAMPLES_PER_SECOND) / 2.0
        
        band_start /= nyquist_frequency
        band_stop /= nyquist_frequency

        assert( band_start >= 0 and band_start <= 1 )
        assert( band_stop >= 0 and band_stop <= 1 )
        assert( band_stop >= band_start )

        passband_edges = []
        stopband_edges = []

        if band_start >= 0.05: # if not, make LPF only
            passband_edges.append( band_start * 1.025 )
            stopband_edges.append( band_start * 0.975 )

        if band_stop <= 0.95: # if not, make HPF only
            passband_edges.append( band_stop * 0.975 )
            stopband_edges.append( band_stop * 1.025 )

        (self.feedforward_taps,
         self.feedback_taps) = iirdesign( passband_edges,
                                          stopband_edges,
                                          0.1,               # max attenuation (dB) in passband
                                          30 )               # min attenuation (dB) in stopband

        self.filter_state = lfiltic( self.feedforward_taps, self.feedback_taps, [] )

def iir_bandstops(fstops, fs, order=4):
    """ellip notch filter
    fstops is a list of entries of the form [frequency (Hz), df, df2]                           
    where df is the pass width and df2 is the stop width (narrower                              
    than the pass width). Use caution if passing more than one freq at a time,                  
    because the filter response might behave in ways you don't expect.
    """
    nyq = 0.5 * fs

    # Zeros zd, poles pd, and gain kd for the digital filter
    zd = np.array([])
    pd = np.array([])
    kd = 1

    # Notches
    for fstopData in fstops:
        fstop = fstopData[0]
        df = fstopData[1]
        df2 = fstopData[2]
        low = (fstop - df) / nyq
        high = (fstop + df) / nyq
        low2 = (fstop - df2) / nyq
        high2 = (fstop + df2) / nyq
        z, p, k = iirdesign([low,high], [low2,high2], gpass=1, gstop=6,
                            ftype='ellip', output='zpk')
        zd = np.append(zd,z)
        pd = np.append(pd,p)

    # Set gain to one at 100 Hz...better not notch there                                        
    bPrelim,aPrelim = zpk2tf(zd, pd, 1)
    outFreq, outg0 = freqz(bPrelim, aPrelim, 100/nyq)

    # Return the numerator and denominator of the digital filter                                
    b,a = zpk2tf(zd,pd,k)
    return b, a

def iir_bandstops(fstops, fs, order=4):
    import numpy as np
    from scipy.signal import iirdesign, zpk2tf, freqz
    nyq = 0.5 * fs
    # Zeros zd, poles pd, and gain kd for the digital filter
    zd = np.array([])
    pd = np.array([])
    kd = 1

    # Notches
    for fstopData in fstops:
        fstop = fstopData[0]
        df = fstopData[1]
        df2 = fstopData[2]
        low = (fstop - df) / nyq
        high = (fstop + df) / nyq
        low2 = (fstop - df2) / nyq
        high2 = (fstop + df2) / nyq
        z, p, k = iirdesign([low,high], [low2,high2], gpass=1, gstop=6,
                            ftype='ellip', output='zpk')
        zd = np.append(zd,z)
        pd = np.append(pd,p)

    # Set gain to one at 100 Hz...better not notch there
    bPrelim,aPrelim = zpk2tf(zd, pd, 1)
    outFreq, outg0 = freqz(bPrelim, aPrelim, 100/nyq)

    # Return the numerator and denominator of the digital filter
    b,a = zpk2tf(zd,pd,k)
    return b, a

    def iir(self):
        """
        Filter the time-series using an IIR filter. Filtering is done back and
        forth (using scipy.signal.filtfilt) to achieve zero phase delay

        """

        #Passband and stop-band are expressed as fraction of the Nyquist
        #frequency:
        if self.ub is not None:
            ub_frac = self.ub / (self.sampling_rate / 2.)
        else:
            ub_frac = 1.0

        lb_frac = self.lb / (self.sampling_rate / 2.)

        wp = [lb_frac, ub_frac]

        #Make sure to not exceed the interval 0-1:
        ws = [np.max([lb_frac - 0.1, 0]),
              np.min([ub_frac + 0.1, 1.0])]

        b, a = signal.iirdesign(wp, ws, self._gpass, self._gstop,
                                ftype=self._ftype)

        return self.filtfilt(b, a)

  def update_design(self):
    ax = self.ax
    ax.cla()
    ax2 = self.ax2
    ax2.cla()

    wp = self.wp
    ws = self.ws
    gpass = self.gpass
    gstop = self.gstop
    b, a = ss.iirdesign(wp, ws, gpass, gstop, ftype=self.ftype, output='ba')
    self.a = a
    self.b = b
    #b = [1,2]; a = [1,2]
    #Print this on command line so we can use it in our programs
    print 'b = ', pylab.array_repr(b)
    print 'a = ', pylab.array_repr(a)

    my_w = pylab.logspace(pylab.log10(.1*self.ws[0]), 0.0, num=512)
    #import pdb;pdb.set_trace()
    w, h = freqz(b, a, worN=my_w*pylab.pi)
    gp = 10**(-gpass/20.)#Go from db to regular
    gs = 10**(-gstop/20.)
    self.design_line, = ax.plot([.1*self.ws[0], self.ws[0], wp[0], wp[1], ws[1], 1.0], [gs, gs, gp, gp, gs, gs], 'ko:', lw=2, picker=5)
    ax.semilogx(w/pylab.pi, pylab.absolute(h),lw=2)
    ax.text(.5,1.0, '{:d}/{:d}'.format(len(b), len(a)))
    pylab.setp(ax, 'xlim', [.1*self.ws[0], 1.2], 'ylim', [-.1, max(1.1,1.1*pylab.absolute(h).max())], 'xticklabels', [])

    ax2.semilogx(w/pylab.pi, pylab.unwrap(pylab.angle(h)),lw=2)
    pylab.setp(ax2, 'xlim', [.1*self.ws[0], 1.2])
    ax2.set_xlabel('Normalized frequency')

    pylab.draw()

    def __init__(self, wp, ws, gpass, gstop):

        b,a = iirdesign(wp, ws, gpass, gstop)
        print b,a
        input_weights = b[::-1]
        feedback_weights = a[::-1]
        LinearFilter.__init__(self, input_weights=input_weights, feedback_weights=feedback_weights)

    def test_filter(self):
        if not self.TEST_FILTER:
            return


        def getSin(hz, sampling, data):
            return sin(hz*2*pi*data/sampling)
        mgr = get_fake_manager(3, 1000)
        x=r_[0:1000]
        y1 = getSin(2, 128.0, x)
        y2 = getSin(15, 128.0, x)
        y3 = getSin(30, 128.0, x)
        mgr.get_samples()[0,:] = y1
        mgr.get_samples()[1,:] = y2
        mgr.get_samples()[2,:] = y3

        b,a = signal.iirdesign([9/128.0, 16/128.0], [4/128.0, 25/128.0],0.2, 70.0, ftype='cheby2')
        y = signal.lfilter(b, a, mgr.get_samples())


        e = ch.Filter([9/128.0, 16/128.0], [4/128.0, 25/128.0], 0.2, 70.0, ftype='cheby2')
        new_mgr = e.process([mgr])[0]
        ch.Plot('CH1').process([new_mgr])

        self.assertAlmostEqual(y[2,5], new_mgr.get_samples()[2,5])
        self.assertAlmostEqual(y[0,100], new_mgr.get_samples()[0,100])
        self.assertNotAlmostEqual(mgr.get_samples()[0,100], new_mgr.get_samples()[0,100])
        LOGGER.info("FILTER tested!")

	def set_params(self, **kwargs):
		f = kwargs['cut_off_freq']
		normalized_pass = f/(RATE*.5)
		normalized_stop = (f+.3*f)/(RATE*.5)
		(self.b, self.a) = signal.iirdesign(normalized_pass, normalized_stop, 
											LowpassFilter.pass_band_loss, 
											LowpassFilter.stop_band_loss)

def lab6_ex4():
    # set parameters of system    
    fs = 2000 #sampling frequency (Hz)
    fn = fs/2.0 #nyqvist frequency
    wp = array([200/fn,400/fn]) #band-pass in normalized corner frequencies
    ws = array([1/fn,100/fn,500/fn,5000/fn]) #band-stop in normalized corner frequencies
    gstop = 40 #minimum attenuation in stopband
    gpass = 1 #maximum loss in passband
    
    # create remez filter, frequency response, and plot
    b,a = iirdesign(wp,ws,gpass,gstop,analog=0)
    w,h = freqz(b,a)
    plot(w/pi, 20*log10(abs(h)),'b-')
    title('IIR of given frequency response')
    xlabel('normalized frequency (1 = fn)')
    ylabel('magnitude (dB scale)')
    grid()
    show()

    print("\nBoth designs achieve the desired response at the desired attenuation, but the IIR design's response attenuates and gains more quickly than the 46th-order remez FIR design.")

    z,p,k = tf2zpk(b,a)
    zplane(z,p)
    title('zplane of IIR of given frequency response')
    show()

    print('\nZ-plane analysis shows that the IIR design uses an 8th degree polynomial, while the remez FIR uses a 46th order polynomial (and thus far more components) to achieve the same desired response. This is because the IIR filter is able to use both zeros and poles (i.e. feedback) to achieve the desired frequency response at the desired attenuation and loss, whereas the remez FIR must use only zeros (no feedback).')

 def _design_filter(self, FS):
     if not FS in self._coefs_cache:
         wp, ws = self.fp*2/FS, self.fs*2/FS
         b,a = signal.iirdesign(wp=wp, ws=ws, gstop=self.gstop, gpass=self.gpass, ftype=self.ftype)
         self._coefs_cache[FS]=(b,a)
     else:
         b,a = self._coefs_cache[FS]
     return b, a  

 def getFilter(self):
     if self.__filter__ == None:
         passband = [calcFreq(self.A, -10)/sampleRate,
                     calcFreq(self.A, 1)/sampleRate]
         stopband = [calcFreq(self.A, -11)/sampleRate,
                     calcFreq(self.A, 2)/sampleRate]
         self.__filter__ = iirdesign(passband, stopband, 5, 10, ftype="cheby1")
     return self.__filter__

 def designIIR(self):
     """Design an iir filter with a passband from .3 to .35 using iirdesign from scipy
     """
     #note - the default filter returned with the default ftype for these
     #parameters was UNSTABLE leading to random unit test failures with nans 
     #cheby2 returns a stable filter
     b, a = iirdesign([.3, .35], [.1, .5],.1,20, ftype='cheby2')
     return b.tolist(), a.tolist()

def butterfilt(finname, foutname, fmt, fs, fl=5.0, fh=100.0, gpass=1.0, gstop=30.0, ftype='butter', buffer_len=100000, overlap_len=100, max_len=-1):
  """Given sampling frequency, low and high pass frequencies design a butterworth filter and filter our data with it."""
  fso2 = fs/2.0
  wp = [fl/fso2, fh/fso2]
  ws = [0.8*fl/fso2,1.4*fh/fso2]
  import pdb; pdb.set_trace()
  b, a = iirdesign(wp, ws, gpass=gpass, gstop=gstop, ftype=ftype, output='ba')
  y = filtfiltlong(finname, foutname, fmt, b, a, buffer_len, overlap_len, max_len)
  return y, b, a

	def get_effected_signal(self, sig):
		if self.it >= 1000:
			self.it = 0
		freq = self.f + self.sin_wave[self.it]
		self.it += 1
		normalized_pass = freq/(RATE*.5)
		normalized_stop = (freq+.3*freq)/(RATE*.5)
		(a, b) = signal.iirdesign(normalized_pass, normalized_stop, 1, 30)
		out = signal.lfilter(b, a, sig)
		return out / np.max(out)

def continuous_bandpass_filter(left, right, freq_low, freq_high, framerate):
	pass_region = [0.075, 0.2]
	stop_region = [0.05, 0.8]
	pass_max_damp = 1
	stop_min_damp = 80
	b, a = signal.iirdesign(pass_region, stop_region, pass_max_damp, stop_min_damp)
	w, h = signal.freqz(b, a)
	new_left = signal.lfilter(b, a, left)
	new_right = signal.lfilter(b, a, right)
	return new_left, new_right

def butterworth_bandpass_filter(left_channel, right_channel, pass_region, stop_region, framerate):
	pass_region = 2*pass_region/framerate
	stop_region = 2*stop_region/framerate
	#print(pass_region, stop_region)
	pass_max_damp = 2
	stop_min_damp = 20
	b, a = signal.iirdesign(pass_region, stop_region, pass_max_damp, stop_min_damp)
	new_left_channel = signal.lfilter(b, a, left_channel)
	new_right_channel = signal.lfilter(b, a, right_channel)
	return new_left_channel, new_right_channel

	def getFramePitchAdvanced(self):
		if len(self.speechSegment)==0 or self.isBlank==True:
			return
		
		self.p = []
		b, a = signal.iirdesign([60.0*2/self.sampleRate,950.0*2/self.sampleRate],[50.0*2/self.sampleRate,1000.0*2/self.sampleRate],2,40)
		for frame in self.frame:
			filt = signal.lfilter(b,a,frame)
			minAMDF = utils.AMDF(filt)
			pitch = self.sampleRate/minAMDF
			self.p.append(pitch)

def naive_dsp(left, right, freq_low, freq_high, framerate):
	low_pass = 0.05
	low_stop = 0.1
	high_pass = low_stop
	high_stop = low_pass
	pass_max_damp = 2
	stop_min_damp = 40

	b1, a1 = signal.iirdesign(low_pass, low_stop, pass_max_damp, stop_min_damp)
	b2, a2 = signal.iirdesign(high_pass, high_stop, pass_max_damp, stop_min_damp)

	left = signal.lfilter(b1, a1, left)
	right = signal.lfilter(b1, a1, right)

	new_left = left - right
	new_right = left - right

	left = signal.lfilter(b2, a2, new_left)
	right = signal.lfilter(b2, a2, new_right)

	return left, right

def makeFIRCoefs(samplingFreq, IFFreq, bandwidth, numPoints, numTaps):
    IFFreq = IFFreq * 1e06
    (b, a) = signal.iirdesign(IFFreq / 2, IFFreq, 3, 30, ftype="butter")
    timeStep = 1.0 / samplingFreq
    timePts = np.arange(0, numPoints * timeStep, timeStep)
    stepImpulse = [0 if i < len(timePts) / 2 else 1 for i in range(len(timePts))]
    unitImpulse = [1 if i == 0 else 0 for i in range(len(timePts))]
    nIFFreq = IFFreq / (samplingFreq / 2)
    nbandwidth = bandwidth / (samplingFreq / 2)
    print "Cuttoff frequency: ", nIFFreq, "Hz\n"
    print "building FIR coefficients with {} taps...\n".format(numTaps)
    firCoefs = signal.firwin(
        float(numTaps),
        0.09 * samplingFreq,
        width=None,
        window="hamming",
        pass_zero=True,
        scale=True,
        nyq=samplingFreq / 2,
    )

    # step impulse for polyphase filter input test:
    with open("/home/nick/firFilter/stepFunction" + str(numPoints) + ".txt", "w") as stepfid:
        stepfid.write("signal myStepImpulse : vector_array :=(")
        stepfid.write(",".join(['x"{0:04X}"'.format(x) for x in stepImpulse]) + ");")
    with open("/home/nick/firFilter/unitImpulse" + str(numPoints) + ".txt", "w") as impfid:
        impfid.write("signal myUnitImpulse : vector_array :=(")
        impfid.write(",".join(['x"{0:04X}"'.format(x) for x in unitImpulse]) + ");")

        # for i in range(len(firCoefs)):
        # 	print 'Reload index {} = {number:.{width}f}'.format(i,number=firCoefs[i],width=32)

        # two's complement filter coefficients for VHDL:
    with open("/home/nick/firFilter/firCoefs" + str(int(numTaps)) + "taps.coe", "w") as fid:
        fid.write("radix=10;\ncoefdata=")
        fid.write(",".join(["{number:.{width}f}".format(number=f, width=32) for f in firCoefs]) + ";")
    with open("/home/nick/firFilter/firCoefs" + str(IFFreq / 1.0e06) + "MHzVHDL.txt", "w") as fid2:
        fid2.write(
            "type vector_arrayTaps is array ("
            + str(numTaps - 1)
            + " downto 0) of std_logic_vector(15 downto 0);\nsignal myReload : vector_arrayTaps:=("
        )
        fid2.write(
            ",".join(['x"{0:04X}"'.format((~int(np.abs(x * pow(2, 16 - 1))) & 0xFFF) + 1) for x in firCoefs]) + ");"
        )

        # IIR filter coefficients written to VHDL:
    with open("/home/nick/firFilter/iirCoefs" + str(int(IFFreq / 1.0e06)) + "MHzVHDL.txt", "w") as iirfid:
        iirfid.write("signal loadIIRCoefsA : vector_array :=(")
        iirfid.write(",".join(['x"{0:04X}"'.format((~int(np.abs(x * pow(2, 64 - 1))) & 0xFFFF) + 1) for x in a]) + ");")
        iirfid.write("\nsignal loadIIRCoefsB : vector_array :=(")
        iirfid.write(",".join(['x"{0:04X}"'.format((~int(np.abs(x * pow(2, 64 - 1))) & 0xFFFF) + 1) for x in b]) + ");")