def makeFiltered(self):
     filtered = np.zeros(len(self.data))
     b1, a1 = butter(1, 0.0003, 'lowpass')
     filtered = filtfilt(b1, a1, self.data)
     filtered = self.data - filtered
     b2, a2 = butter(1, 0.025, 'lowpass')
     self.filt_data = filtfilt(b2, a2, filtered)

def filter_feepmt(feep):
    """
    input: an instance of class FeePmt
    output: buttersworth parameters of the equivalent FEE filter
    """

    # high pass butterswoth filter ~1/RC
    b1, a1 = signal.butter(1, feep.freq_LHPFd, 'high', analog=False)
    b2, a2 = signal.butter(1, feep.freq_LHPFd, 'low', analog=False)

    b0 = b2*(feep.Zin/(1+feep.C1/feep.C2))+b1*(feep.Zin*feep.R1/(feep.Zin+feep.R1))
    a0 = a1

    # LPF order 1
    b1l, a1l = signal.butter(1, feep.freq_LPF1d, 'low', analog=False)
    # LPF order 4
    b2l, a2l = signal.butter(4, feep.freq_LPF2d, 'low', analog=False)
    # convolve HPF, LPF1
    a_aux = np.convolve(a0, a1l, mode='full')
    b_aux = np.convolve(b0, b1l, mode='full')
    # convolve HPF+LPF1, LPF2
    a = np.convolve(a_aux, a2l, mode='full')
    b_aux2 = np.convolve(b_aux, b2l, mode='full')
    b = feep.A2*b_aux2

    return b, a

	def __init__(self,window_size, low, high):
		fs_Hz = 250;
		fn = fs_Hz/2
		self.filtered_data = np.array((window_size,1))
		

		#######################################
		# Filter Creation
		# -------------------------------------
		#
		# Create a filter using the scipy module,
		# based on specifications suggested by
		# Pan-Tompkins (bandpass from 5-15Hz)
		#
		#
		# 1) Establish constants:
		#		a) filter_order = 2
		#		b) high pass cutoff = 15Hz
		#		c) low pass cutoff = 5Hz
		# 2) Calculate the coefficients, store in variables

		filter_order = 2
		f_high = high
		f_low = low
		self.high_pass_coefficients = signal.butter(filter_order,f_low/fn, 'high')
		self.low_pass_coefficients = signal.butter(filter_order,f_high/fn, 'low')

    def __init__(self,source, nchannels, order, fc, btype='low'):
        Wn = np.asarray(np.atleast_1d(fc)).copy() #Scalar inputs are converted to 1-dimensional arrays
        self.samplerate = source.samplerate
        Wn= Wn/float(self.samplerate)*2    # wn=1 corresponding to half the sample rate

        if btype=='low' or btype=='high':
            self.filt_b=np.zeros((nchannels,order+1))
            self.filt_a=np.zeros((nchannels,order+1))
            if len(Wn)==1:     #if there is only one Wn value for all channel just repeat it
                self.filt_b, self.filt_a = signal.butter(order, Wn, btype=btype)
                self.filt_b=np.kron(np.ones((nchannels,1)),self.filt_b)
                self.filt_a=np.kron(np.ones((nchannels,1)),self.filt_a)
            else:               #else make nchannels different filters
                for i in range((nchannels)):
                    self.filt_b[i,:], self.filt_a[i,:] = signal.butter(order, Wn[i], btype=btype)
        else:
            self.filt_b=np.zeros((nchannels,2*order+1))
            self.filt_a=np.zeros((nchannels,2*order+1))
            if Wn.ndim==1:     #if there is only one Wn pair of values for all channel just repeat it
                self.filt_b, self.filt_a = signal.butter(order, Wn, btype=btype)
                self.filt_b=np.kron(np.ones((nchannels,1)),self.filt_b)
                self.filt_a=np.kron(np.ones((nchannels,1)),self.filt_a)
            else:   
                for i in range((nchannels)):
                    self.filt_b[i,:], self.filt_a[i,:] = signal.butter(order, Wn[:,i], btype=btype)   
                    
                    
        self.filt_a=self.filt_a.reshape(self.filt_a.shape[0],self.filt_a.shape[1],1)
        self.filt_b=self.filt_b.reshape(self.filt_b.shape[0],self.filt_b.shape[1],1)  
        self.nchannels = nchannels    
        LinearFilterbank.__init__(self,source, self.filt_b, self.filt_a) 

    def __init__(self, source, nchannels, order, fc, btype="low"):
        Wn = fc.copy()
        Wn = atleast_1d(Wn)  # Scalar inputs are converted to 1-dimensional arrays
        self.samplerate = source.samplerate
        try:
            Wn = Wn / self.samplerate * 2 + 0.0  # wn=1 corresponding to half the sample rate
        except DimensionMismatchError:
            raise DimensionMismatchError("Wn must be in Hz")

        if btype == "low" or btype == "high":
            self.filt_b = zeros((nchannels, order + 1))
            self.filt_a = zeros((nchannels, order + 1))
            if len(Wn) == 1:  # if there is only one Wn value for all channel just repeat it
                self.filt_b, self.filt_a = signal.butter(order, Wn, btype=btype)
                self.filt_b = kron(ones((nchannels, 1)), self.filt_b)
                self.filt_a = kron(ones((nchannels, 1)), self.filt_a)
            else:  # else make nchannels different filters
                for i in xrange((nchannels)):
                    self.filt_b[i, :], self.filt_a[i, :] = signal.butter(order, Wn[i], btype=btype)
        else:
            self.filt_b = zeros((nchannels, 2 * order + 1))
            self.filt_a = zeros((nchannels, 2 * order + 1))
            if Wn.ndim == 1:  # if there is only one Wn pair of values for all channel just repeat it
                self.filt_b, self.filt_a = signal.butter(order, Wn, btype=btype)
                self.filt_b = kron(ones((nchannels, 1)), self.filt_b)
                self.filt_a = kron(ones((nchannels, 1)), self.filt_a)
            else:
                for i in xrange((nchannels)):
                    self.filt_b[i, :], self.filt_a[i, :] = signal.butter(order, Wn[:, i], btype=btype)

        self.filt_a = self.filt_a.reshape(self.filt_a.shape[0], self.filt_a.shape[1], 1)
        self.filt_b = self.filt_b.reshape(self.filt_b.shape[0], self.filt_b.shape[1], 1)
        self.nchannels = nchannels
        LinearFilterbank.__init__(self, source, self.filt_b, self.filt_a)

 def butter_high_low_pass(lowcut, highcut, sampling_rate, order=5):
     nyq_freq = sampling_rate*0.5
     lower_bound = lowcut/nyq_freq
     higher_bound = highcut/nyq_freq
     b_high, a_high = butter(order, lower_bound, btype='high')
     b_low, a_low = butter(order, higher_bound, btype='low')
     return b_high, a_high, b_low, a_low

  def filter_data(self,eeg_data):
    #FILTER CONSTANTS
    fs = self.fs
    fn = self.fn
    filter_order = 2   #2nd order filter
    f_high = 50
    f_low = 5
    wn = [59,61]       #Nyquist filter window

    [b,a] = signal.butter(filter_order,f_high/fn, 'low')
    [b1,a1] = signal.butter(filter_order,f_low/fn, 'high')
    [bn,an] = signal.butter(4,[x/fn for x in wn], 'stop')

    filtered_eeg = []
    spectogram = []
    notched = []
    high_passed = []
    low_passed = []
    print(eeg_data)
    for i in range(len(eeg_data[0])):
      channel =  eeg_data[:,i]
      high_passed = signal.filtfilt(b1,a1,channel);        # high pass filter
      low_passed = signal.filtfilt(b,a,high_passed);                # low pass filter
      y = signal.filtfilt(bn,an,low_passed);        # notch filter
      filtered_eeg.append(y);
    self.filtered_eeg = filtered_eeg
    return filtered_eeg


