def test_frequency(self):
        import scipy.signal as ss
        nsteps = 10000
        dt = 0.1

        t,w = self.regular_integrator.run(self.regular_w0, dt=dt, nsteps=nsteps)
        f,fft = fft_orbit(t, w)

        peak_ix = ss.find_peaks_cwt(fft[:,0], widths=np.linspace(dt*2, dt*100, 10))
        print(peak_ix)

        plt.clf()
        plt.axvline(self.regular_par[1]/(2*np.pi), linewidth=3., alpha=0.35, color='b')
        plt.axvline(1/(2*np.pi), linewidth=3., alpha=0.35, color='r')
        plt.semilogx(f[:,0], fft[:,0], marker=None)
        plt.savefig(os.path.join(plot_path,"pend_fft_regular.png"))

        # ----------------------------------------------------------------------
        t,w = self.chaotic_integrator.run(self.chaotic_w0, dt=dt, nsteps=nsteps)
        f,fft = fft_orbit(t, w)

        peak_ix = ss.find_peaks_cwt(fft[:,0], widths=np.linspace(dt*2, dt*100, 10))
        print(peak_ix)

        plt.clf()
        plt.axvline(self.chaotic_par[1]/(2*np.pi), linewidth=3., alpha=0.35, color='b')
        plt.axvline(1/(2*np.pi), linewidth=3., alpha=0.35, color='r')
        plt.semilogx(f[:,0], fft[:,0], marker=None)
        plt.savefig(os.path.join(plot_path,"pend_fft_chaotic.png"))

  def fft_analysis(self,position=None):
    '''
    Performes Fast fourier analysis on data sets and returns either five most significant frequencies or numpy arrays of 
    frequencies and power spectrum when specific position number is passed as parameter.
    '''

    if position is None:
      freqs_5=np.array([])
      
      for dp, pr in zip(self.data_points, self.prumery):
        time_step=(dp.data_time[-1]-dp.data_time[0])/dp.data.size
        ps=np.abs(np.fft.rfft(dp.data-pr))**2 #Power spectrum
      
        freqs = np.fft.rfftfreq(dp.data.size, time_step)
        peakind = signal.find_peaks_cwt(ps, np.arange(1,50))
        freqs_5 = np.append(freqs_5,freqs[peakind[:5]])
        
      return freqs[peakind[:5]] #Returns first five most significant frequencies
    
    else: #Returns arrays of frequencies and power spectrum and five most significant frequencies.
      for dp,pr in zip(self.data_points, self.prumery):
        if position == dp.position_no:
          time_step=(dp.data_time[-1]-dp.data_time[0])/dp.data.size
          ps=np.abs(np.fft.rfft(dp.data-pr))**2 #Power spectrum
          #ps=np.angle(np.fft.rfft(dp.data-pr)) #Phase spectrum
          #ps=20*np.log10(np.abs(np.fft.rfft(dp.data-pr))) #Amplitude spectrum in [dB]
        
          freqs = np.fft.rfftfreq(dp.data.size, time_step)
          peakind = signal.find_peaks_cwt(ps, np.arange(1,50))
          
          
          return freqs, ps, freqs[peakind[:5]]

def find_threshold(filename, base_offset=0, pos_pol=True, **kws):
    with trace_gen(filename, base_offset=base_offset, **kws) as gen:
        if pos_pol:
            min_data = [event[100:-100].max() for event in gen]
        else:
            min_data = [event[100:-100].min() for event in gen]
    hist = histogram(min_data, bins=2048)
    # second = 2 if pos_pol else -2
    # index = argrelmax(hist[0][hist[0] != 0], order=5)[0][second]
    second = 1 if pos_pol else -1
    try:
        index = find_peaks_cwt(hist[0][hist[0] != 0], arange(10, 40))[second]
    except IndexError:
        index = find_peaks_cwt(hist[0][hist[0] != 0], arange(10, 40))[0]
    return hist[1][hist[0] != 0][index] * 3 / 4, len(min_data)

