def getDisplacements2D(self, Z=None, window=False):
     """
     Use phase correlation to find the relative displacement between
     each time step
     """
     if Z is None:
         Z = self.getNbPixelsPerFrame()/self.getNbPixelsPerSlice()/2
     shape = np.asarray(self.get2DShape())
     if window:
         ham = np.hamming(shape[1])*np.atleast_2d(np.hamming(shape[0])).T
     else:
         ham = 1.0
     displs = np.zeros((self.getNbFrames(),2))
     a = rfft2(self.get2DSlice(T=0, Z=Z)*ham)
     for t in range(1,self.getNbFrames()):
         b = rfft2(self.get2DSlice(T=t, Z=Z)*ham)
         #calculate the normalized cross-power spectrum
         #R = numexpr.evaluate(
         #    'a*complex(real(b), -imag(b)/abs(a*complex(real(b), -imag(b))))'
         #    )
         R = a*b.conj()
         Ra = np.abs(a*b.conj())
         R[Ra>0] /= Ra[Ra>0]
         r = irfft2(R)
         #Get the periodic position of the peak
         l = r.argmax()
         displs[t] = np.unravel_index(l, r.shape)
         #prepare next step
         a = b
     return np.where(displs<shape/2, displs, displs-shape)

def lcn_mauch(X, kernel=None, rho=0):
    """Apply a version of local contrast normalization (LCN), inspired by
    Mauch, Dixon (2009), "Approximate Note Transcription...".

    Parameters
    ----------
    X : np.ndarray, ndim=2
        Input representation.
    kernel : np.ndarray
        Convolution kernel (should be roughly low-pass).
    rho : scalar
        Scalar applied to the final output for heuristic range control.

    Returns
    -------
    Z : np.ndarray
        The processed output.
    """
    if kernel is None:
        dim0, dim1 = 15, 37
        dim0_weights = np.hamming(dim0 * 2 + 1)[:dim0]
        dim1_weights = np.hamming(dim1)
        kernel = dim0_weights[:, np.newaxis] * dim1_weights[np.newaxis, :]

    kernel /= kernel.sum()
    Xh = convolve2d(X, kernel, mode='same', boundary='symm')
    V = hwr(X - Xh)
    S = np.sqrt(
        convolve2d(np.power(V, 2.0), kernel, mode='same', boundary='symm'))
    S2 = np.zeros(S.shape) + S.mean()
    S2[S > S.mean()] = S[S > S.mean()]
    if S2.sum() == 0.0:
        S2 += 1.0
    return V / S2**rho

def parse_ICA_results(ICA, buffer_window): #time
	signals = {}
	signals["id"] = "ICA"
	signals["bufferWindow"] = buffer_window

	# ** for 3 channels with ICA**
	one = np.squeeze(np.asarray(ICA[:, 0])).tolist()
	two = np.squeeze(np.asarray(ICA[:, 1])).tolist()
	three = np.squeeze(np.asarray(ICA[:, 2])).tolist()
	
	one = (np.hamming(len(one)) * one)
	two = (np.hamming(len(two)) * two)
	three = (np.hamming(len(three)) * three)

	one = np.fft.irfft(one).astype(float).tolist()
	two = np.fft.irfft(two).astype(float).tolist()
	three = np.fft.irfft(three).astype(float).tolist()

	power_ratio = [0, 0, 0]
	power_ratio[0] = np.sum(one)/np.amax(one)
	power_ratio[1] = np.sum(two)/np.amax(two)
	power_ratio[2] = np.sum(three)/np.amax(three)

	if np.argmax(power_ratio) == 0:
		signals["array"] = one
	elif np.argmax(power_ratio) == 1:
		signals["array"] = two
	else:
		signals["array"] = three

	print power_ratio
	print signals
	return signals

 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 show_raw(item):
    global audio_x, audio_y, freq, curItem, filt
    curItem = item
    if l_rms is not None:
        p1.removeItem(l_rms)
        p2.removeItem(l_rms2)
    if l_fdlp is not None:
        p1.removeItem(l_fdlp)
        p2.removeItem(l_fdlp2)
    if l_tae is not None:
        p1.removeItem(l_tae)
        p2.removeItem(l_tae2)
    fn = item.text()
    w = ef.load_audio(fn)
    freq = w[2]
    p1.clear()
    p2.clear()
    x_r = np.arange(len(w[0]))/float(w[2])
    audio_x = x_r
    audio_y = w[0]/max(abs(w[0]))
    if filt is not None:
        ind = int(filt * len(audio_y)/freq)
        if ind > 1024:
            filter = np.append(np.zeros(ind-1024),np.hamming(2048),np.zeros(len(audio_y)-ind-1024))
        else:
            filter = np.append(np.hamming(2*ind),np.zeros(len(audio_y)-2*ind))
        audio_y = np.real(np.fft.ifft(np.fft.fft(audio_y)*filter))
        audio_y = audio_y/max(abs(audio_y))
    p1.plot(audio_x,audio_y,pen=(1,4))
    lr.setBounds([x_r[0],x_r[-1]])
    lr.setRegion([x_r[0],x_r[-1]])
    p1.addItem(lr)
    p2.plot(audio_x,np.abs(audio_y),pen=(1,4))