# if __name__ == '__main__':
#   main()

    def filter(self):
        path = os.getcwd()+'/trialGraspEventDetection_dataFiles'
        self.Fgr = np.sum(self.values[:, 9:15], axis=1)  # SAI
        self.Fgl = np.sum(self.values[:, 0:7], axis=1)  # SAI

        # can use this to plot in matlab graspeventdetection_plot.m
        np.savetxt(path+'/SAI_Fgr.txt', self.Fgr)
        # can use this to plot in matlab
        np.savetxt(path+'/SAI_Fgl.txt', self.Fgl)

        # 0.55*pi rad/samples
        b1, a1 = signal.butter(1, 0.55, 'high', analog=False)
        self.f_acc_x = signal.lfilter(b1, a1, self.acc_x, axis=-1, zi=None)
        self.f_acc_y = signal.lfilter(b1, a1, self.acc_y, axis=-1, zi=None)
        self.f_acc_z = signal.lfilter(b1, a1, self.acc_z, axis=-1, zi=None)
        # self.f_eff = signal.lfilter(b1, a1, self.eff, axis=-1, zi=None)
        # type(eff)
        self.FAII = np.sqrt(np.square(self.f_acc_x) +
                            np.square(self.f_acc_y) +
                            np.square(self.f_acc_z))
        # can use this to plot in matlab
        np.savetxt(path+'/FAII.txt', self.FAII)

        # subtract base values from the values array
        self.values1 = self.values - self.values.min(axis=0)
        # pass the filter for each sensor
        self.fvalues1 = np.zeros(self.values1.shape)
        # 0.48*pi rad/samples
        b, a = signal.butter(1, 0.48, 'high', analog=False)
        for i in range(16):
            self.fvalues1[:, i] = signal.lfilter(b, a, self.values1[:, i],
                                                 axis=-1, zi=None)
        self.FAI = np.sum(self.fvalues1, axis=1)
        # can use this to plot in matlab
        np.savetxt(path+'/FAI.txt', self.FAI)

