def graph_spectrogram(wav_file):
    sound_info, frame_rate = get_wav_info(wav_file)
    pylab.figure(num=None, figsize=(19, 12))
    pylab.subplot(111)
    pylab.title('spectrogram of %r' % wav_file)
    pylab.specgram(sound_info, Fs=frame_rate)
    pylab.savefig('spectrogram - %s' % os.path.splitext(wav_file)[0])

def spectrogram(Y, Fs):

    s = Y.shape

    if len(Y.shape) > 1:
        ch = s[1]
    else:
        ch = 1

    # Ym = numpy.sum( Y , 1 ) / float(ch)

    for j in numpy.arange(ch):
        pyplot.subplot(210 + j + 1)
        pylab.specgram(Y[:, j], NFFT=Fs, Fs=Fs, cmap="gnuplot2")

        mean_x = numpy.array([0, s[0] * (1 / Fs)])
        mean_x = numpy.array([mean, mean])

        pyplot.axis([0, s[0] * (1 / Fs), 0, Fs / 2.0])

        pyplot.title("Channel %d" % j)
        pyplot.xlabel("Time [sec]")
        pyplot.ylabel("Freq. [Hz]")
        pyplot.grid(True)

        pyplot.title("Channel: %d" % int(j + 1))

    pyplot.show()

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 graph_spectrogram(wav_file):
    sound_info, frame_rate = get_wav_info(wav_file)
    pylab.figure(num=None, figsize=(19, 12))
    pylab.subplot(111)
    pylab.title('spectrogram of %r' % wav_file)
    pylab.specgram(sound_info, Fs=frame_rate)
    pylab.savefig('D:\spectrogram.png')

def spectrogram(t, x, Fs, NFFT=256):
    ax1 = pylab.subplot(211)
    pylab.plot(t, x)

    pylab.subplot(212, sharex=ax1)
    pylab.specgram(x, NFFT=NFFT, Fs=Fs, noverlap=NFFT/2,
                   cmap=pylab.cm.gist_heat)

def graph_spectrogram(wav_file):
    sound_info, frame_rate = get_wav_info(wav_file)
    pylab.figure(num=None, figsize=(19, 12))
    pylab.axes(frameon = False)
    pylab.specgram(sound_info, Fs=frame_rate)
    pylab.savefig('pngs/%s.png' % os.path.splitext(wav_file)[0], bbox_inches = 'tight')
    sleep(60)

def frequency_plot(channel_one):
    pylab.specgram(channel_one, Fs=44100)
    pylab.xlabel("Time (s)  (equivalent to distance)")
    pylab.ylabel("Frequency")
    ax=pylab.gca()
    ax.xaxis.set_major_formatter(FormatStrFormatter('%0.4f'))
    ax.yaxis.set_major_formatter(FormatStrFormatter('%0.0f'))
    pylab.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.1)

def spektrogramy():
	sred_r,sred_s=np.average(o1_bez,axis=0),np.average(o1_sty,axis=0)
	py.subplot(121)
	py.specgram(sred_r, Fs=fs, scale_by_freq=True)
	py.title('ref')
	py.subplot(122)
	py.specgram(sred_s, Fs=fs, scale_by_freq=True)
	py.title('stym')
	py.show()

def graph_spectrogram(wav_file):
    sound_info, frame_rate = get_wav_info(wav_file)
    pylab.figure(num=None, figsize=(8, 6))
    pylab.subplot(111)
    pylab.title('spectrogram of %r' % wav_file)
    pylab.ylabel('Frequency')
    pylab.xlabel('Time')
    pylab.specgram(sound_info, Fs=frame_rate)
    pylab.savefig('spectrogram.png')

