def noiseFilter(data_in):

	N = int(np.ceil((4 / b)))
	if not N % 2: N += 1  # Make sure that N is odd.
	n = np.arange(N)
	 
	# Compute a low-pass filter with cutoff frequency fH.
	hlpf = np.sinc(2 * fH * (n - (N - 1) / 2.))
	hlpf *= np.blackman(N)
	hlpf = hlpf / np.sum(hlpf)
	 
	# Compute a high-pass filter with cutoff frequency fL.
	hhpf = np.sinc(2 * fL * (n - (N - 1) / 2.))
	hhpf *= np.blackman(N)
	hhpf = hhpf / np.sum(hhpf)
	hhpf = -hhpf
	hhpf[(N - 1) / 2] += 1
	 
	# Convolve both filters.
	h = np.convolve(hlpf, hhpf)
	s = np.convolve(data_in, hlpf)

	fig, ax = plt.subplots()
	ax.plot(data_in)
	plt.show()
	fig1, ax1 = plt.subplots()
	ax1.plot(s)
	plt.show()
	return s

 def applyWindow(self, window="hanning", ww=0, cf=0):
     '''
     Apply window function to frequency domain data
     cf: the frequency the window is centered over [Hz]
     ww: the window width [Hz], if ww equals 0 the window covers the full range
     '''
     self.info("Applying %s window ..." % window)
     if window == "hanning":
         if ww == 0:
             w = np.hanning(self.numfreq)
         else:
             pos = int((cf - self.lowF) / self.deltaF)
             halfwidth = int(ww / (2.0 * self.deltaF))
             w = np.zeros(self.numfreq)
             w[pos - halfwidth:pos + halfwidth] = np.hanning(2 * halfwidth)
     elif window == "hamming":
         if ww == 0:
             w = np.hamming(self.numfreq)
         else:
             pos = int((cf - self.lowF) / self.deltaF)
             halfwidth = int(ww / (2.0 * self.deltaF))
             w = np.zeros(self.numfreq)
             w[pos - halfwidth:pos + halfwidth] = np.hamming(2 * halfwidth)
     elif window == "blackman":
         if ww == 0:
             w = np.blackman(self.numfreq)
         else:
             pos = int((cf - self.lowF) / self.deltaF)
             halfwidth = int(ww / (2.0 * self.deltaF))
             w = np.zeros(self.numfreq)
             w[pos - halfwidth:pos + halfwidth] = np.blackman(2 * halfwidth)
     self.data = self.data * w
     self.done()

    def test_calculate_lanczos_kernel(self):
        """
        Tests the kernels implemented in C against their numpy counterpart.
        """
        x = np.linspace(-5, 5, 11)

        values = calculate_lanczos_kernel(x, 5, "hanning")
        np.testing.assert_allclose(
            values["only_sinc"], np.sinc(x), atol=1E-9)
        np.testing.assert_allclose(
            values["only_taper"], np.hanning(len(x)), atol=1E-9)
        np.testing.assert_allclose(
            values["full_kernel"], np.sinc(x) * np.hanning(len(x)),
            atol=1E-9)

        values = calculate_lanczos_kernel(x, 5, "blackman")
        np.testing.assert_allclose(
            values["only_sinc"], np.sinc(x), atol=1E-9)
        np.testing.assert_allclose(
            values["only_taper"], np.blackman(len(x)), atol=1E-9)
        np.testing.assert_allclose(
            values["full_kernel"], np.sinc(x) * np.blackman(len(x)),
            atol=1E-9)

        values = calculate_lanczos_kernel(x, 5, "lanczos")
        np.testing.assert_allclose(
            values["only_sinc"], np.sinc(x), atol=1E-9)
        np.testing.assert_allclose(
            values["only_taper"], np.sinc(x / 5.0), atol=1E-9)
        np.testing.assert_allclose(
            values["full_kernel"], np.sinc(x) * np.sinc(x / 5.0),
            atol=1E-9)

 def __init__(self, fl=1200, fs=240, fftl=2048):
     self.fl = fl
     self.fs = fs
     self.fftl = fftl
     winpower = np.sqrt(np.sum(np.square(np.blackman(self.fl).astype(np.float32))))
     self.window = {'ana' : np.blackman(self.fl).astype(np.float32) / winpower,
                    'syn' : winpower / 0.42 / (self.fl / self.fs)}

    def _update_mode(self):
        if self.allocation.width <= 0:
            return

        if self.fft_show:
            max_freq = (self.freq_div * self.get_ticks())
            wanted_step = 1.0 / max_freq / 2 * self._input_freq
            self.input_step = max(floor(wanted_step), 1)

            self.draw_interval = 5.0

            self.set_max_samples(
                ceil(self.allocation.width / \
                        float(self.draw_interval) * 2) * self.input_step)

            # Create the (blackman) window
            self.fft_window = blackman(
                ceil(self.allocation.width / float(self.draw_interval) * 2))

            self.draw_interval *= wanted_step / self.input_step
        else:
            # Factor is just for triggering:
            time = (self.time_div * self.get_ticks())
            if time == 0:
                return
            samples = time * self._input_freq
            self.set_max_samples(samples * self.max_samples_fact)

            self.input_step = max(ceil(samples\
                                           / (self.allocation.width / 3.0)), 1)
            self.draw_interval = self.allocation.width\
                                           / (float(samples) / self.input_step)

            self.fft_window = None