def butter_bandpass(fs, lowcut=None, highcut=None, order=10):
    """
    Parameters
    ----------
    fs : float
        Sample Rate
    lowcut : float
        Low pass frequency cutoff
    highcut : float
        High pass frequency cutoff
    order : int
        Butterworth filter order [Default = 10, i.e., 10th order Butterworth filter]

    Reference
    ---------
    http://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
    """
    from scipy.signal import butter
    nyq = 0.5 * fs
    if lowcut and highcut:
        high = highcut / nyq
        low = lowcut / nyq
        sos = butter(order, [low, high], analog=False, btype='band', output='sos')
    elif highcut:
        high = highcut / nyq
        sos = butter(order, high, btype='low', output='sos')
    elif lowcut:
        low = lowcut / nyq
        sos = butter(order, low, btype='high', output='sos')
    else:
        print('Error -- must supply lowcut, highcut or both')
        sos = None

    return sos

def temporally_filter(spatial, tfiltN, flash):
    print("Temporally Filtering")
    padlen=0
    if tfiltN[2]==1: #%if only high pass temporal filter
        [b,a]= butter(tfiltN[0],tfiltN[1],btype='highpass')
        padlen=3
    else: #%if band pass temporal filter
        [b,a]=butter(tfiltN[0],[tfiltN[1],tfiltN[2]], btype='bandpass')
        padlen=6
    temporal=np.array(spatial, dtype = np.float32)

    if len(flash) == 2:
        temporal[flash[0]:flash[1]+1,:,:]=np.mean(spatial[4:flash[0]+1, :, :], 0)
        bg=np.mean(temporal[:flash[0] + 1, :, :],0)
    else:
        bg = np.min(temporal, 0)

    for t in np.arange(len(temporal)):
        temporal[t, :, :] -= bg

    #temporal[0, :, :] = 0
    p = 0.0
    for y in range(np.size(temporal, 1)):
        for x in range(np.size(temporal, 2)):
            temporal[:, y, x]=filtfilt(b,a, temporal[:, y, x], padlen=padlen)
        if (100 * y / np.size(temporal, 1)) > p:
            p = (100 * y / np.size(temporal, 1))
            print("  %d%%\r" % (100 * y / np.size(temporal, 1))),
    print('Temporally Filtered')
    return temporal

    def __init__(self, channels, sample_rate, leds=None):
        self.leds = leds # Not needed if just plotting
        self.channels = channels
        self.sample_rate = sample_rate
        self.nyquist = float(sample_rate) / 2

        # Filter order - higher the order the sharper
        # the curve
        order = 3

        # Cut off frequencies:
        # Low pass filter
        cutoff = 200 / self.nyquist
        # Numerator (b) and denominator (a)
        # polynomials of the filter. 
        b, a = butter(order, cutoff, btype='lowpass')
        self.low_b = b
        self.low_a = a

        # High pass filter
        cutoff = 4000 / self.nyquist
        b, a = butter(order, cutoff, btype='highpass')
        self.high_b = b
        self.high_a = a
        
        # Keep track of max brightness for each
        # colour
        self.max = [0.0, 0.0, 0.0]
        # Make different frequencies fall faster
        # bass needs to be punchy.
        self.fall = [15.0, 2.5, 5.0]