def find_peaks(x, y, widthrange, rel_threshold=0.1):
    """Peak-finding in a 2d dataset.
    
    Parameters
    ----------
    x, y : array_like
        Input arrays.
    widthrange : tuple
        Lower and upper limit of peak widths to find.
    rel_threshold : float, optional
        Peaks with a height lower than this value times the height of the
        maximum in 'y' are ignored.

    Returns
    -------
    list
        Array indices of where the peaks were found in 'y'.

    See Also
    --------
    scipy.signal.find_peaks_cwt : Peak-finding using a continous wavelet
    transform technique.
    """
    dx = abs(x[1] - x[0])
    minwidth, maxwidth = widthrange
    widths = np.arange(floor(minwidth/dx), ceil(maxwidth/dx))
    peakpos = find_peaks_cwt(y, widths)
    maxy = max(y)
    return [pos for pos in peakpos if y[pos] >= rel_threshold*maxy]

def main(data):
    start = time.time()
    T = 1.0/SAMPLING_RATE # sampling interval
    Fs = 1.0 / T
    ir_data, bpm_data, avg_bpm_data = clean_data(data)
    hr_data = zero_mean(ir_data)

    cutoff_bpm = [50.0, 200.0]
    cutoff_hz = [x/60 for x in cutoff_bpm] # cutoff frequency in HZ
    cutoff = [x/(Fs/2) for x in cutoff_hz]
    [b, a] = signal.butter(2, cutoff, 'bandpass') # 2nd order butterworth filter

    is_valid = check_data_good(ir_data)
    if (is_valid == False):
        return "Heart Rate: Please Place Sensor on Feet"
    else:
        # filter out the noise from signal
        hr_filt = signal.lfilter(b, a, hr_data)

        pks = signal.find_peaks_cwt(hr_filt, np.arange(3, 10))
        num_pks = len(pks)
        beats_from_peaks = num_pks/2
        bpm_from_peaks = beats_from_peaks*60/TIME_SEC
        print("HR Found from Beats is = " + str(bpm_from_peaks) + " BPM")
        time_btw_peaks = sum(np.array(pks[1:num_pks]) - np.array(pks[0:-1]))/(num_pks - 1)
        bpm_from_peaks = SAMPLING_RATE*60/time_btw_peaks/2
        print("HR Found from Time is = " + str(bpm_from_peaks) + " BPM")
        end = time.time()
        print("total time = " + str(end - start))
        return bpm_from_peaks

def find_peak_data(data, n, spec=False):
    indices = find_peaks_cwt(data, arange(30, 80))
    values = [data[i] for i in indices]
    params = [val for i in range(n) for val in [values[i], indices[i], 30]]
    if spec:
        params.extend((0.2, 1))
    return params

    def execute_analysis(self):
        """
        The continuous wavelet transform peak identification algorithm from scipy.signal.find_peaks_cwt().
        In openmsi it is renamed "findpeaks_cwt" with one fewer underscore.
        """

        from scipy import signal as sig

        msidata = self['msidata']    #now in memory as hdf5 cube or np.ndarray
        widths = self['widths']      #should already be numpy ndarray
        min_snr = self['min_snr']

        shape_x = msidata.shape[0]
        shape_y = msidata.shape[1]

        mzindices = []

        for xi in xrange(0, shape_x):
            for yi in xrange(0, shape_y):
                print xi, yi

                # Load the spectrum
                m = msidata[xi, yi, :]

                # find indices of m where peaks are
                peak_indices = sig.find_peaks_cwt(m, widths=widths, wavelet=sig.ricker, min_snr=min_snr)
                mzindices.append(peak_indices)
        return mzindices

 def get_peaks(self, data, domain):
     thresh = 1.1  # TODO: Make this adjustable - issue#2
     peakrange = np.arange(1, 100)
     rawpeaks = np.array(sc_sg.find_peaks_cwt(data, peakrange))
     rawpeaks = rawpeaks[rawpeaks > rawpeaks.min()]  # Remove leftmost peak (false positive at input)
     threshpeaks = rawpeaks[data[rawpeaks] > thresh]
     return np.array([domain[threshpeaks], data[threshpeaks]]).T