def plot_spectrogram(input_data, windowfn=None, units='kHz', channel_number=0, filename=None, coloraxis=None, noverlap=0,NFFT=None, **kwargs):
    import pylab as pl
    
    if windowfn == None: windowfn=pl.window_hanning

    # look in the config file section Plots for NFFT = 1234
    # Dave - how about a method to allow this in one line
    # e.g. pyfusion.config.numgetdef('Plots','NFFT', 2048)
    # usage:  
    # if (NFFT==None): NFFT = pyfusion.config.numgetdef('Plots','NFFT', 2048)
    # 
    # also nice to have pyfusion.config.re-read()
    if NFFT == None:
        try:
            NFFT=(int(pyfusion.config.get('Plots','NFFT')))
        except:
            NFFT = 2048

    print(NFFT)        
    if units.lower() == 'khz': ffact = 1000.
    else: ffact =1.        
    xextent=(min(input_data.timebase),max(input_data.timebase))

    pl.specgram(input_data.signal.get_channel(channel_number), NFFT=NFFT, noverlap=noverlap, Fs=input_data.timebase.sample_freq/ffact, window=windowfn, xextent=xextent, **kwargs)
    #accept multi or single channel data (I think?)
        
    if coloraxis != None: pl.clim(coloraxis)
    else:
        try:
            pl.clim(eval(pyfusion.config.get('Plots','coloraxis')))
        except:
            pass

    # look in the config file section Plots for a string like 
    # FT_Axis = [0,0.08,0,500e3]   don't quote
    try:
        #pl.axis(eval(pyfusion.config.get('Plots','FT_Axis')))
        # this is clumsier now we need to consider freq units.
        axt = eval(pyfusion.config.get('Plots','FT_Axis'))
        pl.axis([axt[0], axt[1], axt[2]/ffact, axt[3]/ffact])
    except:
        pass
    # but override X if we have zoomed in bdb
    if 'reduce_time' in input_data.history:
        pl.xlim(np.min(input_data.timebase),max(input_data.timebase))
        
    try:
        pl.title("%d, %s"%(input_data.meta['shot'], input_data.channels[channel_number].name))
    except:
        pl.title("%d, %s"%(input_data.meta['shot'], input_data.channels.name))
        
    if filename != None:
        pl.savefig(filename)
    else:
        pl.show()

 def spectrogram(self, inputSignal=None, tempo=120):
     if inputSignal is None:
         inputSignal = self.data
     NFFT = int(8.0 * self.samplingRate / float(tempo))
     noverlap = NFFT >> 2
     if (len(inputSignal.shape) == 2):
         Pxx, freqs, bins, im = pylab.specgram(inputSignal[:, 0], NFFT=NFFT, Fs=self.samplingRate, noverlap=noverlap)
         return [Pxx, freqs, bins]
     else:
         Pxx, freqs, bins, im = pylab.specgram(inputSignal, NFFT=NFFT, Fs=self.samplingRate, noverlap=noverlap)
         return [Pxx, freqs, bins]

    def plot(self, noverlap=0, cmap=cm.binary, window_length=20):
        """
        Function to give Praat-style waveform + spectrogram + textgrid plot
        time is given in sec units
        """
        # plot waveform
        subplot(3, 1, 1)
        plot(self.time, self.data, color='black', linewidth=0.2)

        # plot spectrogram
        subplot(3, 1, 2)
        nfft = int(float((window_length * self.rate)) / 1000)
        specgram(self.data, NFFT=nfft, noverlap=noverlap, cmap=cmap)

 def plotPartialSpectrogram(self, timeStart, timeEnd, inputSignal=None, tempo=120):
     if inputSignal is None:
         inputSignal = self.data
     startIndex = timeStart * self.samplingRate
     stopIndex = timeEnd * self.samplingRate
     NFFT = int(8.0 * self.samplingRate / float(tempo))
     noverlap = NFFT >> 2
     if (len(inputSignal.shape) == 2):
         Pxx, freqs, bins, im = pylab.specgram(inputSignal[startIndex:stopIndex, 0], NFFT=NFFT, Fs=self.samplingRate, noverlap=noverlap)
         pylab.show()
     else:
         Pxx, freqs, bins, im = pylab.specgram(inputSignal[startIndex:stopIndex], NFFT=NFFT, Fs=self.samplingRate, noverlap=noverlap)
         pylab.show()