def prepare_audio_filters():
    tf_rangel = 100000
    tf_rangeh = 170000

    # audio filters
    tf = SysParams["audio_lfreq"]
    N, Wn = sps.buttord(
        [(tf - tf_rangel) / (freq_hz / 2.0), (tf + tf_rangel) / (freq_hz / 2.0)],
        [(tf - tf_rangeh) / (freq_hz / 2.0), (tf + tf_rangeh) / (freq_hz / 2.0)],
        5,
        15,
    )
    Faudl = filtfft(sps.butter(N, Wn, btype="bandpass"))

    tf = SysParams["audio_rfreq"]
    N, Wn = sps.buttord(
        [(tf - tf_rangel) / (freq_hz / 2.0), (tf + tf_rangel) / (freq_hz / 2.0)],
        [(tf - tf_rangeh) / (freq_hz / 2.0), (tf + tf_rangeh) / (freq_hz / 2.0)],
        5,
        15,
    )
    Faudr = filtfft(sps.butter(N, Wn, btype="bandpass"))

    N, Wn = sps.buttord(0.016 / (afreq / 2.0), 0.024 / (afreq / 2.0), 5, 15)
    FiltAPost = filtfft(sps.butter(N, Wn))

    N, Wn = sps.buttord(3.1 / (freq / 2.0), 3.5 / (freq / 2.0), 1, 20)
    SysParams["fft_audiorf_lpf"] = Faudrf = filtfft(sps.butter(N, Wn, btype="lowpass"))

    SysParams["fft_audiolpf"] = FiltAPost  # * FiltAPost * FiltAPost

    SysParams["fft_audio_left"] = Faudrf * Faudl * fft_hilbert
    SysParams["fft_audio_right"] = Faudrf * Faudr * fft_hilbert

def butter_bandpass(lowcut, highcut, samplingrate, order=4):
    nyq = 0.5 * samplingrate
    low = lowcut / nyq
    high = highcut / nyq
    print high, low
    if high >=1. and low == 0.:
        b = np.array([1.])
        a = np.array([1.])

    elif high < 0.95 and low > 0. :
        wp = [1.05*low,high-0.05]
        ws = [0.95*low,high+0.05]
 
        print wp,ws
        order,wn = buttord(wp,ws,0., 30.)
        b, a = butter(order, wn, btype='band')
    
    elif high>= 0.95:
        print 'highpass',low,1.2*low,0.8*low
        order,wn = buttord( 15*low,0.05*low,gpass=0.0, gstop=10.0)
        print order,wn
        b, a = butter(order, wn, btype='high')
    elif low <= 0.05:
        print 'lowpass',high
        order,wn = buttord( high-0.05,high+0.05,gpass=0.0, gstop=10.0)
        b, a = butter(order, wn, btype='low')

    return b, a

    def __init__(self, fs, filter_type):

        self.filter_type = filter_type # Of two types. Alpha and Bandpass
        self.fs = fs # sample frequency: 220hz
        
        if self.filter_type == 'alpha':
            # butter = butterworth filter function
            # An optimal (in one way -- smoothness) for designing a filter. Simple and popular.
            # signal.butter is from scipy. Is there an equivalent in Java? Could be developed in Python and copied
            # and pasted into Java. The values of these arrays is really only dependent on sampling frequency and thus
            # likely to be the same for all uses of the app
            # Returns two arrays of floats. Could b
            self.b, self.a = signal.butter(5# order of filter. More coefficients = better filtration at cost of speed
                , np.array([8.,12.])# boundaries of filter (8-12hz for Alpha
                /(self.fs/2.) # /2 because this function needs normalized frequency
                , 'bandpass')
        elif self.filter_type == 'bandpass':
            self.b, self.a = signal.butter(5, np.array([2., 36.] # bandpass filter has different freq cutoffs
                )/(self.fs/2.), 'bandpass')
        else:
            print('Filter type ''%s'' not supported.'%self.filter_type)
            return

        self.nb = len(self.b)
        self.na = len(self.a)