def checkerboard_matrix_filtering(similarity_matrix, kernel_width, peak_range):
    """
    Moving the checkerboard matrix over the main diagonal of the similarity matrix one sample at a time.

    :param similarity_matrix:
    :param peak_range: the number of samples in which the peak detection algorithms finds a peak (TODO:clarify)
    :param kernel_width: the size of one quarter of the checkerboard matrix
    :return: peaks and convolution values
    """

    checkerboard_matrix = get_checkerboard_matrix(kernel_width)

    # The values calculated in this step are starting from the 'kernel_width' position and ending
    # at length - kernel_width
    d = []
    for i in range(0, similarity_matrix.shape[0] - 2 * kernel_width):
        base = similarity_matrix[i:i + kernel_width * 2, i:i + kernel_width * 2]
        d.append(np.sum(np.multiply(base, checkerboard_matrix)))

    # The missing values from 0 to kernel_width are calculated here
    top_left_d = []
    for i in range(0, kernel_width):
        base = similarity_matrix[0:i + kernel_width, 0:i + kernel_width]
        top_left_d.append(np.sum(np.multiply(base, checkerboard_matrix[kernel_width - i:, kernel_width - i:])))

    # The missing kernel_width values at the bottom right are set to 0
    convolution_values = top_left_d + d + [0 for i in range(0, kernel_width)]

    peaks = find_peaks_cwt(convolution_values, np.arange(1, peak_range))
    peaks = [0] + peaks + [len(convolution_values)-1]
    return peaks, convolution_values

def call_boundary_peaks(Pb):
    n = len(Pb)
    peaks = find_peaks_cwt(
        Pb, 
        np.array([0.5]), 
        wavelet=None, 
        max_distances=None, 
        gap_thresh=None, 
        min_length=None, 
        min_snr=1, 
        noise_perc=10)

    regions = []
    for i in peaks:
        top = Pb[i]
        x = [l for l in range(1, min(i, 5)) if Pb[i-l] > 0.4*top]
        y = [r for r in range(1, min(n-i, 5)) if Pb[i+r] > 0.4*top]
        start = (i - x[-1]) if len(x) else i
        end = (i + y[-1]) if len(y) else i
        regions.append([start-1, end+1, top])

    if len(regions):
        starts, ends, values = zip(*regions)
    else:
        starts, ends, values = [], [], []

    return np.array(starts), np.array(ends), np.array(values)

def dih_plotter3(dirname,savename,numplot):
    inlist = dih_tablereader(dirname)
    plotlist = inlist[0:numplot]
    colors = iter(cm.rainbow(np.linspace(0,1,len(plotlist)))) #creates color table
    for memberlist in plotlist:
        x = memberlist[0] #x coordinate data
        y = memberlist[1] #y coordinate data
        ysmooth = dih_boxcar(y)
        xshort = x[0:len(ysmooth)]
        peaklist =signal.find_peaks_cwt(ysmooth, np.arange(1,30))#continuous wavelet transformation
        plt.plot(xshort,ysmooth,color = next(colors))
        for num in peaklist:
            plt.plot(xshort[num],ysmooth[num],'gD')#places markers on peaks
        peak = max(ysmooth)
        peaklist2 = [i for i, j in enumerate(ysmooth) if j == peak]
        for num in peaklist2:
            plt.plot(xshort[num],ysmooth[num],'rD')

#finish up plot characteristics
    plt.title('Super Most Awesome Graph!')
    plt.ylabel('Flux')
    plt.xlabel('Time')       
    pylab.ylim([-5,5])
    pylab.xlim([0,6.3])
    plt.savefig(savename)#saves postscript file
    return plotlist

 def curPeakdetect(y, x):
     if ds["type"] == "cwt":
         _maxidx = find_peaks_cwt(y, ds["widths"])
         _max = zip(x[_maxidx], y[_maxidx])
     else:
         _max, _ = peakdetect(y, x, lookahead=ds["lookahead"], delta=ds["delta"])
     return _max

def calculatepeak(step, jump, dwell):
	totalpeak = np.array([])
	totalinterval = np.array([])
	totalfilterpeak = np.array([])
	totalfilterinterval = np.array([])
	for window in step:
		if window/jump % 2 == 0:
			window += jump
		ma = convertma(dwell, window, jump)
		interval = window / jump
		#Calculate the moving average of the result
		std = movstd(ma, interval)
		#Calculate the moving standard deviation of the result
		oripeak = signal.find_peaks_cwt(std, np.arange(interval * percentage, interval, interval*0.05))
		oristdamp = std[oripeak]
		#Obtain the standard deviation peak and the amplitude of the peak

		if len(oripeak) == 0:
			pass
		else:
			oripopendiffpeak = calpopendiffpeak(dwell, oripeak)
			#Calculate the difference in Popen from the raw data using the time of the peak
			oripopendiffstdamp = oristdamp / np.sqrt(1.0/12)
			#Calculate the difference in Popen based on the value of the standard deviation peak
			filterpeak, oripopendiffstdamp, oripopendiffpeak = filterpeaktime(dwell, std, oripeak, oripopendiffstdamp, oripopendiffpeak)

			totalpeak = np.append(totalpeak, oripeak)
			totalinterval = np.append(totalinterval, [interval] * len(oripeak))
			totalfilterpeak = np.append(totalfilterpeak, filterpeak)
			totalfilterinterval = np.append(totalfilterinterval, [interval] * len(filterpeak))
			for iprint in range(len(filterpeak)):
				print 'peak', filterpeak[iprint], 'Popen(avg)', oripopendiffpeak[iprint], 'Popen(std)', oripopendiffstdamp[iprint]

	np.savez(pathfilename+str(step[0])+'_'+str(step[-1]), totalpeak=totalpeak, totalinterval=totalinterval, totalfilterpeak=totalfilterpeak, totalfilterinterval=totalfilterinterval)

