def TakeBF(self, time, show=False, saveIm=False):
        """ Call signature example: BF = M.TakeBF(70)
            -----------------------------------------

            Return the beat frequencia of signal to the given instant (time).
            Also is considered a mean of power spectrum of ten sweeps, to
            avoid some resolution problems. """
        S_mean_K  = 0.;
        S_mean_KA = 0.;
        i1 = self.NearSweep(time);
        i2 = i1 + int(self.SD.st * 1E-6 * self.SD.rate);
        # time elapsed to next sweep in milliseconds
        T_elapse = (self.SD.st + self.SD.si) * 1E-3;
        # mean of 10 spectrograms.
        for j in range(10):
            # K band
            S, f, t = mlab.specgram(self.SD.K[i1:i2], NFFT=self.nfft,
                    Fs=self.SD.rate, noverlap=self.nfft-self.fft_step, 
                    pad_to=2**12);
            S_mean_K += S / 10;
            # Ka band
            S, f, t = mlab.specgram(self.SD.KA[i1:i2], NFFT=self.nfft,
                    Fs=self.SD.rate, noverlap=self.nfft-self.fft_step, 
                    pad_to=2**12);
            S_mean_KA += S / 10;
            # update indexes to next sweep.
            i1 = self.NearSweep(time + T_elapse);
            i2 = i1 + int(self.SD.st * 1E-6 * self.SD.rate);
            T_elapse += (self.SD.st + self.SD.si) * 1E-3;
        # avoid useles frequency depends on the case
        if time < 10: 
            limSup = np.where(f > 1.65E7)[0].min();
            limInf = np.where(f < 0.65E7)[0].max();
        else:
            limSup = np.where(f > 1.55E7)[0].min();
            limInf = np.where(f < 0.35E7)[0].max();
        # Take max line in spectrum.
        freqIndexK  = S_mean_K[limInf:limSup].argmax(axis=0) + limInf;
        freqIndexKA = S_mean_KA[limInf:limSup].argmax(axis=0) + limInf;
        # return to default image inferior limite
        limSup = np.where(f > 1.55E7)[0].min();
        limInf = np.where(f < 0.35E7)[0].max();
        # Take frequency of max line.
        BF_K  = f[freqIndexK];
        BF_KA = f[freqIndexKA];
        # remove overlap if it has.
        BF = np.concatenate([BF_K, BF_KA]);
        if (saveIm or show):
            S_K  = S_mean_K[limInf:limSup];
            S_KA = S_mean_KA[limInf:limSup];
            f1 = f[limInf]; # Inferior figure limit
            f2 = f[limSup]; # superior figure limit
            self.__DrawImage(S_K, S_KA, time, f1, f2, BF, show, saveIm);
        return BF;

    def test_specgram(self):
        for y, fstims in zip(self.y, self.fstimsall):
            Pxx1, freqs1, t1 = mlab.specgram(y, NFFT=self.NFFT,
                                             Fs=self.Fs,
                                             noverlap=self.noverlap,
                                             pad_to=self.pad_to,
                                             sides='default')
            Pxx1m = np.mean(Pxx1, axis=1)
            np.testing.assert_array_equal(freqs1, self.freqss)
            np.testing.assert_array_equal(t1, self.t)
            # since we are using a single freq, all time slices should be
            # about the same
            np.testing.assert_allclose(np.diff(Pxx1, axis=1).max(), 0,
                                       atol=1e-08)
            for fstim in fstims:
                i = np.abs(freqs1 - fstim).argmin()
                self.assertTrue(Pxx1m[i] > Pxx1m[i+1])
                self.assertTrue(Pxx1m[i] > Pxx1m[i-1])

            Pxx2, freqs2, t2 = mlab.specgram(y, NFFT=self.NFFT,
                                             Fs=self.Fs,
                                             noverlap=self.noverlap,
                                             pad_to=self.pad_to,
                                             sides='onesided')
            Pxx2m = np.mean(Pxx2, axis=1)
            np.testing.assert_array_equal(freqs2, self.freqss)
            np.testing.assert_array_equal(t2, self.t)
            np.testing.assert_allclose(np.diff(Pxx2, axis=1).max(), 0,
                                       atol=1e-08)
            for fstim in fstims:
                i = np.abs(freqs2 - fstim).argmin()
                self.assertTrue(Pxx2m[i] > Pxx2m[i+1])
                self.assertTrue(Pxx2m[i] > Pxx2m[i-1])

            Pxx3, freqs3, t3 = mlab.specgram(y, NFFT=self.NFFT,
                                             Fs=self.Fs,
                                             noverlap=self.noverlap,
                                             pad_to=self.pad_to,
                                             sides='twosided')
            Pxx3m = np.mean(Pxx3, axis=1)
            np.testing.assert_array_equal(freqs3, self.freqsd)
            np.testing.assert_array_equal(t3, self.t)
            np.testing.assert_allclose(np.diff(Pxx3, axis=1).max(), 0,
                                       atol=1e-08)
            for fstim in fstims:
                i = np.abs(freqs3 - fstim).argmin()
                self.assertTrue(Pxx3m[i] > Pxx3m[i+1])
                self.assertTrue(Pxx3m[i] > Pxx3m[i-1])