def generate_noise(D,N):
  """Generate data for the changepoint detection. Data can either be of type 0
  or type 1, but when it's a combination fo both, we define a target label
  Input
  - D,N Dimenstionality arguments D dimensions over N samples
  Output
  - Data in format
  X is a matrix in R^{N x D}
  y is a matrix in R^{N,} not to donfuse with {N,1}"""
  #Check if we have even D, so we can split the array in future
  assert D%2 == 0, 'We need even number of dimensions'
  ratioP = 0.5   #balance of targets
  X = np.random.randn(N,D)
  y = np.zeros(N)
  mark = np.zeros(N)
  #Generate two filter cofficients
  filters = {}
  filters['b1'],filters['a1'] = signal.butter(4,2.0*cutoff1/fs,btype='lowpass')
  filters['b0'],filters['a0'] = signal.butter(4,2.0*cutoff0/fs,btype='lowpass')
  for i in xrange(N):
    if np.random.rand() > 0.5:	#Half of the samples will have changepoint, other half wont
      Dcut = np.random.randint(pattern_len,D-pattern_len)
      signalA = signal.filtfilt(filters['b1'],filters['a1'],X[i])
      signalB = signal.filtfilt(filters['b0'],filters['a0'],X[i])
      X[i] = np.concatenate((signalA[:Dcut],signalB[Dcut:]),axis=0)    #Concatenate the two signals
      if True:  #Boolean: do you want to introduce a pattern at the changepoint?
        Dstart = int(Dcut - pattern_len/2)
        X[i,Dstart:Dstart+pattern_len] = pattern
      y[i] = 1		#The target label
      mark[i] = Dcut
    else:
      mode = int(np.random.rand()>ratioP)
      X[i] = signal.filtfilt(filters['b'+str(mode)],filters['a'+str(mode)],X[i])
      y[i] = 0		#The target label
  return X,y,mark   

 def __init__(self, fn, channels, to_freq, montage_channels, montage, idx=1, csp_time=[0, 0.3], a_time=0.5):
     self.data = sp.signalParser(fn)
     self.tags = self.data.get_p300_tags(idx=idx, samples=False)
     signal = self.data.prep_signal(to_freq, channels, montage_channels, montage)
     signal2 = np.zeros(signal.shape)
     self.fs = to_freq
     signal2 = np.zeros(signal.shape)
     self.wrong_tags = self.data.get_p300_tags(idx=idx, rest=True, samples=False)
     artifacts_data = np.zeros([len(channels), self.fs * a_time, len(self.tags)])
     for i, p in enumerate(self.tags):
         artifacts_data[..., i] = signal[
             :, p * self.data.sampling_frequency : (p + a_time) * self.data.sampling_frequency
         ]
     self.a_features, self.bands = artifactsCalibration(artifacts_data, self.data.sampling_frequency)
     signal2 = np.zeros(signal.shape)
     self.fs = to_freq
     self.channels = channels
     self.idx = idx
     b, a = ss.butter(3, 2 * 1.0 / self.fs, btype="high")
     b_l, a_l = ss.butter(3, 20.0 * 2 / self.fs, btype="low")
     for e in xrange(len(self.channels)):
         tmp = filtfilt(b, a, signal[e, :])
         signal2[e, :] = filtfilt(b_l, a_l, tmp)
     self.signal_original = signal2
     self.t1, self.t2 = self.show_mean(csp_time, "Cz", dont_plot=False)
     P, vals = self.train_csp(signal2, [self.t1, self.t2])
     self.P = P
     self.signal = np.dot(P[:, 0], signal2)

def test_fir():
    print "Testing dsp-fir"

    #ref = ones(512)
    ref = (2.0 * random.rand(512)) - 1.0

    #test short mono fir
    writeaudio(ref)
    h = signal.firwin(21, 0.4)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    writeaudio(expected, 'expected.wav')
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(expected, readaudio(), 1e-6)

    #test long mono fir
    writeaudio(ref)
    h = signal.firwin(312, 0.4)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(expected, readaudio(), 1e-6)

    #test short stereo fir, mono coeffs
    writeaudio(transpose([ref,-ref]))
    h = signal.firwin(21, 0.4)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(transpose([expected, -expected]), readaudio(), 1e-6)

    #test long stereo fir, mono coeffs
    writeaudio(transpose([ref,-ref]))
    h = signal.firwin(312, 0.4)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(transpose([expected, -expected]), readaudio(), 1e-6)

    #test asymmetric mono fir
    writeaudio(ref)
    impulse = concatenate(([1], zeros(499)))
    b, a = signal.butter(2, 500.0/24000, 'low')
    h = signal.lfilter(b, a, impulse)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(expected, readaudio(), 1e-6)

    #test asymmetric stereo fir
    writeaudio(transpose([ref,-ref]))
    impulse = concatenate(([1], zeros(499)))
    b, a = signal.butter(2, 500.0/24000, 'low')
    h = signal.lfilter(b, a, impulse)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(transpose([expected, -expected]), readaudio(), 1e-6)

    os.remove('test_coeffs.txt')