def call_spec():
    global y,NFFT,Fsamp,Fcentre,foverlap,detrend,_window, _type, fmod, chan_name, diag_name, hold
    print len(y), NFFT,foverlap, _type, fmod
    ax = pl.subplot(111)
    z=_window(y)
    if _type=='F': 
        shot=callback.get_shot()
        print("shot={s}".format(s=shot))
        data = device.acq.getdata(shot, diag_name)    
        if chan_name=='':
            try:
                ch=data.channels
                print("Choosing from", [chn.name for chn in ch])
                name=ch[channel_number].name
            except:
                print "Failed to open channel database - try mirnov_1_8"
                name='mirnov_1_8'
                name='mirnov_linear_2'
        else:        
            name=chan_name

#        data = pyfusion.load_channel(shot,name)
#        data = pyfusion.acq.getdata(shot_number, diag_name)    
        if type(data)==type(None): return(False)
        
        if _window==local_none: windowfn=pl.window_none
#        else: windowfn=pl.window_hanning
        elif _window==local_hanning: windowfn=pl.window_hanning
        else: windowfn=_window(arange(NFFT))
        clim=(-70,0)   # eventually make this adjustable
        if hold==0:  pl.clf()
# colorbar commented out because it keeps adding itself
        data.plot_spectrogram(NFFT=NFFT, windowfn=windowfn, noverlap=foverlap*NFFT, 
                              channel_number=channel_number)
#                         colorbar=True, clim=clim)
#        colorbar() # used to come up on a separate page, fixed, but a little clunky - leave for now

        return(True)
    elif _type == 'T':
# some matplotlib versions don't know about Fc
        pl.specgram(z*y, NFFT=NFFT, Fs=Fsamp, detrend=detrend,
#                 window = _window
                 noverlap=foverlap*NFFT, cmap=cmap)
    elif _type == 'L':
        pl.plot(20*log10(abs(fft.fft(y*z))))
    elif _type == 'W':
        pl.plot(z)
    elif _type =='C':
        pl.plot(hold=0)
    else: raise ' unknown plot type "' + _type +'"'

def graph_spectrogram(wav_file):
    sound_info, frame_rate = get_wav_info(wav_file)
    pylab.figure(num=None, figsize=(8, 6))
    pylab.subplot(111)
    pylab.title('spectrogram of %r' % wav_file)
    pylab.specgram(sound_info, Fs=frame_rate,
                   NFFT=4096, noverlap=4000,
                   )
    pylab.grid(True)
    pylab.axis((0.1, 1.85, 600, 2600))
    pylab.yticks([ 740.0,  830.6,  932.6, 1108.7, 1244.5,
                  1480.0, 1661.2, 1864.66, 2217.46, 2489.0])
    file_name, file_ext = os.path.splitext(wav_file)
    pylab.savefig('%s.%s' % (file_name, 'png'))

def read_word(word):
    try:
        data = W.read('words/%s.wav' % word)
    except:
        return
    freq = int(data[0])  # 16,000
    if not freq == 16000:
        raise Exception('shit wrong freq')

    freq = freq / DOWN_SAMPLE
    wav_data = down_sample_vect(data[1])
    if len(wav_data) > LIMIT or len(wav_data) < LIMIT * 0.9:
        return

    data = (data[0], center_vect(wav_data))
    overlap_size = int(0.01 * freq)  # 10ms ==> 160
    Pxx, freqs, bins, im = P.specgram(x=data[1], NFFT = 256, Fs=freq, 
                                    noverlap=overlap_size)
    #P.plot(data[1])
    #import pdb; pdb.set_trace()
    print Pxx.shape
    
    # For now, cut the spectrogram to make it fit
    # TODO: should use PCA instead
    if not Pxx.shape == (129, 4):
        raise Exception('shape mismatched')

    global counter
    counter += 1
    print len(data[1])

    return (data[1], Pxx)

