def vad_callback(data):
	global big_data; global WINSIZE; global vad_pub; global stream; global SAMPLE_RATE
	global counter; global background_noise; global FILTER
	signal = np.array(data.data, dtype=np.int16)
	#stream.write(np.asarray(signal))
	print "recieved = " + str(len(signal)) + " frames = " + str(float(len(signal))/SAMPLE_RATE) + " seconds"
	#signal = np.asarray(big_data)
	
	window = sp.hanning(WINSIZE)
	ltsd = LTSD(WINSIZE,window,5)
	res =  ltsd.compute(signal)
	start, end = fence(res, len(signal))
	final = np.array(signal[start:end],dtype=np.float32)
	print 'start = ' + str(start)
	print 'end   = ' + str(end)
	if end - start > SAMPLE_RATE/2:
		#there is speech activity in the sample
		#print signal
		print "FOUND ACTIVITY - " + str(max(final))
		
		if FILTER and len(background_noise) > 0: #if activity is grater than half a sec:
			#take the last bg_noise in the list for better filtering
			f = cocktail(signal, background_noise[len(background_noise)-1])
			vad_pub.publish(np.array(f[0], dtype=np.float32))
		else:
			vad_pub.publish(np.array(signal, dtype=np.float32))
	else:
		if FILTER:
			background_noise.append(signal)
			if len(background_noise) > 5:
				background_noise = []
				background_noise.append(signal)

def stft(x, fs, framesz, hop):
    framesamp = int(framesz*fs)
    hopsamp = int(hop*fs)
    w = scipy.hanning(framesamp)
    X = scipy.array([scipy.fft(w*x[i:i+framesamp])
                     for i in range(0, len(x)-framesamp, hopsamp)])
    return X

def test(filename=None):
    import random, os
    import matplotlib.pyplot as plt
    from sys import argv
    #signal, params = read_signal(sound,WINSIZE)
    scenario=None
    if filename != None:
        scene = os.path.basename(filename)[0]
    else:
        filename = random.choice([x for x in os.listdir("tmp/") if os.path.splitext(x)[1] == ".flac"])
        scene = filename[0]
        filename = "tmp/"+filename
    print(filename)
    truths = vad.load_truths()
    signal,rate = speech.read_soundfile(filename)
    seconds = float(len(signal))/rate
    winsize = librosa.time_to_samples(float(WINMS)/1000, rate)[0]
    window = sp.hanning(winsize)
    ltsd = LTSD(winsize,window,5)
    res, threshold,nstart,nend =  ltsd.compute(signal)
    segments = ltsd.segments(res, threshold)
    #print(float(len(signal))/rate, librosa.core.frames_to_time(len(res), 8000, winsize/2))
    segments = librosa.core.frames_to_time(segments, rate, winsize/2)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    #ax.plot((signal/np.max(signal))*np.mean(res)+np.mean(res))
    ax.plot(np.linspace(0,seconds, len(res)), res)
    ax.plot([0, seconds], [threshold, threshold])
    vad.plot_segments(truths[scene]['combined'], segments, ax)
    n1 = float(nstart)/rate
    n2 = float(nend)/rate
    ax.vlines([n1,n2], -20,20)
    plt.show()