def determine_step_process(dwell, percentage, start, end, jump):
	attempt = np.array([])
	firstpeak = np.array([])
	lastpeak = np.array([])
	suminterval = np.array([])
	for window in np.arange(start, (end+jump*2), jump*2):
		#print 'testing interval:', interval,
		#window must be an odd number
		if window/jump % 2 == 0:
			window += jump

		ma = convertma(dwell, window, jump)
		interval = window / jump
		std = movstd(ma, interval)
		oripeak = signal.find_peaks_cwt(std, np.arange(interval * percentage, interval, interval*0.05))
		if oripeak != []:
			peaklength = np.append(oripeak, np.cumsum(dwell)[-1]/jump-1) - np.append(0, oripeak)
			if np.amin(peaklength) > interval:
				attempt = np.append(attempt, 1)
			else:
				attempt = np.append(attempt, 0)
			firstpeak = np.append(firstpeak, oripeak[0])
			lastpeak = np.append(lastpeak, oripeak[-1])
			suminterval = np.append(suminterval, interval)

	np.savez(pathfilename+'_'+str(start)+'_'+str(end), attempt = attempt, firstpeak = firstpeak, lastpeak = lastpeak, suminterval = suminterval)

def sort_spikes(sweep):
    filtered = np.array(sweep) - medfilt(sweep, 51)
    pks = find_peaks_cwt(filtered, np.arange(9, 20))
    offsets = np.zeros_like(pks)
    datamtrx = np.zeros(shape=(80, len(pks[3:])))

    for i, pk in enumerate(pks[3:]):
        offset = np.argmax(filtered[pk - 40 : pk + 40]) - 40
        datamtrx[:, i] = filtered[pk - 40 + offset : pk + 40 + offset]
        offsets[i + 3] = offset
        # plb.plot(filtered[pk-40+offset:pk+20+offset],color = 'k', alpha = 0.1)

    from scipy.linalg import svd

    U, s, Vt = svd(datamtrx, full_matrices=False)
    V = Vt.T

    ind = np.argsort(s)[::-1]
    U = U[:, ind]
    s = s[ind]
    V = V[:, ind]

    features = V
    es, idx = kmeans2(features[:, 0:4], 4, iter=50)
    colors = [([1, 0, 1], [1, 0, 0], [0, 0, 1], [0, 1, 1])[i] for i in idx]

def autocorrelation(x, y):

    # normalise so range is 2 - no idea if this is the right thing to do...
    y = 2*y/(max(y)-min(y))
    y = y-np.median(y)

    # calculate autocorrelation fn
    lags, acf, lines, axis = pb.acorr(y, maxlags = len(y)/2.)

    # halve acf and find peaks
    acf = acf[len(acf)/2.:]
    pks = find_peaks_cwt(acf, np.arange(10, 20)) # play with these params,

    #they define peak widths
    peaks = pks[1:] # lose the first peak
    period =  peaks[0]*cadence
    print 'acf period = ', period

    pl.clf()
    pl.subplot(2,1,1)
    pl.plot(x[:4000], y[:4000], 'k.')
    pl.title('Period = %s' %period)
    pl.subplot(2,1,2)
    pl.plot(np.arange(5000)*cadence, acf[:5000])
#     pl.plot(np.arange(len(acf))*cadence, acf)
    [pl.axvline(peak*cadence, linestyle = '--', color = 'r') for peak in peaks]
    pl.xlim(0, 5000*cadence)
    pl.savefig('/Users/angusr/Python/george/acf/%sacf' %int(KID))