def extract_features(fname, bdir, sox, htk_mfc, mfc_extension, stereo_wav,
        gammatones, spectrograms, filterbanks):
#def extract_features(fname, bdir):
    if fname[-4:] != '.wav':
        return
    rawfname = bdir+'/'+fname[:-4]+'.rawaudio'
    wavfname = bdir+'/'+fname
    tempfname = bdir+'/'+fname[:-4]+'_temp.wav'
    # temp fname with .wav for sox
    mfccfname = bdir+'/'+fname[:-4]+mfc_extension
    if sox:
        shutil.move(wavfname, tempfname)
        call(['sox', tempfname, wavfname])
        #call(['sox', '-G', tempfname, '-r 16k', wavfname])
        # w/o headers, sox uses extension
        shutil.move(tempfname, rawfname)
    if htk_mfc:
        call(['HCopy', '-C', 'wav_config', wavfname, mfccfname])
    srate = 16000
    #srate, sound = wavfile.read(wavfname)
    sound, srate = readwav(wavfname)
    if stereo_wav and len(sound.shape) == 2: # in mono sound is a list
        sound = 0.5 * (sound[:, 0] + sound[:, 1])
        # for stereo wav, sum both channels
    if gammatones:
        gammatonefname = bdir+'/'+fname[:-4]+'_gamma.npy'
        tmp_snd = loadsound(wavfname)
        gamma_cf = erbspace(20*Hz, 20*kHz, N_GAMMATONES_FILTERS)
        gamma_fb = Gammatone(tmp_snd, gamma_cf)
        with open(gammatonefname, 'w') as o_f:
            npsave(o_f, gamma_fb.process())
    if spectrograms:
        powerspec, _, _, _ = specgram(sound, NFFT=int(srate
            * SPECGRAM_WINDOW), Fs=srate, noverlap=int(srate
                * SPECGRAM_OVERLAP)) # TODO
        specgramfname = bdir+'/'+fname[:-4]+'_specgram.npy'
        with open(specgramfname, 'w') as o_f:
            npsave(o_f, powerspec.T)
    if filterbanks:
        # convert to Mel filterbanks
        fbanks = Spectral(nfilt=N_FBANKS,      # nb of filters in mel bank
                     alpha=0.97,               # pre-emphasis
                     do_dct=False,             # we do not want MFCCs
                     compression='log',
                     fs=srate,                 # sampling rate
                     lowerf=50,                # lower frequency
                     frate=FBANKS_RATE,        # frame rate
                     wlen=FBANKS_WINDOW,       # window length
                     nfft=1024,                # length of dft
                     do_deltas=False,          # speed
                     do_deltasdeltas=False     # acceleration
                     )
        sound /= np.abs(sound).max(axis=0)  # TODO put that as option
        fbank = fbanks.transform(sound)
        fbanksfname = bdir+'/'+fname[:-4]+'_fbanks.npy'
        with open(fbanksfname, 'w') as o_f:
            npsave(o_f, fbank)
    # TODO wavelets scattergrams / scalograms
    print "dealt with file", wavfname