def hpf():

    """
    Ref:
        https://tomroelandts.com/articles/how-to-create-a-simple-high-pass-filter

    Usage:
        h = hpf()
        s = np.convolve(s, h)
    """

    fc = 200.0 / 16000
    b = 150.0 / 16000
    N = int(np.ceil((4 / b)))
    if not N % 2: N += 1 # N, odd
    n = np.arange(N)

    # compute a low-pass filter
    h = np.sinc(2 * fc * (n - (N - 1) / 2.))
    w = np.blackman(N)
    h = h * w
    h = h / np.sum(h)

    # spectral inversion, make high-pass filter
    h = -h
    h[(N - 1) / 2] += 1

    return h

    def __init__(self, sample_rate, fft_peak=7.0, sample_format=None, search_size=1, verbose=False, signal_width=40e3):
        self._sample_rate = sample_rate
        self._fft_size=int(math.pow(2, 1+int(math.log(self._sample_rate/1000,2)))) # fft is approx 1ms long
        self._bin_size = float(self._fft_size)/self._sample_rate * 1000 # How many ms is one fft now?
        self._verbose = verbose
        self._search_size = search_size
        self._fft_peak = fft_peak

        if sample_format == "rtl":
            self._struct_elem = numpy.uint8
            self._struct_len = numpy.dtype(self._struct_elem).itemsize * self._fft_size *2
        elif sample_format == "hackrf":
            self._struct_elem = numpy.int8
            self._struct_len = numpy.dtype(self._struct_elem).itemsize * self._fft_size *2
        elif sample_format == "float":
            self._struct_elem = numpy.complex64
            self._struct_len = numpy.dtype(self._struct_elem).itemsize * self._fft_size
        else:
            raise Exception("No sample format given")

        self._window = numpy.blackman(self._fft_size)
        self._fft_histlen=500 # How many items to keep for moving average. 5 times our signal length
        self._data_histlen=self._search_size
        self._data_postlen=5
        self._signal_maxlen=1+int(30/self._bin_size) # ~ 30 ms
        self._fft_freq = numpy.fft.fftfreq(self._fft_size)
        self._signal_width=signal_width/(self._sample_rate/self._fft_size) # Area to ignore around an already found signal in Hz
        
        if self._verbose:
            print "fft_size=%d (=> %f ms)"%(self._fft_size,self._bin_size)
            print "calculate fft once every %d block(s)"%(self._search_size)
            print "require %.1f dB"%(10*math.log(self._fft_peak,10))
            print "signal_width: %d (= %.1f Hz)"%(self._signal_width,self._signal_width*self._sample_rate/self._fft_size)

 def apply_filter(self, array):
     if len(array.shape)!=1:
         return False # need 2D data
     nx = array.shape[0]
     # Apodization function
     if self.type == 1:
         ft = np.bartlett(nx)
     elif self.type == 2:
         ft = np.blackman(nx)
     elif self.type == 3:
         ft = np.hamming(nx)
     elif self.type == 4:
         ft = np.hanning(nx)
     elif self.type == 5:
         # Blackman-Harris window
         a0 = 0.35875
         a1 = 0.48829
         a2 = 0.14128
         a3 = 0.01168
         x = np.linspace(0,1,nx)
         ft = a0 - a1*np.cos(2*np.pi*x) + a2*np.cos(4*np.pi*x) + a3*np.cos(6*np.pi*x) 
     elif self.type == 6:
         # Lanczos window
         x = np.linspace(-1,1,nx)
         ft = np.sinc(x)
     elif self.type == 7:
         x = np.linspace(-1,1,nx)
         exec('x='+self.custom)
         ft = x
     else:
         ft = np.ones(nx)
     
     array.data = array.data*ft
     return True