#     np.savetxt('/Users/angusr/Python/george/acf/%sacf_per.txt'%int(KID), period)

    return period

def detect_lines(w_arr, f_arr, sigma=3, bsigma=None, niter=5, mask=None,
                 kern=default_kernal, center=False):
    """Detect lines goes through a 1-D spectra and detect peaks

      w_arr--xaxis array (pixels, wavelength, etc)
      f_arr--yaxis array (flux, counts, etc)
      sigma--Threshold for detecting sources
      bsigma--Threshold for determining background statistics
      niter--iterations to determine background
      center--return centroids and not pixels
      mask--Pixels not to use
    """
    # set up the variables
    if bsigma is None:
        bsigma = sigma

    if mask:
        f_arr = f_arr[mask]
        w_arr = w_arr[mask]

    # find all peaks
    xp = signal.find_peaks_cwt(f_arr, np.array([sigma]))
    xp = np.array(xp)
  
    # set the output values
    if center:
        xdiff = int(0.5 * len(kern) + 1)
        xp = xp * 1.0
        for i in range(len(xp)):
            xp[i] = mcentroid(w_arr, f_arr, kern=kern, xdiff=xdiff, xc=w_arr[xp[i]])

    return xp

    def wavelet_peaks(vector, *args, **kwargs):
        """
        This is the function that will be mapped by multiprocess. This is a wrapper around the scipy function.
        It uses a parameter - wavelet_widths that is configured outside this function.

        Parameters
        ----------
        vector : 1D numpy array
            Feature vector containing peaks

        Returns
        -------
        peak_indices : list
            List of indices of peaks within the prescribed peak widths
        """
        try:
            peak_width_bounds = kwargs.get('peak_widths')
            kwargs.pop('peak_widths')
            peak_width_step = kwargs.get('peak_step', 20)
            kwargs.pop('peak_step')
            # The below numpy array is used to configure the returned function wpeaks
            wavelet_widths = np.linspace(peak_width_bounds[0], peak_width_bounds[1], peak_width_step)

            peak_indices = find_peaks_cwt(np.abs(vector), wavelet_widths, **kwargs)

            return peak_indices

        except KeyError:
            warn('Error: Please specify "peak_widths" kwarg to use this method')

def mean_hr_bpm(filename):
    """module to take user input for time scale, and analyze ECG input in \
    that time scale

    :param filename: the name of a file located in the /test_data folder \
    entered as a string

    :returns heartrate: heartrate during a specific period as a float
    :raises IOError: raised if user tries to input value not accepted by \
    program
    :raises ValueError: raised if the generally accepted values fall outside \
    of the signal time range
    """

#   time_input = input("Please input time (10 sec or 20 sec): ")
    time_input = "10 sec"
    time_vector = extract_time_data(filename)
    if np.max(time_vector) >= float(time_input[:-4]):
        if str(time_input) == "10" + " sec":
            ind = np.where(time_vector == 10)[0]
            df = import_data(filename)
            values = df.values
            trimmed = values[np.arange(0, ind), 1]
            trim_norm = trimmed - np.mean(trimmed)
            template = pd.read_csv("test_data/template.csv", header=None)
            norm_template = extract_template_data(template)
            corr = np.correlate(norm_template, trim_norm, mode="full")
            peaks = signal.find_peaks_cwt(corr, np.arange(1, 300))
            heartrate = len(peaks) / (10/60)
        elif str(time_input) == "20" + " sec":
            ind = np.where(time_vector == 20)[0]
            df = import_data(filename)
            values = df.values
            trimmed = values[np.arange(0, ind), 1]
            trim_norm = trimmed - np.mean(trimmed)
            template = pd.read_csv("test_data/template.csv", header=None)
            norm_template = extract_template_data(template)
            corr = np.correlate(norm_template, trim_norm, mode="full")
            peaks = signal.find_peaks_cwt(corr, np.arange(1, 300))
            heartrate = len(peaks) / (20/60)
        else:
            raise IOError("Invalid input. Try Again (Make sure to include "
                          "sec)")
    else:
        raise ValueError("Attempted input outside signal range")
    return heartrate

def FindPeaksInSpectrum (S):
  "Find peaks using scipy"

  Peaks = []

  peakind = signal.find_peaks_cwt(S, numpy.arange(1, 10), min_length=16)

  return peakind