def LFPPowerPerTrial_SingleBand_PerChannel_Timestamps(lfp,timestamps,Avg_Fs,channels,t_start,t_before,t_after,freq_window):
	'''
	This method computes the power in a single band defined by power_band over sliding windows in the range 
	[window_start_times-time_before,window_start_times + time_after]. This is done per trial and then a plot of 
	power in the specified band is produce with time in the x-axis and trial in the y-axis. This is done per channel.

	Main difference with LFPPowerPerTrial_SingleBand_PerChannel is that t_start is in seconds and we don't assume
	a fixed time between samples.

	Inputs:
		- lfp: lfp data arrange in an array of data x channel
		- timestamps: time stamps for lfp samples, which may occur at irregular intervals
		- channels: list of channels for which to apply this method
		- Fs: sampling rate of data in lfp_data array in Hz
		- t_start: window start times in units of s 
		- t_before: length of time in s to look at before the alignment times in window_start_times
		- t_after: length of time in s to look at after the alignment times 
		- freq_window: list defining the window of frequencies to look at, should be of the form power_band = [f_min,f_max]
	'''

	# Set up plotting variables
	# Set up matrix for plotting peak powers

	for chann in channels:
		trial_power = []
		trial_times = []
		for i, time in enumerate(t_start):
			lfp_snippet = []
			lfp_snippet = [lfp[ind,chann] for ind in range(0,len(timestamps)) if (timestamps[ind] <= time + t_after)&(time - t_before <= timestamps[ind])]
			lfp_snippet = np.array(lfp_snippet)
			
			Sxx,f,t, fig = specgram(lfp_snippet,NFFT = 128,Fs=Avg_Fs,noverlap=56,scale_by_freq=False)
			f_band_ind = [ind for ind in range(0,len(f)) if (f[ind] <= freq_window[1])&(freq_window[0] <= f[ind])]
			f_band = np.sum(Sxx[f_band_ind,:],axis=0)
			f_band_norm = f_band/np.sum(f_band)
			trial_power.append(f_band_norm)
			trial_times.append(t)
		
		power_mat = np.zeros([len(t_start),len(trial_power[0])])
		
		for ind in range(0,len(t_start)):
			if len(trial_power[ind]) == len(trial_power[0]):
				power_mat[ind,:] = trial_power[ind]
			else:
				power_mat[ind,:] = np.append(trial_power[ind],np.zeros(len(trial_power[0])-len(trial_power[ind])))
		
		dx, dy = float(t_before+t_after)/len(t), 1
		y, x = np.mgrid[slice(0,len(t_start),dy),
			slice(-t_before,t_after,dx)]
		cmap = plt.get_cmap('RdBu')
		plt.figure()
		plt.title('Channel %i - Power in band [%f,%f]' % (chann,freq_window[0],freq_window[1]))
		plt.pcolormesh(x,y,power_mat,cmap=cmap)
		plt.ylabel('Trial num')
		plt.xlabel('Time (s)')
		plt.axis([x.min(),x.max(),y.min(),y.max()])
		plt.show()
		
	return Sxx, f, t

def decoderDiagnostics(waveform=None):
    if waveform is None:
        waveform = badPacketWaveforms[-1]
    Fs_eff = Fs/upsample_factor
    ac = wifi.autocorrelate(waveform)
    ac_t = np.arange(ac.size)*16/Fs_eff
    synch = wifi.synchronize(waveform)/float(Fs_eff)
    pl.figure(2)
    pl.clf()
    pl.subplot(211)
    pl.specgram(waveform, NFFT=64, noverlap=64-1, Fc=Fc, Fs=Fs_eff, interpolation='nearest', window=lambda x:x)
    pl.xlim(0, waveform.size/Fs_eff)
    yl = pl.ylim(); pl.vlines(synch, *yl); pl.ylim(*yl)
    pl.subplot(212)
    pl.plot(ac_t, ac)
    yl = pl.ylim(); pl.vlines(synch, *yl); pl.ylim(*yl)
    pl.xlim(0, waveform.size/Fs_eff)

def graph_spectrogram(wav_file, Specgram_path, dir_path):
    try:
        sound_info, frame_rate = get_wav_info(wav_file)
    except Exception as e:
        print e
        return

    pylab.figure(num=None, figsize=(19, 12))
    pylab.subplot(111)
    pylab.title('spectrogram of %r' % os.path.basename(wav_file).strip('.WAV').strip('.wav'))
    pylab.specgram(sound_info, Fs=frame_rate, cmap=cm.Greys)

    Specgram_dir = Specgram_path + dir_path.strip('.') + '/'
    Specgram_file_name = Specgram_dir + os.path.basename(wav_file).strip('.WAV').strip('.wav')+'.png'

    if not os.path.exists(Specgram_dir):
        os.makedirs(Specgram_dir)
    pylab.savefig(Specgram_file_name)