def fingerprint(
    channel_samples,
    Fs=DEFAULT_FS,
    wsize=DEFAULT_WINDOW_SIZE,
    wratio=DEFAULT_OVERLAP_RATIO,
    fan_value=DEFAULT_FAN_VALUE,
    amp_min=DEFAULT_AMP_MIN,
    plot=False,
):
    """
    FFT the channel, log transform output, find local maxima, then return
    locally sensitive hashes.
    """
    # FFT the signal and extract frequency components
    # http://matplotlib.org/api/mlab_api.html#matplotlib.mlab.specgram
    # return of specgram is (spectrum, freqs, t)
    arr2D, freqs, times = mlab.specgram(
        channel_samples, NFFT=wsize, Fs=Fs, window=mlab.window_hanning, noverlap=int(wsize * wratio)
    )

    # apply log transform since specgram() returns linear array
    arr2D = 10 * np.log10(arr2D)
    arr2D[arr2D == -np.inf] = 0  # replace infs with zeros

    # find local maxima
    local_maxima = get_2D_peaks(arr2D, plot=plot, amp_min=amp_min, freqs=freqs, times=times)

    # return hashes
    return generate_hashes(local_maxima, fan_value=fan_value)

 def calculate_specgram(self, nfft, noverlap, **kwargs):
     spec = mlab.specgram(self.audio, NFFT=nfft, Fs=self.framerate, noverlap=noverlap, **kwargs)
     self.specgram_nfft = nfft
     self.specgram_overlap = noverlap
     self.specgram = spec[0]
     self.specgram_freqs = spec[1]
     self.specgram_bins = spec[2]

def fingerprint(channel_samples, Fs=DEFAULT_FS, wsize=DEFAULT_WINDOW_SIZE, wratio=DEFAULT_OVERLAP_RATIO, fan_value=DEFAULT_FAN_VALUE, amp_min=DEFAULT_AMP_MIN):
    """
    FFT the channel, loklsg transform output, find local maxima, then return
    locally sensitive hashes.
    """


    # FFT the signal and extract frequency components
    arr2D = mlab.specgram(
        channel_samples,
        NFFT=wsize,
        Fs=Fs,
        window=mlab.window_hanning,
        noverlap=int(wsize * wratio))[0]

    # apply log transform since specgram() returns linear array
    np.seterr(all='ignore')
    arr2D = 10 * np.log10(arr2D)
    arr2D[arr2D == -np.inf] = 0  # replace infs with zeros

    # find local maxima
    local_maxima = get_2D_peaks(arr2D, plot=False, amp_min=amp_min)

    # return hashes
    return generate_hashes(local_maxima, fan_value=fan_value)

def transform_one(row):
    row -= MIN
    row /= MAX
    spec = mlab.specgram(row, NFFT=256, Fs=16384)[0]
    spec = np.log(spec)#.ravel()
    #spec = ndimage.gaussian_filter(spec, 1)
    return spec

def psd_eeg(data, **kwargs):

    for arg in kwargs:
        if arg == "NFFT":
            NFFT = int(kwargs[arg])
        if arg == "dt":
            Fs = 1.0 / float(kwargs[arg])
        if arg == "noverlap":
            noverlap = int(kwargs[arg])

    px_list = []

    print "PSD Computing..."
    for sample in range(data.shape[0]):
        for ch in range(data.shape[1]):
            eeg = data[sample, ch, :]
            [Pxx, freq, t] = specgram(eeg, NFFT=NFFT, noverlap=noverlap, Fs=Fs)
            px_list.append(Pxx)

    shape = px_list[0].shape
    pot = np.array(px_list)

    pot = pot.reshape(data.shape[0], data.shape[1], shape[0], -1)

    del px_list

    return [pot, freq]