def smavg():
    # 权重设置
    weights_exp=np.exp(np.linspace(-1,0,N))
    weights_exp/=weights_exp.sum()
    weights_sim=np.ones(N)
    weights_sim/=weights_sim.sum()
    weights_blw=np.blackman(N)
    weights_blw/=weights_blw.sum()

    # 平滑操作
    sm_exp=np.convolve(weights_exp,data)[N-1:-N+1]
    sm_sim=np.convolve(weights_sim ,data)[N-1:-N+1]
    sm_blw=np.convolve(weights_blw ,data)[N-1:-N+1]

    #np.savetxt('result',expa)
    exp_mean=1.01*np.average(sm_exp)
    sim_mean=1.02*np.average(sm_sim)

    #sub mean
    #expa=expa-exp_mean
    #sima=sima-sim_mean

    # 绘图
    #plt.plot(t,data[N-1:],lw=1.0)
    #plt.figure()
    plt.plot(t,data[N-1:],'k',label='origin')
    #plt.plot(t,sm_exp,'r',label='exp avg') #指数移动平均
    #plt.plot(t,sm_sim,'g',label='sim avg') #简单移动平均
    plt.plot(t,sm_blw,'y',lw=2.0,label='blackman win avg')  #布莱克曼窗函数
    plt.legend(loc='upper center')

    #plt.axhline(y=exp_mean,lw=2,color='y')
    #plt.axhline(y=sim_mean,lw=2,color='b')
    plt.show()

	def blackman_window(self, show=0):
		apod = N.blackman(self.data_points)
		for i in range(2):
			self.the_result.y[i] = self.the_result.y[i]*apod
		if show == 1:
			return self.the_result
		return self

    def Execute(self,data):

       #compute power in both components
       pReal = (data['real']/2.0**self.bits)
       pImag = (data['imag']/2.0**self.bits)

       # create empty complex array to combine pReal and pImag
       data = np.empty(pReal.shape, dtype=complex)
       data.real = pReal;
       data.imag = pImag;
       
       # subtract average dc from each range
       dc = np.mean(data,axis=0)
       data = np.subtract(data[:,],dc)
       data = np.clip(data,1e-15,np.max(data))

     
       coeff = np.blackman(64)
       for i in range(0,data.shape[1]):
          data[:,i] = np.convolve(data[:,i],coeff,'same')


       # compute complex 256-point for each range cell
       dMap = np.fft.fft(data,256, axis=0)
       dMap = np.fft.fftshift(dMap,axes=(0,))

       # convert values to dB
       dMap = 20*np.log10(np.abs(dMap));
    
       return dMap

def single_taper_spectrum(data, delta, taper_name=None):
    """
    Returns the spectrum and the corresponding frequencies for data with the
    given taper.
    """
    length = len(data)
    good_length = length // 2 + 1
    # Create the frequencies.
    # XXX: This might be some kind of hack
    freq = abs(np.fft.fftfreq(length, delta)[:good_length])
    # Create the tapers.
    if taper_name == 'bartlett':
        taper = np.bartlett(length)
    elif taper_name == 'blackman':
        taper = np.blackman(length)
    elif taper_name == 'boxcar':
        taper = np.ones(length)
    elif taper_name == 'hamming':
        taper = np.hamming(length)
    elif taper_name == 'hanning':
        taper = np.hanning(length)
    elif 'kaiser' in taper_name:
        taper = np.kaiser(length, beta=14)
    # Should never happen.
    else:
        msg = 'Something went wrong.'
        raise Exception(msg)
    # Detrend the data.
    data = detrend(data)
    # Apply the taper.
    data *= taper
    spec = abs(np.fft.rfft(data)) ** 2
    return spec, freq