def bandpass_filter(series, low=200.0, high=4000.0):
    """
    Filters given series with a 2nd order bandpass filter with default
    cutoff frequencies of 200 Hz and 4 kHz.

    Parameters
    ----------
    series : pandas.Series
        The data series to filter

    Returns
    -------
    filtered_series : pandas.Series
    """
    dt = series.index[1] - series.index[0]
    fs_nyquist = (1.0/dt) / 2.0
    if low < 0.1:
        # Lowpass filter only.
        bf, af = signal.butter(2, high/fs_nyquist, btype='lowpass')
    elif high > 10000:
        # Highpass filter only.
        bf, af = signal.butter(2, low/fs_nyquist, btype='highpass')
    else:
        bf, af = signal.butter(2, (low/fs_nyquist, high/fs_nyquist),
                               btype='bandpass')
    return pd.Series(signal.filtfilt(bf, af, series).astype(np.float32),
                     index=series.index)

def decode_efm():
    """ Decode EFM from STDIN, assuming it's a 28Mhz 8bit raw stream  """
    datao = np.fromstring(sys.stdin.read(SAMPLES), dtype=np.uint8).astype(np.int16)
    datao = sps.detrend(datao, type='constant')  # Remove DC

    datao = auto_gain(datao, 10000, 'pre-filter')  # Expand before filtering, since we'll lose much of signal otherwise

    low_pass = sps.butter(4, 1.75 / FREQ_MHZ, btype='lowpass')  # Low pass at 1.75 Mhz
    datao = sps.lfilter(low_pass[0], low_pass[1], datao)

    high_pass = sps.butter(4, 0.01333 / FREQ_MHZ, btype='highpass')  # High pass at 13.333 khz
    datao = sps.lfilter(high_pass[0], high_pass[1], datao)

    # This is too slow, need to work out a way to do it in scipy
    de_emphasis_filter = biquad_filter(-1.8617006585639506, 0.8706642683920058, 0.947680874725466, -1.8659578411373265, 0.9187262110931641)
    datao = np.fromiter(run_filter(de_emphasis_filter, datao), np.int16)  # De-emph - 26db below 500khz

    # Could tie edge_pll and run_filter together as generators, but we want to see the filter output

    bit_gen = edge_pll(datao, EFM_PIXEL_RATE)  # This is a ultra-naive PLL that returns a bit-stream of 1 = edge, 0 = no-edge
    try:
        while 1:
            run_until_start_code(bit_gen)
            eat_three_bits(bit_gen)
            process_efm_frame(bit_gen, 31)  # 31 14 bit EFM codes in a frame
    except StopIteration:
        printerr('Hit the end of the bitstream')

    datao = np.clip(datao, 0, 255).astype(np.uint8)
    sys.stdout.write(datao.tostring())

 def bandfilter(data, fa=None, fb=None, Fs=1000.0, order=4, zerophase=True, bandstop=False):
     N = len(data)
     assert len(shape(data)) == 1
     padded = zeros(2 * N, dtype=data.dtype)
     padded[N / 2 : N / 2 + N] = data
     padded[: N / 2] = data[N / 2 : 0 : -1]
     padded[N / 2 + N :] = data[-1 : N / 2 - 1 : -1]
     if not fa == None and not fb == None:
         if bandstop:
             b, a = butter(order, array([fa, fb]) / (0.5 * Fs), btype="bandstop")
         else:
             b, a = butter(order, array([fa, fb]) / (0.5 * Fs), btype="bandpass")
     elif not fa == None:
         # high pass
         b, a = butter(order, fa / (0.5 * Fs), btype="high")
         assert not bandstop
     elif not fb == None:
         # low pass
         b, a = butter(order, fb / (0.5 * Fs), btype="low")
         assert not bandstop
     else:
         assert 0
     if zerophase:
         return filtfilt(b, a, padded)[N / 2 : N / 2 + N]
     else:
         return lfilter(b, a, padded)[N / 2 : N / 2 + N]
     assert 0