def stft(x, fftsize=64, overlap_pct=.5):   
    hop = int(fftsize * (1 - overlap_pct))
    w = scipy.hanning(fftsize + 1)[:-1]    
    raw = np.array([np.fft.rfft(w * x[i:i + fftsize]) for i in range(0, len(x) - fftsize, hop)])
    return raw[:, :(fftsize // 2)]

    import matplotlib.pyplot as plt

def stft(x, chunk_size, hop, w=None):
    """
    Takes the short time fourier transform of x.

    Args:
      x: samples to window and transform.
      chunk_size: size of analysis window.
      hop: hop distance between analysis windows
      w: windowing function to apply. Must be of length chunk_size

    Returns:
      STFT of x (X(t, omega)) hop size apart with windows of size chunk_size.

    Raises:
      ValueError if window w is not of size chunk_size
    """

    if not w:
        w = sp.hanning(chunk_size)
    else:
        if len(w) != chunk_size:
            raise ValueError("window w is not of the correct length {0}.".format(chunk_size))
    X = sp.array([sp.fft(w*x[i:i+chunk_size])
                     for i in range(0, len(x)-chunk_size, hop)])/np.sqrt(((chunk_size/hop)/2))
    return X

def stft(data, fs, framesize = 0.075, hopsize = 0.0625):

    # data = a numpy array containing the signal to be processed
    # fs = a scalar which is the sampling frequency of the data

    objType = type(data).__name__.strip()
       
    if objType <> "ndarray":
        raise Exception('data argument is no instance of numpy.array')

    size = len(data)
    if (size < 1):
        raise Exception('data array is empty')  

    frameSamp = int(framesize * fs)
    hopSamp = int(hopsize * fs)
    window = scipy.hanning(frameSamp)

    threshold = numpy.mean(numpy.absolute(data))*0.20
    
    X = numpy.array([numpy.absolute(scipy.fft(window * data[i : (i + frameSamp)])) for i in xrange(0, len(data) - frameSamp, hopSamp) if numpy.mean(numpy.absolute(data[i : (i + frameSamp)])) > threshold])

    # Deleting the second half of each row
    # Fourier Transform gives Hermite-symmetric result for real-valued input
    X = numpy.array([X[i][: numpy.ceil((X.shape[1] + 1.0) / 2)] for i in xrange(0, X.shape[0])])
    
    return X

def stft(x, fftsize=1024, overlap=4):
    """fftsize is in samples
    """

    hop = fftsize / overlap
    w = scipy.hanning(fftsize + 1)[:-1]  # better reconstruction with this trick +1)[:-1]
    return np.array([np.fft.rfft(w * x[i:i + fftsize]) for i in range(0, len(x) - fftsize, hop)])

def analyze_whole_waveform(waveform):
    """
    niquist_freq = framerate / 2
    precision = niquist_freq / window_size

    Want precision to be within 5% of target pitches or "5 cent".
    (+-600Hz @ 12KHz to +-10Hz @ 220Hz)

    window_size = framerate / 2 / precision

    Gives window sizes in the range of:
    - 400 Frames at 8K Frames/sec
    - 2205 Frames at 44.1K Frames/sec
    """
    desired_precision = 10  # Hz
    window_size = int(waveform.framerate / 2 / desired_precision)
    hanning_window = hanning(window_size)
    spectrum = OrderedDict()
    for start_frame in range(0, len(waveform.frames),
                             int((len(hanning_window) / 2) - 1)):
        window = zeros(len(hanning_window))
        # Do I need to add a first frame case to start with half a window to
        # match the half window at the end of stream?
        for frame in range(len(window)):
            if start_frame + frame < len(waveform.frames):
                window[frame] = (hanning_window[frame] *
                                 waveform.frames[start_frame + frame])
            else:
                window[frame] = 0
        spectrum[start_frame] = analyze_window(Waveform(window))
    return spectrum

    def fft(self, window="hanning", nfft=None):
        from numpy.fft.fftpack import fft as npfft
        from numpy.fft import fftfreq as npfftfreq
        from scipy import hamming, hanning

        sig = self.get_data()
        n = sig.shape[0]

        if window == "hamming":
            win = hamming(n)
        elif window == "hanning":
            win = hanning(n)
        elif window == "square":
            win = 1
        else:
            raise StandardError("Windows is not %s" % (window,))

        #: FFT, 折り返しこみ
        if nfft is None:
            nfft = n

        spec = npfft(sig * win, n=nfft)

        #: Freq, 折り返しこみ
        freq = npfftfreq(nfft, d=1. / self.get_fs())

        # : 折り返しを削除して返却
        se = round(nfft / 2)
        spectrum = SpectrumData(data=spec[:se], xdata=freq[:se], name=self.name)
        spectrum.set_fs(self.get_fs())

        return spectrum

def mmse_stsa(infile, outfile, noise_sum):
    signal, params = read_signal(infile, WINSIZE)
    nf = len(signal)/(WINSIZE/2) - 1
    sig_out=sp.zeros(len(signal),sp.float32)

    G = sp.ones(WINSIZE)
    prevGamma = G
    alpha = 0.98
    window = sp.hanning(WINSIZE)
    gamma15=spc.gamma(1.5)
    lambdaD = noise_sum / 5.0
    percentage = 0
    for no in xrange(nf):
        p = int(math.floor(1. * no / nf * 100))
        if (p > percentage):
            percentage = p
            print "{}%".format(p),

        y = get_frame(signal, WINSIZE, no)
        Y = sp.fft(y*window)
        Yr = sp.absolute(Y)
        Yp = sp.angle(Y)
        gamma = Yr**2/lambdaD
        xi = alpha * G**2 * prevGamma + (1-alpha)*sp.maximum(gamma-1, 0)
        prevGamma = gamma
        nu = gamma * xi / (1+xi)
        G = (gamma15 * sp.sqrt(nu) / gamma ) * sp.exp(-nu/2) * ((1+nu)*spc.i0(nu/2)+nu*spc.i1(nu/2))
        idx = sp.isnan(G) + sp.isinf(G)
        G[idx] = xi[idx] / (xi[idx] + 1)
        Yr = G * Yr
        Y = Yr * sp.exp(Yp*1j)
        y_o = sp.real(sp.ifft(Y))
        add_signal(sig_out, y_o, WINSIZE, no)
    
    write_signal(outfile, params, sig_out)

def stft(x, width):
    """Short time fourier transform of a real sequence.

    This method performs a discrete short time Fourier transform. It
    uses a sliding window to perform discrete Fourier transforms on the
    data in the Window. The results are returned in an array.

    This method uses a Hanning window on the data in the window before
    calculating the Fourier transform.

    The sliding windows are overlapping by ``width / 2``.

    Parameters
    ----------
    x : ndarray
    width: int
        the width of the sliding window in samples

    Returns
    -------
    fourier : 2d complex array
        the dimensions are time, frequency; the frequencies are evenly
        binned from 0 to f_nyquist

    See Also
    --------
    spectrum, spectrogram, scipy.hanning, scipy.fftpack.rfft

    """
    window = sp.hanning(width)
    fourier = np.array([sp.fftpack.rfft(x[i:i+width] * window) for i in range(0, len(x)-width, width//2)])
    fourier *= (2 / width)
    return fourier

def stft(x, fftsize, overlap):
    '''Computes the Short Time Fourier Transform with sensible defaults : Hanning window, window length is a power of 2 
    '''
    
    hop = fftsize // overlap
    w = scipy.hanning(fftsize+1)[:-1]      # better reconstruction with this trick +1)[:-1]  
    return numpy.array([numpy.fft.rfft(w*x[i:i+fftsize]) for i in range(0, len(x)-fftsize, hop)])

def istft(X, chunk_size, hop, w=None):
    """
    Naively inverts the short time fourier transform using an overlap and add
    method. The overlap is defined by hop

    Args:
      X: STFT windows to invert, overlap and add. 
      chunk_size: size of analysis window.
      hop: hop distance between analysis windows
      w: windowing function to apply. Must be of length chunk_size

    Returns:
      ISTFT of X using an overlap and add method. Windowing used to smooth.

    Raises:
      ValueError if window w is not of size chunk_size
    """

    if not w:
        w = sp.hanning(chunk_size)
    else:
        if len(w) != chunk_size:
            raise ValueError("window w is not of the correct length {0}.".format(chunk_size))

    x = sp.zeros(len(X) * (hop))
    i_p = 0
    for n, i in enumerate(range(0, len(x)-chunk_size, hop)):
        x[i:i+chunk_size] += w*sp.real(sp.ifft(X[n]))
    return x

def stft(x, fs, framesz = .05, hop = .025):
    framesamp = int(framesz*fs)
    hopsamp = int(hop*fs)
    w = scipy.hanning(framesamp)
    X = scipy.array([np.fft.rfft(w*x[i:i+framesamp]) 
                     for i in range(0, len(x)-framesamp, hopsamp)])
    return np.real(X)

def tfplots(data, Fs = 44100, color = 'b', fract=3):

	octbin = 100.
	FFTSIZE = 2**18

	logfact = 2**(1./octbin)
	LOGN = np.floor(np.log(Fs/2)/np.log(logfact))
	# logarithmic scale from 1 Hz to Fs/2
	logscale = np.power(logfact, np.r_[:LOGN]) 

	# creating a half hanning window
	WL = data.size
	hann = sp.hanning(WL*2)
	endwin = hann[WL:2*WL]
	tf = fft(data*endwin, FFTSIZE)

	magn = np.abs(tf[:FFTSIZE/2])
	compamp = tf[:FFTSIZE/2]

	# creating 100th octave resolution log. spaced data from the lin. spaced FFT data
	logmagn = np.empty(LOGN)
	fstep = Fs/np.float64(FFTSIZE)
	
	for k in range(logscale.size):
		start = np.round(logscale[k]/np.sqrt(logfact)/fstep)
		start = np.maximum(start,1)
		start = np.minimum(start, FFTSIZE/2)
		stop = np.round(logscale[k]*np.sqrt(logfact)/fstep)
		stop = np.maximum(stop,1)
		stop = np.minimum(stop, FFTSIZE/2)
		# averaging the power
		logmagn[k] = np.sqrt(np.mean(np.power(magn[start-1:stop],2))) 

	# creating hanning window
	# fractional octave smoothing
	HL = 2 * np.round(octbin/fract)
	hh = sp.hanning(HL)
	
	L = logmagn.size
	logmagn[L-1:L+HL] = 0

	# Smoothing the log. spaced data by convonvling with the hanning window
	tmp = fftfilt(hh, np.power(logmagn,2))
	smoothmagn = np.sqrt(tmp[HL/2:HL/2+L]/hh.sum(axis=0))

	# plotting
	plt.semilogx(logscale, 20*np.log10(smoothmagn), color)

def stft(x, fftsize=1024, overlap=4):
    hop = fftsize / overlap
    w = scipy.hanning(fftsize+1)[:-1]      # better reconstruction with this trick +1)[:-1]
    l = []
    for i in range(0, len(x)-fftsize, hop):
        v = np.fft.rfft(w*x[i:i+fftsize])
        l.append(np.abs(v)**2/np.max(np.abs(v)**2))
    return np.array(l)

def stft(x, fs, hop, width):
	"""
	Compute the Short Time Fourier Transform of 'x', with sample rate 'fs' (Hz), window width 'width' (samples), and hop length 'hop' (samples)
	Ideally, width is even (this works better with the FFT)
	"""
	window = sp.hanning(width)
	out = sp.array([ao.fft(window*x[i:i+width]) for i in range(0, len(x)-width, hop)])
	times = np.arange(width/float(2*fs), len(x)/float(fs)-width/float(2*fs), hop/float(fs))
	freqs = fs*sp.array([i/float(width) for i in range(0, width/2+1)])
	return {'stft' : out, 'times' : times, 'frequencies' : freqs}

def reduce_noise(signal, noisy_signal, winsize=2**10, window=sp.hanning(2**10)):
    """ Reduce noise """
    method = SpectralSubtraction(winsize, window)

    out = sp.zeros(len(signal), sp.float32)
    power = sig.welch(noisy_signal, window=window, return_onesided=False, scaling='spectrum')[1] * window.sum()**2
    nf = len(signal)/(winsize/2) - 1
    for no in xrange(nf):
        s = get_frame(signal, winsize, no)
        add_signal(out, method.compute_by_noise_pow(s, power), winsize, no)
    return out

def istft(X, overlap=1):   
    fftsize=(X.shape[1]-1)*2
    hop = fftsize / overlap
    w = scipy.hanning(fftsize+1)[:-1]
    x = scipy.zeros(X.shape[0]*hop)
    wsum = scipy.zeros(X.shape[0]*hop) 
    for n,i in enumerate(range(0, len(x)-fftsize, hop)): 
        x[i:i+fftsize] += scipy.real(np.fft.irfft(X[n])) * w   # overlap-add
        wsum[i:i+fftsize] += w ** 2.
    pos = wsum != 0
    x[pos] /= wsum[pos]
    return x

def stft(x, fs, framesz, hop):
    #print("STFT got", x, fs, framesz, hop)
    framesamp = int(framesz*fs)
    hopsamp = int(hop*fs)
    w = scipy.hanning(framesamp)
    def do_fft(w,x,i,framesamp):
        #print("Running FFT for ", i, framesamp)
        return fft(w*x[i:i+framesamp])
    X = scipy.array([do_fft(w,x,i,framesamp) 
                     for i in range(0, len(x)-framesamp, hopsamp)])
    #print("X SHAPE IS", len(X), len(X[0]))
    return X