def get_freq2(sample, rate):
    sample_len = 22100
    frequencies = (np.arange(0, sample_len/2+1)).astype('float')*rate/sample_len
    
    
    avg_sample = np.convolve(sample, np.ones((5,)))
    avg_sample = avg_sample[2:-2]
    window = np.blackman(len(sample))
    fft_data4 = np.log(abs(np.fft.rfft(avg_sample*window, sample_len))**2)

    w = 50
    t = 3

    peaks = []    
    for j in range(0,len(fft_data4)-w):
        mx_i = np.argmax(fft_data4[j:j+w])
        strength = fft_data4[j+mx_i]-np.median(fft_data4[j:j+w])
        if mx_i == w/2 and strength > t:
           d_x = frequencies[j+mx_i]
           y0,y1,y2 = fft_data4[j+mx_i-1:j+mx_i+2]
           x1 = (y2 - y0) * .5 / (2 * y1 - y2 - y0)
           peaks.append([d_x, x1, strength])
           
    peaks = np.array(peaks)
    (freq, new_peaks) = freq_from_harmonics(peaks)
    return get_note(freq)  

    def get_wave_data(self, file_path):
        with wave.open(file_path, 'rb') as wf:
            sample_width = wf.getsampwidth()
            frame_rate = wf.getframerate()
            num_frames = wf.getnframes()

            length_in_seconds = num_frames / (frame_rate * 1.0)

            num_iterations = int(length_in_seconds * 2)
            iteration_size = int(num_frames / num_iterations)

            wave_data = []

            for fragment_num in range(0, num_iterations):
                data = wf.readframes(iteration_size)

                window_size = len(data) / sample_width
                window = np.blackman(window_size)

                fragment_data = wave.struct.unpack("%dh" % window_size, data) * window
                fft_data = abs(np.fft.rfft(fragment_data)) ** 2

                wave_data = np.hstack((wave_data, fft_data[:100]))  # take only 100 frequences

        return wave_data

 def _sweep_harmonics_responses(self, response, N, freqlist, shift=True):
     g = self.input_signal.ptp() / 2
     di = self._fft_convolve(self.sweep_inverse_signal, response / g)  # deconvolved impulse response
     imp_indices = numpy.zeros(
         N + 1, dtype=int
     )  # the indices of each impulse (1st linear, 2nd the first harm. dist. etc.)
     imp_indices[0] = len(self.sweep_inverse_signal) - 1
     for n in range(1, N + 1):
         imp_indices[n] = imp_indices[0] - int(round(math.log(n + 1) * self.sweep_rate))
     imps = []
     for n in range(N):
         dl = (imp_indices[n] - imp_indices[n + 1]) // 2
         lo = imp_indices[n] - dl  # the low limit where to cut
         hi = imp_indices[n] + dl  # the high limit where to cut
         di_w = di[lo:hi]
         win = numpy.blackman(len(di_w))
         if len(di.shape) == 2:
             win = win.reshape((len(win), 1))
         di_w *= win  # smoothed version
         w = freqlist
         if shift:
             w = w * (n + 1)
         h = scipy.signal.freqz(di_w, worN=w)[1]
         off = imp_indices[n] - lo
         h *= numpy.exp(off * 1j * freqlist)
         imps.append(h)
     return imps

def main():
    if len(sys.argv) > 1:
        wf = wave.open(sys.argv[1],'rb')
        firstTone = False
    else:
        inp = raw_input('Press enter to start record')
        while inp:
            inp = raw_input('Press enter to start record')
        rec = Recorder(channels = 1)
        with rec.open('Test.wav','wb') as recfile:
            recfile.start_recording()
            print 'Start recording...'
            inp = raw_input('Press enter to stop')
            while inp:
                inp = raw_input('Press enter to stop')
        recfile.stop_recording()
        wf = wave.open('Test.wav', 'rb')
        firstTone = True
        ScaleList = getScale(0, Tones)
        curPos = 0
    swidth = wf.getsampwidth()
    RATE = wf.getframerate()
    window = np.blackman(chunk)
    p = pyaudio.PyAudio()
    stream = p.open(format =
                p.get_format_from_width(wf.getsampwidth()),
                channels = wf.getnchannels(),
                rate = RATE,
                output = True)

    data = wf.readframes(chunk)
    values = []
    while len(data) == chunk*swidth:
        stream.write(data)
        if ( firstTone == False):
            thefreq = getHz( swidth, window, RATE, data )
            posTone = binaryNoteSearch(thefreq, Tones)
            if posTone != -1:
                ScaleList = getScale( posTone, Tones )
                firstTone = True
                curPos = 0
                print "The freq is %f Hz. Note: %s" % ( round(thefreq), Tones[posTone].NoteName )
        else:
            thefreq = getHz( swidth, window, RATE, data )
            if( thefreq > ScaleList[curPos].HighFreq or ScaleList[curPos].LowFreq > thefreq ):
                posTone = binaryNoteSearch(thefreq, ScaleList)
                if posTone != -1:
                    if posTone != 16:
                        values.append(posTone)
                    if len(values) == 2:
                        Text = getString(values)
                        sys.stdout.write(Text)
                        values = []
                        sys.stdout.flush()
                    curPos = posTone;
        data = wf.readframes(chunk)

    Text = getString(values)
    print Text
    p.terminate()

    def receiving(self, y):
        received = ''

        for x in range(0, len(y), len(generator.time_segment)):
            segment = y[x: x+len(generator.time_segment)]
            fft_segment = numpy.fft.fft(segment*numpy.blackman(len(segment)))
            fft_segment /= max(abs(fft_segment))
            freqs = numpy.fft.fftfreq(fft_segment.size, d=1 / generator.sampling_frequency)
            peak_value_index = numpy.argmax(fft_segment[0:fft_segment.size / 4])
            """
            code alphabet:
            f1 = 00
            f2 - 01
            f3 - 10
            f4 - 11
            """
            if freqs[peak_value_index] == generator.frequency:
                received += "00"
            elif freqs[peak_value_index] == generator.frequency2:
                received += "01"
            elif freqs[peak_value_index] == generator.frequency3:
                received += "10"
            elif freqs[peak_value_index] == generator.frequency4:
                received += "11"
        return received

 def doPitchDetect(self):
     self.frames = []
     # play stream and find the frequency of each chunk
     swidth = 2
     window = np.blackman(self.CHUNK)
     data = self.micStream.read(self.CHUNK)
     self.frames.append(data)
     # unpack the data and times by the hamming window
     indata = np.array(wave.struct.unpack("%dh"%(len(data)/swidth),
         data))*window
     # Take the fft and square each value
     fftData=abs(np.fft.rfft(indata))**2
     # find the maximum
     which = fftData[1:].argmax() + 1
     # use quadratic interpolation around the max
     if which != len(fftData)-1:
         y0,y1,y2 = np.log(fftData[which-1:which+2:])
         x1 = (y2 - y0) * .5 / (2 * y1 - y2 - y0)
         # find the frequency
         theFreq = (which+x1)*self.RATE/self.CHUNK/2
         if (theFreq > 300 and theFreq < 715):
             self.checkFreq(theFreq)
     else:
         theFreq = which*self.RATE/self.CHUNK/2
         if (theFreq > 300 and theFreq < 715):
             self.checkFreq(theFreq)