def plot_specgram(ax, data, fs, nfft=256, noverlap=128, window='hann',
                  cmap='jet', interpolation='bilinear', rasterized=True):

    if window not in SPECGRAM_WINDOWS:
        raise ValueError("Window not supported")

    elif window == "boxcar":
        mwindow = signal.boxcar(nfft)
    elif window == "hamming":
        mwindow = signal.hamming(nfft)
    elif window == "hann":
        mwindow = signal.hann(nfft)
    elif window == "bartlett":
        mwindow = signal.bartlett(nfft)
    elif window == "blackman":
        mwindow = signal.blackman(nfft)
    elif window == "blackmanharris":
        mwindow = signal.blackmanharris(nfft)

    specgram, freqs, time = mlab.specgram(data, NFFT=nfft, Fs=fs,
                                          window=mwindow,
                                          noverlap=noverlap)
    specgram = 10 * np.log10(specgram[1:, :])
    specgram = np.flipud(specgram)

    freqs = freqs[1:]
    halfbin_time = (time[1] - time[0]) / 2.0
    halfbin_freq = (freqs[1] - freqs[0]) / 2.0
    extent = (time[0] - halfbin_time, time[-1] + halfbin_time,
              freqs[0] - halfbin_freq, freqs[-1] + halfbin_freq)

    ax.imshow(specgram, cmap=cmap, interpolation=interpolation,
                            extent=extent, rasterized=rasterized)
    ax.axis('tight')

    def _energy_func(self, x, **kwargs):
        from matplotlib.mlab import specgram

        rval = sp.zeros_like(x)
        ns, nc = x.shape
        ov_samples = 0
        offset = 0
        if self.overlap == 1:
            ov_samples = self.nfft * 0.5
            offset = self.nfft / 4
        elif self.overlap == 2:
            ov_samples = self.nfft * 0.75
            offset = self.nfft * 0.375
        step = self.nfft - ov_samples

        for c in xrange(nc):
            psd_arr, freqs, times = specgram(x[:, c], NFFT=self.nfft, Fs=self.srate, noverlap=ov_samples)
            mask = freqs < self.cutoff_hz
            for b in xrange(len(times)):
                bin_s = b * step + offset
                bin_e = bin_s + step
                
                if self.en_func == 'mean_coeff':
                    rval[bin_s:bin_e, c] = psd_arr[mask == True, b].mean() / psd_arr[mask == False, b].mean()
                elif self.en_func == 'max_coeff':
                    rval[bin_s:bin_e, c] = psd_arr[mask == True, b].max() / psd_arr[mask == False, b].max()
                elif self.en_func == 'max_normed':
                    rval[bin_s:bin_e, c] = psd_arr[mask == True, b].max() / psd_arr[:, b].sum(axis = 0)
                else:
                    raise RuntimeError('Energy function does not exist!')

        return rval

def get_fingerprint_from_data(wave_data):
    # pxx[freq_idx][t] - мощность сигнала
    pxx, _, _ = mlab.specgram(
        wave_data,
        NFFT=cf.WINDOW_SIZE,
        noverlap=cf.WINDOW_OVERLAP,
        Fs=cf.SAMPLE_RATE)

    # 300-2870 | delta = 256 * 10 = 8 * 32 * 10
    matrix = pxx[30*2:300*2].transpose()

    cnt1, cnt2 = 0, 0
    arr = []
    for time, timeline in enumerate(matrix):
        if time == 0:
            continue

        hash64, pow2 = 0, 1
        for j in range(1, 65):
            energy1 = energy(matrix, time, j) - energy(matrix, time, j - 1)
            energy2 = energy(matrix, time - 1, j) - energy(matrix, time - 1, j - 1)
            if energy1 - energy2 > 0:
                hash64 += pow2
                cnt1 += 1
            else:
                cnt2 += 1
            pow2 *= 2
        arr.append(hash64)

    print('Done fingerprinting...', cnt1, cnt2)

    return [arr]

    def _plot_spectrum(self, data, sampling_frequency):
        (Pxx, freqs, bins) = mlab.specgram(data,
                                           Fs=sampling_frequency,
                                           NFFT=self.NFFT,
                                           noverlap=self.noverlap)
                    
        if numpy.any(Pxx[0,0] == 0):
            self._log("SpectrumPlot::Instance has power 0 in a frequency band, skipping...")
            return

        # Update minimal and maximal value
        self.min_value = min(self.min_value, min(Pxx.flatten()))
        self.max_value = max(self.max_value, max(Pxx.flatten()))
        
        Z = numpy.flipud(Pxx)
        extent = 0, numpy.amax(bins), freqs[0], freqs[-1]
        
        pylab.imshow(Z, None, extent=extent, vmin=self.min_value,
                     vmax=self.max_value)
        pylab.axis('auto')
        pylab.xlabel("Time(s)")
        pylab.ylabel("Frequency(Hz)")
        
        if self.colorbar and not self.physiological_arrangement:
            pylab.colorbar()
        
        return (Pxx, freqs, bins)