def average_energy(audio, fs=44100, n=1024):

    Ew = np.sum(np.hamming(n)**2)
    result = np.empty(len(audio)/n)
    for i in range(0,len(audio)/n):
        result[i] = (np.sum(np.absolute(np.hamming(n)*audio[i*n:(i+1)*n]))/(float(n)*Ew))
    t = np.arange(len(result)) * (float(n)/fs)      
    return result, t

def process_patch(X):
    win = np.outer(
        np.hamming(X.shape[0]), np.hamming(X.shape[1])
    )
    if np.any(np.iscomplex(X)):
        return np.abs(np.fft.fftn(X*win))**0.5
    else:
        return np.abs(np.fft.rfftn(X*win))**0.5

 def getDispl2DImage(self, t0=0, t1=1, Z=0):
     ham = np.hamming(self.get2DShape()[1])*np.atleast_2d(np.hamming(self.get2DShape()[0])).T
     a = rfft2(self.get2DSlice(T=t0, Z=Z)*ham)
     b = rfft2(self.get2DSlice(T=t1, Z=Z)*ham)
     R = numexpr.evaluate(
         'a*complex(real(b), -imag(b)/abs(a*complex(real(b), -imag(b))'
         )
     return irfft2(R)

def align(frames, template):
    """
    Warp each slice of the 3D array frames to align it to *template*.

    """
    if frames.shape[:2] != template.shape:
        raise ValueError('Template must be same shape as one slice of frame array')

    # Calculate xs and ys to sample from one frame
    xs, ys = np.meshgrid(np.arange(frames.shape[1]), np.arange(frames.shape[0]))

    # Calculate window to use in FFT convolve
    w = np.outer(np.hamming(template.shape[0]), np.hamming(template.shape[1]))

    # Calculate a normalisation for the cross-correlation
    ccnorm = 1.0 / fftconvolve(w, w)

    # Set border of normalisation to zero to avoid overfitting. Borser is set so that there
    # must be a minimum of half-frame overlap
    ccnorm[:(template.shape[0]>>1),:] = 0
    ccnorm[-(template.shape[0]>>1):,:] = 0
    ccnorm[:,:(template.shape[1]>>1)] = 0
    ccnorm[:,-(template.shape[1]>>1):] = 0

    # Normalise template
    tmpl_min = template.min()
    norm_template = template - tmpl_min
    tmpl_max = norm_template.max()
    norm_template /= tmpl_max

    warped_ims = []
    for frame_idx in xrange(frames.shape[2]):
        logging.info('Aligning frame {0}/{1}'.format(frame_idx+1, frames.shape[2]))
        frame = frames[:,:,frame_idx]

        # Normalise frame
        norm_frame = frame - tmpl_min
        norm_frame /= tmpl_max

        # Convolve template and frame
        conv_im = fftconvolve(norm_template*w, np.fliplr(np.flipud(norm_frame*w)))
        conv_im *= ccnorm

        # Find maximum location
        max_loc = np.unravel_index(conv_im.argmax(), conv_im.shape)

        # Convert location to shift
        dy = max_loc[0] - template.shape[0] + 1
        dx = max_loc[1] - template.shape[1] + 1
        logging.info('Offset computed to be ({0},{1})'.format(dx, dy))

        # Warp image
        warped_ims.append(dtcwt.sampling.sample(frame, xs-dx, ys-dy, method='bilinear'))

    return np.dstack(warped_ims)

def apply_hamming(frames, inv=False):
    """
    Computes either the hamming window or its inverse and applies
    it to a sequence of frames.

    :param frames: Frames with dimension num_frames x num_elements_per_frame
    :param inv: Indicates if the window should be inversed.
    :return:
    """
    M = frames.shape[1]
    win = np.hamming(M)**(-1) if inv else np.hamming(M)
    return frames * win

def notSoRandomWalk(shape, std=1, trendFilterLength=32, lpfLength=16):
	"""bandpass filter a random walk so that the low-frequency trend /
	drift is eliminated and the high-frequency noise is attenuated"""
	walk = randwalk(shape, std=std)
	filt = np.hamming(trendFilterLength)
	filt /= np.sum(filt)
	whichAxis = len(walk.shape) > 1 # 0 iff 1d, else 1
	# subtract baseline drift, roughly
	trend = filters.convolve1d(walk, weights=filt, axis=whichAxis, mode='reflect')
	walk -= trend
	# subtract noisey spikes
	walk = filters.convolve1d(walk, weights=np.hamming(lpfLength), axis=whichAxis, mode='reflect')
	return walk