def plot_spectrogram(x, fs, filename,nfft):
    cf=0.0
    plt.figure(figsize=(18,20))
    fig, ax1 = plt.subplots(nrows=1)
    from matplotlib import pylab
    #wrapper = lambda n: np.kaiser(n,40)
    wrapper = lambda n : np.blackman(n)
    #wrapper = lambda n : np.hamming(n)
    #wrapper = lambda n : np.hanning(n)
    #Pxx, freqs, time, im = ax1.specgram(x, NFFT=nfft,  Fs=fs, detrend=pylab.detrend_none, noverlap=nfft/2, window=wrapper(nfft))
    Pxx, freqs, time, im = ax1.specgram(x, Fs=fs, NFFT=nfft)
    ax1.set_title('specgram spectgm'+'NFFT= %d'%nfft)
    #ax3.axis('tight')
    #labels = [item.get_text() for item in ax1.get_xticklabels()]
    ax1.set_ylabel('frequency(Hz)')
    ax1.set_xlabel('time')
    #ax1.set_ylim(40000,80000.0)
    fig.colorbar(im, ax=ax1).set_label("Amplitude (dB)")
    #ax1.set_xlim(0,7)
    #c=ax1.get_yticks().tolist()
    #print c
    #a=[str(float((i+cf)/1000.0)) for i in c ]
    #print a
    #ax1.set_yticklabels(a)
    #ax1.set_ylim(30,40)
    plt.savefig(filename+'.pdf')

 def play_(self):
     wf = wave.open('output.wav', 'rb')
     chunk = 2048
     swidth = wf.getsampwidth()
     RATE = wf.getframerate()
     window = np.blackman(chunk)
     p = pyaudio.PyAudio()
     stream = p.open(format=p.get_format_from_width(wf.getsampwidth()), channels=wf.getnchannels(), rate=RATE, output=True)
     data = wf.readframes(chunk)
     thefreq = []
     while len(data) == chunk * swidth:
         stream.write(data)
         indata = np.array(wave.struct.unpack("%dh" % (len(data) / swidth), data)) * window
         fftData = abs(np.fft.rfft(indata)) ** 2
         which = fftData[1:].argmax() + 1
         if which != len(fftData) - 1:
             y0, y1, y2 = np.log(fftData[which - 1:which + 2:])
             x1 = (y2 - y0) * .5 / (2 * y1 - y2 - y0)
             thefreq.append((which + x1) * RATE / chunk)
         else:
             thefreq.append(which * RATE / chunk)
         data = wf.readframes(chunk)
     if data:
         stream.write(data)
     stream.close()
     p.terminate()
     self.olahFrequency(thefreq)