def simple_spectogram(df, meta):
    import numpy as np
    import pandas as pd
    from matplotlib.mlab import specgram
    
    
    def merge_two_dicts(x, y):
        z = x.copy()
        z.update(y)
        return z
    
    # extract numpy array from pandas.DataFrame
    data = np.asarray(df)
    
    spectrum, freqs, t = specgram(data.squeeze())    
    
    # Return value must be of a sequence of pandas.DataFrame
    df = pd.DataFrame(spectrum)
    

    # updata meta
    new_meta = {'spectogram_out': {'frequencies': freqs.tolist(),
                                   'time_points': t.tolist()}}
    
    d = merge_two_dicts(meta.to_dict(), new_meta)
    meta = pd.DataFrame(d)   
    
    return df, meta

def load_spec(path):
	'''
	Loads the spectrogram file that was saved by the function sess_spectrogram
	Input:
		path : either the absolute path for the spectrogram file or a list containing the
				cage ID, block number, and sess
	Output:
		P : MxN spectrogram array (M frequencies, N time bins)
		F : vector of length M indicating the frequencies included in the spectrogram
	'''

	exten = os.path.splitext(path)[-1]
	
	if exten == '.npz':
		tmp = np.load(path)
		P = np.array(tmp['P'].min().todense())
		F = tmp['F']
		T = tmp['T']
		tmp.close()
		
	elif exten == '.h5':
		NFFT = 512
		noverlap = 256

		f = h5py.File(path, 'r+')
		s = f['s'].value
		fs = f['fs'].value
		f.close()

		(P, F, T) = specgram(s, Fs = fs, NFFT = NFFT, noverlap = noverlap)

	return P, F, T

def get_fingerprints(X, 
                    ws=WINDOW_SIZE,
                    sr=SAMPLE_RATE,
                    si=SLIDE_INTERVAL,
                    ns=NEIGHBOR_SIZE,
                    thres=MIN_AMPLITUDE,
                    fo=FAN_OUT):
    """
    Calculate the fingerprints (hash values) of a given signal X.

    Return the fingerprints in term of (sha256, diff time).
    """

    # apply fft to X and retrieve the spectrogram of X
    sg, freqs, t = specgram(X, 
                            NFFT=ws,
                            Fs=sr,
                            noverlap=si)

    # apply log transform to the spectrogram
    sg = 10 * log10(sg)

    # find the local maxima in neighborhood
    peaks = find_peaks(sg, ns, thres)

    # generate final fingerprints in form of (sha256, diff time)
    fingerprints = generate_fingerprints(peaks, fo)

    return fingerprints

    def _do_fft(self):
        t = self.getp('t')
        v = self.getp('v')
        try:
            Fs = 1/(t[1] - t[0]) 
            NFFT = self.getp('NFFT')
            detrend = detrends[self.getp('detrend')]
            window = windows[self.getp('window')]
            noverlap = self.getp('noverlap')
            sides = self.getp('sides')

            pad_to = self.getp('pad_to')
            if pad_to == 'none': pad_to = None
            z, y, x = mlab.specgram(v, NFFT=NFFT, 
                               Fs=Fs, 
                               detrend = detrend,
                               window = window,
                               noverlap = noverlap,
                               pad_to=pad_to, 
                               sides = sides, 
                               scale_by_freq=None)
            self.setp("x", x + t[0])
            self.setp("y", y)
            self.setp("z", z)
        except:
            traceback.print_exc()
            self.setp("x", None)
            self.setp("y", None)
            self.setp("z", None)
            return False
        return True

def calc_power_spectrum(lfp):
	
	npts, ntrials = lfp.shape
	Fs = 384.
	X = np.zeros()
	x = np.zeros((129, 7, ntrials))
	for i in range(ntrials):
		x[:, :, i] = specgram(lfp[:127, i], Fs = Fs)[0]