def get_image_data(filename):
    im = pygame.image.load(filename)
    sz = im.get_size()
    im = pygame.transform.scale(im, (sz[0]/SCALE_FACTOR, sz[1]/SCALE_FACTOR))
    im2 = im.convert(8)
    a = pygame.surfarray.array2d(im2)
    hw1 = numpy.hamming(a.shape[0])
    hw2 = numpy.hamming(a.shape[1])
    a = a.transpose()
    a = a*hw1
    a = a.transpose()
    a = a*hw2
    return a

def cfrequency(data, fs, smoothie, fk):
    """
    Central frequency of a signal.

    Computes the central frequency of the given data which can be windowed or
    not. The central frequency is a measure of the frequency where the
    power is concentrated. It corresponds to the second moment of the power
    spectral density function.

    The central frequency is returned.

    :type data: :class:`~numpy.ndarray`
    :param data: Data to estimate central frequency from.
    :param fs: Sampling frequency in Hz.
    :param smoothie: Factor for smoothing the result.
    :param fk: Coefficients for calculating time derivatives
        (calculated via central difference).
    :return: **cfreq[, dcfreq]** - Central frequency, Time derivative of center
        frequency (windowed only).
    """
    nfft = util.nextpow2(data.shape[1])
    freq = np.linspace(0, fs, nfft + 1)
    freqaxis = freq[0:nfft / 2]
    cfreq = np.zeros(data.shape[0])
    if np.size(data.shape) > 1:
        i = 0
        for row in data:
            Px_wm = welch(row, np.hamming(len(row)), util.nextpow2(len(row)))
            Px = Px_wm[0:len(Px_wm) / 2]
            cfreq[i] = np.sqrt(np.sum(freqaxis ** 2 * Px) / (sum(Px)))
            i = i + 1
        cfreq = util.smooth(cfreq, smoothie)
        #cfreq_add = \
        #        np.append(np.append([cfreq[0]] * (np.size(fk) // 2), cfreq),
        #        [cfreq[np.size(cfreq) - 1]] * (np.size(fk) // 2))
        # faster alternative
        cfreq_add = np.hstack(
            ([cfreq[0]] * (np.size(fk) // 2), cfreq,
             [cfreq[np.size(cfreq) - 1]] * (np.size(fk) // 2)))
        dcfreq = signal.lfilter(fk, 1, cfreq_add)
        #dcfreq = dcfreq[np.size(fk) // 2:(np.size(dcfreq) - np.size(fk) // 2)]
        # correct start and end values of time derivative
        dcfreq = dcfreq[np.size(fk) - 1:np.size(dcfreq)]
        return cfreq, dcfreq
    else:
        Px_wm = welch(data, np.hamming(len(data)), util.nextpow2(len(data)))
        Px = Px_wm[0:len(Px_wm) / 2]
        cfreq = np.sqrt(np.sum(freqaxis ** 2 * Px) / (sum(Px)))
        return cfreq

    def __init__(self, winSecs, soundSecs, sampleRate):
        """

        :param winSecs:
        :param soundSecs:
        :param sampleRate:
        """
        self.sampleRate = sampleRate
        self.winSecs = winSecs
        self.winSamples = int(round(sampleRate*winSecs))
        self.soundSecs = soundSecs
        self.soundSamples = int(round(sampleRate*soundSecs))
        self.startWindow = numpy.hamming(self.winSamples*2)[0:self.winSamples]
        self.endWindow = numpy.hamming(self.winSamples*2)[self.winSamples:]
        self.finalWinStart = self.soundSamples-self.winSamples

def makeimg(wav):
	global callpath
	global imgpath

	fs, frames = wavfile.read(os.path.join(callpath, wav))
	
	pylab.ion()

	# generate specgram
	pylab.figure(1)
	
	# generate specgram
	pylab.specgram(
		frames,
		NFFT=256, 
		Fs=22050, 
		detrend=pylab.detrend_none,
		window=numpy.hamming(256),
		noverlap=192,
		cmap=pylab.get_cmap('Greys'))
	
	x_width = len(frames)/fs
	
	pylab.ylim([0,11025])
	pylab.xlim([0,round(x_width,3)-0.006])
	
	img_path = os.path.join(imgpath, wav.replace(".wav",".png"))

	pylab.savefig(img_path)
	
	return img_path

def iff_filter(sig, scale, plot_show = 0):
    
    order = max(sig.size*scale,90)
    #order = 80
    # Extend signal on both sides for removing boundary effect in convolution
    sig_extend = np.ones(sig.size+int(order/2)*2)
    sig_extend[int(order/2):(sig.size+int(order/2))] = sig
    sig_extend[0:int(order/2)] = sig[(sig.size-int(order/2)):sig.size]
    sig_extend[(sig.size+int(order/2)):sig_extend.size] = sig[0:int(order/2)]
    
    # convolve with hamming window and normalize
    smooth_sig = np.convolve(sig_extend,np.hamming(order),'same')
    smooth_sig = smooth_sig[int(order/2):(sig.size+int(order/2))]
    smooth_sig = np.amax(sig)/np.amax(smooth_sig)*smooth_sig

    # Plot signal for debug
    if(plot_show == 1):
        fig, ax = plt.subplots(ncols=2)
        ax[0].plot(sig)
        ax[0].plot(smooth_sig,'-r')
        ax[0].plot(med_sig,'black')
        ax[1].loglog(rfft(sig))
        ax[1].loglog(rfft(smooth_sig),'-r')
        ax[1].loglog(rfft(med_sig),'black')
        plt.show()
        
    return smooth_sig

def smooth(params, win, mode = _SMOOTH):
    
    """
    gaussian smoothing
    """
    if win >= len(params)-1:
        win = len(params)-1
    if win % 2 != 0:
        win+=1
    
    s = np.r_[params[win-1:0:-1],params,params[-1:-win:-1]]
    w = np.hamming(win)
    
        
    y = np.convolve(w/w.sum(),s,mode='valid')

    if mode == _DETREND:
        
        yy = y[(win/2-1):-(win/2)]
        return params-yy

    elif mode == _TREND:
        return y[(win/2-1):-(win/2)]
    else: 
        return y[(int(round(win/2))-1):-(int(round(win/2)))]

def compute_pitch_hps(x, Fs, dF=None, Fmin=30., Fmax=900., H=5):
    # default value for dF frequency resolution
    if dF == None:
        dF = Fs / x.size
    
    # Hamming window apodization
    x = x.copy()
    x *= np.hamming(x.size)

    # number of points in FFT to reach the resolution wanted by the user
    n_fft = np.ceil(Fs / dF)

    # DFT computation
    X = np.abs(np.fft.fft(x, n=int(n_fft)))
    
    # limiting frequency R_max computation
    R = np.floor(1 + n_fft / 2. / H)

    # computing the indices for min and max frequency
    N_min = np.ceil(Fmin / Fs * n_fft)
    N_max = np.floor(Fmax / Fs * n_fft)
    N_max = min(N_max, R)
    
    # harmonic product spectrum computation
    indices = (np.arange(N_max)[:, np.newaxis] * np.arange(1, H+1)).astype(int)
    P = np.prod(X[indices.ravel()].reshape(N_max, H), axis=1)
    ix = np.argmax(P * ((np.arange(P.size) >= N_min) & (np.arange(P.size) <= N_max)))
    return dF * ix

def window (v, func='hanning', params=None):
    """ applies a windowing function to the 3D volume v (inplace, as reference) """
    
    N = v.shape[0]
    D = v.ndim
    if any( [ d != N for d in list(v.shape) ] ) or D != 3:
        raise Exception("Error: Volume is not Cube.")
    
    def apply_seperable_window (v, w):
        v *= n.reshape(w,(-1,1,1))
        v *= n.reshape(w,(1,-1,1))
        v *= n.reshape(w,(1,1,-1))
    
    if func=="hanning":
        w = n.hanning(N)
        apply_seperable_window(v,w)
    elif func=='hamming':
        w = n.hamming(N)
        apply_seperable_window(v,w)
    elif func=='gaussian':
        raise Exception('Unimplimented')
    elif func=='circle':
        c = gencoords(N,3)
        if params==None:
            r = N/2 -1
        else:
            r = params[0]*(N/2*1)
        v *= (n.sum(c**2,1)  < ( r ** 2 ) ).reshape((N,N,N))
    elif func=='box':
        v[:,0,0] = 0.0
        v[0,:,0] = 0.0
        v[0,0,:] = 0.0
    else:
        raise Exception("Error: Window Type Not Supported")

    def make_spectrogram(self, seg_length, win_flag=True):
        """Computes the spectrogram of the wave.

        seg_length: number of samples in each segment
        win_flag: boolean, whether to apply hamming window to each segment

        returns: Spectrogram
        """
        if win_flag:
            window = np.hamming(seg_length)
        i, j = 0, seg_length
        step = seg_length / 2

        # map from time to Spectrum
        spec_map = {}

        while j < len(self.ys):
            segment = self.slice(i, j)
            if win_flag:
                segment.window(window)

            # the nominal time for this segment is the midpoint
            t = (segment.start + segment.end) / 2
            spec_map[t] = segment.make_spectrum()

            i += step
            j += step

        return Spectrogram(spec_map, seg_length)