def specgram(signal, sampling_frequency, time_resolution, 
             frequency_resolution, bath_signals=[], 
             high_frequency_cutoff=None,  axes=None, logscale=True, **kwargs):
    """
    This function wraps matplotlib.mlab.specgram to provide a more intuitive 
        interface.
    Inputs:
        signal                  : the input signal (a one dimensional array)
        sampling_frequency      : the sampling frequency of signal
        time_resolution         : the desired time resolution of the specgram
                                    this is the guaranteed worst time resolution
        frequency_resolution    : the desired frequency resolution of the 
                                    specgram.  this is the guaranteed worst
                                    frequency resolution.
        --keyword arguments--
        bath_signals            : Subtracts a bath signal from the spectrogram
        axes=None               : If an Axes instance is passed then it will
                                  plot to that.
        **kwargs                : Arguments passed on to 
                                   matplotlib.mlab.specgram
    Returns:
        If axes is None:
            Pxx
            freqs
            bins
        if axes is an Axes instance:
            Pxx, freqs, bins, and im
    """
    if (high_frequency_cutoff is not None 
        and high_frequency_cutoff < sampling_frequency):
        resampled_signal = resample_signal(signal, sampling_frequency, 
                                                    high_frequency_cutoff)
    else:
        high_frequency_cutoff = sampling_frequency
        resampled_signal = signal
    num_data_samples = len(resampled_signal)
    specgram_settings = find_NFFT_and_noverlap(frequency_resolution, 
                                               high_frequency_cutoff, 
                                               time_resolution, 
                                               num_data_samples)
    NFFT     = specgram_settings['power_of_two_NFFT']
    noverlap = specgram_settings['noverlap']
    Pxx, freqs, bins = mlab.specgram(resampled_signal, 
                                                NFFT=NFFT, 
                                                Fs=high_frequency_cutoff, 
                                                noverlap=noverlap, **kwargs)
    plotted_Pxx = Pxx
    if bath_signals:
        bath_signal = numpy.hstack(bath_signals)
        psd_Pxx, psd_freqs = psd(bath_signal, sampling_frequency, 
                                 frequency_resolution,
                                 high_frequency_cutoff=high_frequency_cutoff ) 
        plotted_Pxx = (Pxx.T/psd_Pxx).T

    if axes is not None:
        im = plot_specgram(plotted_Pxx, freqs, bins, axes, logscale=logscale)
        return plotted_Pxx, freqs, bins, im
    return plotted_Pxx, freqs, bins

 def score(self, audio):
     nfft = int(self.window*audio.framerate)
     audio.calculate_specgram(nfft=nfft, noverlap=nfft/2)
     freqs = np.where((audio.specgram_freqs >= self.lower_call_frequency)*(audio.specgram_freqs <= self.upper_call_frequency))
     spec2 = mlab.specgram(mean(log(audio.specgram[freqs[0],]), 0), NFFT=1024, noverlap=512, Fs=2/self.window)
     freqs2 = np.where((spec2[1] >= self.lower_syllable_frequency)*(spec2[1] <= self.upper_syllable_frequency))
     max_kiwi = max(np.max(spec2[0][freqs2[0], :], 0))
     mean_kiwi = np.exp(np.mean(np.mean(np.log(spec2[0][freqs2[0], :]), 0)))
     return max_kiwi/mean_kiwi

def generateSpectogram(channel_samples, Fs=DEFAULT_FS, wsize=DEFAULT_WINDOW_SIZE, wratio=DEFAULT_OVERLAP_RATIO, fan_value=DEFAULT_FAN_VALUE, amp_min=DEFAULT_AMP_MIN):
	# a tuple (Pxx, freqs, t) is said to be the return value of specgram
	arr2D = mlab.specgram(channel_samples, NFFT=wsize, Fs=Fs, window=mlab.window_hanning, noverlap=int(wsize * wratio))[0]

	# apply log transform since specgram() returns linear array
	arr2D = 10 * np.log10(arr2D)
	arr2D[arr2D == -np.inf] = 0 # replace infs with zeros

	# find local maxima
	plotPeaks(arr2D, amp_min=amp_min)

 def get_spectrum_data(self):
     print("Calculating spectrum data...")
     self.spec_PSDperHz, self.spec_freqs, self.spec_t = mlab.specgram(np.squeeze(self.data),
                                    NFFT=self.NFFT,
                                    window=mlab.window_hanning,
                                    Fs=self.fs_Hz,
                                    noverlap=self.overlap
                                    ) # returns PSD power per Hz
     # convert the units of the spectral data
     self.spec_PSDperBin = self.spec_PSDperHz * self.fs_Hz / float(self.NFFT)